membrane proteins play 
 a key role in the structure and function 
 of cells . most membrane proteins have transmembrane ( tm ) domains formed 
 by bundles of -helices that span the lipid bilayer . the dynamics 
 of tm domains are responsible for conformational changes that underlie , 
 e.g. , receptor activation and solute transport by membrane proteins . thus , understanding the underlying free energy landscapes of tm helix 
 interactions provides the key to mechanisms of membrane protein function . a relatively simple and well characterized example of such interactions 
 is provided by tm helix dimers such as those found in receptor tyrosine 
 kinases ( rtks ) . the receptor tyrosine kinases are a major family 
 of membrane receptors 
 controlling many cellular activities . among 
 the rtks , epidermal growth factor receptors ( egfr / erbb1/her1 ) aberrant erbb activity is associated with a number of illnesses , including 
 cancers and neurodegenerative diseases , and therefore 
 erbb receptors are important drug targets . understanding 
 the molecular mechanisms of signal transduction by the egfr is of 
 fundamental importance in providing a paradigm for other rtks and will facilitate development of a range of therapeutic interventions . there has been extensive structural characterization of the 
 extracellular 
 ectodomain and the intracellular kinase domain of the egfr . binding of the ligand ( egf ) 
 to the ectodomain of the receptor induces a conformational change 
 that is transmitted through an allosteric mechanism to the kinase 
 domain . formation of an asymmetric dimer 
 of the kinase domains enables tyrosine phosphorylation , which in turn 
 activates the downstream signal transduction cascade in the cytoplasm . although cryo - electron microscopy has yielded 
 images of the intact egfr , in the absence 
 of a high resolution structure of the complete receptor , including 
 the transmembrane ( tm ) domain , our understanding of its activation 
 mechanism remains incomplete ( figure 1 ) . the tm helix dimer is 
 believed to play a key role in the activation mechanism . mutagenesis experiments have also highlighted 
 the role of the juxtamembrane domain ( jm ) , immediately c - terminal 
 to the tm domain , in regulating the allosteric transition between 
 inactive and active states of egf - bound receptor dimers and in their interactions with lipids . truncation of the ectodomain activates intracellular 
 phosphorylation , demonstrating that the 
 ectodomain has an autoinhibitory function , and that the tm domain 
 can independently activate the receptor . it is suggested that the active state is characterized by dimerization 
 contacts in the n - terminal region of the tm helices ( figure 1c , motif highlighted in red ) 
 along with an antiparallel interaction between jm segments . the nature of the inactive state(s ) remains 
 somewhat less clear but is thought to involve a more c - terminal interaction 
 of the tm helices ( as recently discussed in the literature ) . ( a ) schematic diagram of the structure of a 
 ligand - bound egfr dimer . the egfr monomers are in cyan and ice blue , the egf ligands in red . the model was built using the structures of the ligand - bound ectodomain 
 dimer ( pdb i d : 3njp ) , of the tm helix dimer ( pdb i d : 2m20 ) , and of the asymmetric kinase domain 
 dimer ( pdb i d : 3gop ) . ( c ) sequence 
 of the tm domain and schematic representation of the secondary structure . cationic residues are in purple ; residues from n - terminal amd c - terminal 
 small - x3-small motifs are in red and green , respectively . in this context it is important 
 to elucidate modes of interaction 
 of the helices within the dimeric tm domain . the conformational 
 dynamics of these interactions may be explored through computer simulations , 
 which have been used to provide detailed descriptions of different 
 ensembles of conformations visited by tm helix dimer systems . however , due to the limiting time scale of atomistic molecular dynamics 
 ( md ) simulations only partial description of egfr dynamics have so 
 far been achieved . consequently , models which employ 
 a simplified ( e.g. , coarse - grained : cg ) molecular representation are 
 valuable , as they may be used to investigate 
 tm helix dimerization events beyond the accessible time scale of atomistic 
 simulations . cg simulations have been used to obtain one - dimensional 
 projections 
 of the free energy landscape of tm helix dimerization for glycophorin 
 a , for erbb , and 
 for epha1 receptors . however , in order 
 to link the conformational motions of tm helices to possible activation 
 mechanisms , an accurate and exhaustive ( i.e. , multidimensional ) description 
 of the underlying free energy landscape is essential . this remains 
 beyond the accessible time scale of simple cg simulations , even with 
 the availability of advanced supercomputing resources . these methods accelerate the simulation of the process of interest , 
 allowing capture of long time scale motions at affordable computational 
 cost . here we present an innovative computational protocol that 
 by combining 
 cg simulations with well - tempered metadynamics ( cg - metad ) greatly extends the time scales typically reached by 
 simulations . in particular , cg models reduce the dimensionality of 
 the system under investigation , while 
 metadynamics enhances the phase space exploration , thus allowing practical 
 investigation of long time scale and large - scale motions of biomolecules . pioneering applications of metadynamics to cg models have indeed 
 shown great potential . in the present study we employ 
 380 s of cg - metad simulations to fully inspect both 
 association and dissociation events of the egfr tm 
 helices . this provides an accurate and comprehensive description of 
 the associated free energy landscape , which allows identifying the 
 low free - energy basins that correspond to distinctive tm dimer conformations . the reversible pathways of transition between these dimeric states 
 reveal key functional motions of the tm helices to reach the active 
 dimer state , thus unveiling previously unknown features of the activation 
 mechanism . this information is of paramount relevance to develop an 
 exogenous control of the receptor activity . the conformational stability 
 of the identified low energy dimer states has been confirmed by over 
 3 s of atomistic md simulations , thus providing structural 
 details that explain and complement previous studies . this study demonstrates 
 that cg - metad provides a powerful method to investigate large - scale 
 motions in the tm domains of membrane proteins , thus representing 
 an important advance in simulation of protein conformational dynamics 
 in a membrane environment . the starting structure of the 
 tm domain dimer was a homology model based on the nmr structure of 
 the erbb1/erbb2 heterodimer ( pdb i d : 2ks1 ) , generated 
 using modeler . coarse - graining used an 
 in - house variant of the martini force field . in addition to the standard bonded interactions , harmonic restraints 
 are generally required to conserve secondary structure elements . therefore , 
 in the tm models , harmonic restraints were applied between backbone 
 particles separated by three residues in the -helical region 
 ( from ile to his ) to mimic hydrogen - bonds , 
 while the intrinsic flexibility of the n - loop and c - loop regions were 
 maintained with no additional restraints . in the tm+jm models ( based 
 on the 2m20 nmr structure ) , an elastic network model acts between 
 pairs of backbone particles separated by less than 0.7 nm . the cg - metad simulation was performed 
 with the well - tempered version of metadynamics , implemented in plumed 1.3 . the estimate f(s , t ) of the free - energy surface f(s ) at time t was determined as1where v(s , t ) is the bias potential added to the 
 system at 
 time t as a function of a set of selected cvs , s. t is the difference between the fictitious temperatures ( t + t ) of the selected cvs and the temperature t of the system . the bias potential v(s , t ) is made up by the 
 sum of the gaussians deposited so far . in our case the main 
 simulation was performed at 323 k ( i.e. , the temperature at which 
 the martini force field was parametrized ) with a time step of 40 fs . 
 the lateral separation ( l ) between the center of 
 mass of the helices was chosen as single active collective variable . 
 gaussians of height 0.05 kj / mol and width 0.05 nm were initially used 
 and deposited every 5000 steps with a bias factor ( t + t)/t of 10 . this 
 means that the gaussian height was gradually decreased along the 380 
 s of cg - metad simulation on the basis of the adaptive bias 
 with a t of 2907 k ( supporting information figure s1 ) . additional cg - metad simulations 
 were performed using different 
 parameters : the width of the gaussians was reduced to 0.025 nm in 
 one simulation or alternatively a soft wall potential was used to 
 restrain the lateral helix separation l within 4.5 
 nm ( supporting information figure s2 ) . the results obtained from the different simulations led to a coherent 
 description of the free energy landscape ( supporting information figure s3 ) , demonstrating the robustness of the 
 approach . in contrast we observed both tempering and convergence issues 
 with higher deposition rates ( either with higher gaussians or higher 
 deposition frequencies ) . unless stated in the text , cg - metad results 
 refer to the main simulation . a series of standard cg simulations 
 were performed using gromacs 4.5.5 ( i.e. , 
 185 s of cg - md simulations , see supporting information table s1 ) . a leapfrog algorithm was 
 used to integrate newton s equations of motion , with a time 
 step of 40 fs . standard cg - md and cg - metad simulations were performed 
 in npt and nvt ensembles , respectively . the temperatures of the individual monomers , the solvent and the lipid 
 bilayer were separately maintained at 323 k using the bussi thermostat 
 with relaxation times of 1 ps . in standard 
 cg - md simulations , semi - isotropic pressure coupling was used . the 
 berendsen barostat was employed to keep the average pressure at 1 
 bar , with time constants of 10 ps and 
 compressibilities of 3  10 bar . lennard - jones interactions were shifted to zero between 0.9 and 
 1.2 nm , while electrostatic interactions were shifted to zero between 
 0 and 1.2 nm using a relative dielectric constant of 20 . neighbor 
 lists were updated every 10 steps for the calculation of nonbonding 
 interactions , using a 1.3 nm cutoff . the populations of lowest energy 
 identified by cg - metad were converted to atomistic level using a fragment - based 
 approach . a total of over 3 s 
 of standard atomistic md simulations ( at - md ) were performed using 
 the amber force field ff14sb and lipid14 for protein and lipids , respectively . all the simulations 
 were performed with the pmemd module of the amber md code . short - range van der waals interactions were 
 switched to zero at a cutoff distance of 1 nm . the long - range electrostatic 
 interactions were computed by means of the particle mesh ewald ( pme ) 
 method using a real - space cutoff of 1 
 nm and a 0.1 nm grid spacing in periodic boundary conditions . the 
 shake algorithm was applied to constraint bonds involving hydrogen 
 atoms , and thus an integration time step of 2 fs could be used . harmonic 
 constraints were applied to the protein and gradually released along 
 the equilibration process . the temperatures of the individual monomers , 
 the solvent and the lipid bilayer were maintained at 323 k using the 
 langevin thermostat with collision frequency of 3 ps . production runs were performed in the npt ensemble 
 at 1 atm under anisotropic pressure scaling conditions using a monte 
 carlo barostat and pressure relaxation time of 2 ps . analysis of the simulations 
 was carried out using plumed 1.3 , gromacs 
 4.5.5 and locally written scripts . the 
 probability distribution function of various cvs was reweighted for 
 free energy calculations . the cvs used 
 were : l , the lateral distance between the centers 
 of mass of the tm helices ;  , the crossing angle between the 
 tm helices ; d , calculated as the difference 
 in interhelix distance at the n - termini and the c - termini ; and || , 
 which describes the rotation of the helices along their long axis 
 relative to one another . more detailed definitions are provided in 
 the paper when these cvs are used . however , the crossing angle cv 
 (  ) benefits here from a more complete description . the positions 
 of the backbone particles of the tm helical region were used to determine , 
 by singular value decomposition ( svd ) , the eigenvectors h1 and h2 associated with the helix axes . the absolute value of crossing angle 
 was determined using the definition || = asin(|h1  h2|/(|h1||h2| ) ) . one might also consider atan2(|h1  h2|,h1h2 ) that leads to identical || 
 values . in order to compute the sign of , h1 was first aligned with the positive x - axis , with the same rotation applied to h2 , which provided two new helix vectors h1 and h2. then , the 
 cross product h1 h2 was aligned with the x - axis , and the same rotation 
 applied to the two helix vectors provided h1 and h2. all rotations were done using 
 the origin of the cartesian coordinate system as center . the sign 
 of  was determined by considering two properties : ( i ) the sign 
 of the x - component of the geometric center of the cross product h1  h2 , ( ii ) the 
 comparison of the x - component of the geometric centers of h1 and h2. negative and positive 
 signs of  characterized right- and left - crossings , respectively . all minimum free 
 energy paths ( mfeps ) presented in the paper were calculated in one 
 of the two dimensions of the free energy surfaces obtained by cg - metad 
 calculations . the paths were obtained by scanning of the lowest free 
 energy points of the surface along the ( cv ) dimension of interest . for the determination of the free energy profiles , partition functions 
 were calculated at each value of the cv of interest . thus , contrary 
 to mfeps , all the free energy points obtained in the second dimension 
 of the fes contributed to the free energy profiles . the comparison 
 showed that both are really similar , although the mfeps allowed a 
 more accurate detection of the free energy basins of the fes . therefore , 
 we chose to present our results on the basis of the mfeps . importantly 
 the relative free energies of the basins identified from the different 
 mfeps were similar for each dimer population ( supporting information figures s3 and s4 ) . this guarantees 
 both the validity of the free energy calculations and the suitability 
 of the different cvs to describe the multidimensional conformational 
 landscape of the egfr tm domain ( supporting information table s2 ) . a starting structure 
 for the simulations 
 of the egfr tm helix dimer was based on the nmr structure of the erbb1/erbb2 
 heterodimer ( pdb i d : 2ks1 ) . this structure includes flexible n - terminal ( extracellular ) and 
 c - terminal ( intracellular ) loops , which we will refer to as the n - loop 
 and c - loop , respectively ( figure 1 ) . therefore , in the cg simulations , harmonic restraints 
 were applied ( see methods above ) to maintain 
 the secondary structure of the core -helical region ( from ile to his ) while retaining the intrinsic flexibility 
 of the n - loop and c - loop regions . the n - loop contains residue k618 
 that is thought to interact with the glycolipid gm3 in cell membranes . the c - terminus of the tm helix plus the c - loop contains six basic 
 residues ( r645 to r653 ) that can interact with anionic lipids ( e.g. , 
 pip2 ) of the intracellular membrane leaflet . the tm helix dimer was embedded 
 in a phospholipid bilayer . to describe 
 the motion of one helix with respect to the other we used the lateral 
 separation between the centers of mass of the helices ( l ) as a collective variable ( cv ) . the sampling of tm helix dimerization 
 took a total of 0.4 ms of cg - metad simulation , during which 
 65 interconversion events between the dimeric and monomeric states 
 of the tm helices were observed . thus , the diffusive exploration of 
 the cv space in our cg - metad simulation yielded a well - characterized 
 and converged free energy surface ( figure 2 and supporting information figure s1 ) . ( a ) interconversion of the tm domain between three states : 
 monomeric 
 state ( lateral helix separation , l > 3 nm ) , predimeric 
 ( 1.5 nm < l < 3 nm ) , and dimeric ( l < 1.5 nm ) . ( b ) top - view of the free energy surface ( fes ) as a 
 function of l and of the crossing angle . 
 isoenergy lines are plotted every 5 kj / mol on the fes . conformations corresponding to the free energy 
 basins are drawn with gray tubes for the helix backbone and colored 
 spheres for residues engaged in interhelical contacts . projection 
 of the fes in the l and  dimensions shows 
 isoenergy lines every 5 kj / mol ( middle panel ) . one dimensional free 
 energy profile as a function of the active cv l ( lower 
 panel ) . from the free energy surface 
 ( fes ) thus , the free energy of tm helix dimerization 
 was calculated from the fes as the difference between the dimeric 
 and monomeric states , yielding an estimate of 38  3 
 kj / mol ( error calculated as the standard deviation from the weighted 
 average value of the tm helix dimerization free energy obtained from 
 the last part of the simulation where the calculation is converged ; 
 see also supporting information figure 
 s1c ) . this is in reasonable agreement with an experimental estimate 
 of the dimerization free energy of 31 kj / mol , derived from 
 a kd of 2 m in fret experiments 
 on egfr tm helices in ldao micelles . we 
 do note that these values are larger than some other computational 
 and experimental estimates , of 26 kj / mol and 11 kj / mol respectively . however , considering the sensitivity of the dimerization process to 
 experimental conditions ( e.g. , lipid vs detergent ) , the approximations 
 implicit in the cg force fields used in simulations , and especially 
 the different construct sequences of the egfr tm domain used in the 
 various studies ( supporting information figure s5 ) , it is perhaps unsurprising that we observe differences . indeed , we note that a wide range of values ( from 30 kj / mol 
 to 60 kj / mol ) has been estimated for the dimerization free 
 energy of the canonical glycophorin a helix dimer on the basis of 
 atomistic and cg simulations . using a reweighting protocol 
 we computed the boltzmann distribution of cvs different from those 
 biased in the original calculation . this 
 allows for a more detailed inspection and characterization of the 
 motions observed during the simulation . in this manner we reconstructed 
 the fes as a function of both the original cv ( the interhelix distance l ) and of a cv corresponding to the interhelix crossing 
 angle (  ) the latter is a frequently used structural descriptor 
 of the packing interactions within tm helix dimers . the lower minimum 
 corresponds to right - handed (  = 30 ) tm helix 
 dimer conformations , while the second minimum is 9 kj / mol higher and 
 corresponds to a left - handed dimer (  =  + 12 ; figures 2bc and 3a ) . within the right - handed energy basin of the 
 tm dimer , we refer to the most populated ensemble of conformations 
 as nter+ ( figure 3d ) since it shows the two helices interacting through 
 residues of the n - terminal motif along with additional contacts at 
 the helix core ( gxxgaxxllxxxa ) . we will refer to the left - handed dimer as the lzip population ( figure 3d ) since the helices form leucine zipper - like interactions 
 ( vxxlxxxl ) . the nter+ state is 11  2 kj / mol more stable than the lzip state . schematic definition and associated free energy profiles 
 of the 
 cvs : ( a )  , ( b ) d , and ( c ) || . is the crossing angle between helices ; d cv is calculated as the difference in interhelix distance at the 
 n - termini and the c - termini ( defined as the centers of mass of tgxxga and axxxg motifs , respectively ) ; || describes the rotation of 
 the tm helices about their long axes relative to one another ( see 
 also supporting information figure s4 ) . each free energy profile was calculated as the minimum free energy 
 path ( mfep ) located at low l values of the associated 
 fes ( see methods ) . ( d ) structure of the dimer 
 populations identified in the basins of the free energy landscape 
 of egfr tm domain . colored spheres represent residues mainly engaged 
 in interhelical contacts , while gray tubes show helix backbones . although the helix crossing angle is a metric commonly 
 used to 
 characterize tm helix dimers , it may 
 not allow us to readily distinguish subtly different dimer conformations 
 within the same macro - ensembles ( i.e. , within the right - handed and 
 left - handed energy basins ) . therefore , we undertook a series of local 
 analyses of all dimer conformations . the sequence 
 of the egfr tm domain contains a tandem small - x3-small 
 sequence motif toward the n - terminus ( tgxxga ) and a single small - x3-small motif ( axxxg ) toward the c - terminus ( figure 1 ) . such small - x3-small motifs play a key role in tm helix dimerization in 
 a number of membrane proteins . interconversion 
 between different tm helix packing modes has been implicated in activation 
 of the egfr and related receptors ( see the literature for a recent discussion ) . on this basis , we defined a new 
 cv , d , which corresponds to the difference 
 in interhelix distances between these two motifs , such that a negative 
 value of d indicates an n - terminal helix / helix 
 interaction and a positive value indicates a c - terminal interaction 
 ( figure 3b ) . the minimum 
 free energy path ( mfep ) along the d cv , corresponding 
 to the lowest energy points computed in the d range from 1 to + 1 nm , reveals the presence of four energy 
 minima ( figure 3b ) . the deepest energy minimum ( at d = 0.1 
 nm ) corresponds to the previously described nter+ population . an additional low free - energy minimum , named cter+ , shows the two helices interacting through both the 
 c - terminal small - x3-small motif and some of the core hydrophobic 
 residues ( axxxlxxxaxxxg ) . we note that the nter+ population is 
 slightly more stable than cter+ ( 1.4 kj / mol ) , 
 as suggested by previous experimental studies . the other two free - energy minima , nter and cter , are 6 to 7.5 kj / mol higher than nter+ and correspond to helix dimers in contact primarily 
 through their n - terminal and c - terminal small - x3-small 
 motifs , respectively . these interactions are similar to those found 
 in the nter+ and cter+ conformations , 
 with the lack of additional contacts at the core of the helices explaining 
 the higher free energy values of these populations ( supporting information figure s4a  d ) . from the 
 analyses presented so far however , the left - handed lzip dimer is absent from the fes computed as a function 
 of d ( figure 3b ) . accordingly , we defined an additional cv , || , 
 which is able to further characterize specific conformational differences 
 among the diverse dimer structures . this cv describes the rotation 
 of the tm helices about their long axes relative to one another . the 
 mfep computed as a function of || reveals three out 
 of the five dimer populations previously described ( figure 3c ) . in this representation cg - metad 
 simulations have allowed us to explore the large - scale 
 motions of the tm domain , revealing the existence of five dimeric 
 states ( figure 3d , supporting information table s2 ) . we have used 
 each of these dimeric state models as the starting point for extended 
 atomistic simulations , in order to explore in more detail specific 
 inter - residue interactions between the two helices . thus , each dimeric 
 model of the five lowest energy minima identified by cg - metad was 
 used to initiate a 0.5 s of unrestrained atomistic ( at ) md 
 simulations ( yielding a total simulation time of 2.5 s ) . we 
 found that the dimer conformations were stable during the at - md simulations , 
 with , e.g. , no significant change in helix crossing angles over the 
 course of the 0.5 s simulations ( figure 4a ; also see supporting information figures s6 ) . examination of the structures at the 
 end of the atomistic simulations confirmed that the hydrophobic residues 
 of the tm domain form interactions with the lipid aliphatic chains 
 within the bilayer core , while interhelix contacts are formed by the 
 interaction motifs identified through the cg - metad simulation . we 
 found that a number of h - bonding interhelix interactions involving 
 the n - ter and c - ter loop residues also contribute to stabilize the 
 dimeric states , in addition to electrostatic interactions between 
 cationic residues and lipid headgroups ( figure 4b ; also see supporting information si text ) . atomistic md simulations initiated from each 
 dimer conformation 
 of the five lowest energy minima identified by cg - metad . ( a ) crossing 
 angles as a function of time for the five at - md simulations . ( b ) atomistic 
 models selected from the five egfr dimer populations obtained from 
 at - md simulations . helix - helix interactions at loop regions are shown 
 as insets . h - bonds are displayed as black dashed lines . in the central 
 panel a superimposition of the nter+ , nter and lzip conformations is shown , highlighting the 
 intermediate role played by nter in the lzip - to - nter+ transition . tm domains are represented 
 as gray cartoons , zwitterionic headgroups of phospholipids as yellow 
 spheres , while waters and lipid tails are omitted for clarity . examination of the cg - metad derived 
 free energy landscape reveals 
 that in the unbound state ( i.e. , monomeric tm helices at l > 3 nm ) the surface is not flat along the  dimension ( figure 2c ) . the shallow free - energy 
 dependence on the crossing angle , with favorable free energies at 
   18 and   + 18 , 
 corresponds to the ( monomeric ) helices adopting a tilted orientation , 
 with a lowest free - energy tilt angle of 20 relative 
 to the bilayer normal ( supporting information figure s7 ) , close to values reported for the tilt of synthetic peptides 
 in model lipid bilayers . the bias added to the system 
 during the cg - metad calculation allows acceleration of sampling , and 
 the reconstruction at the end of the simulation of the thermodynamics 
 of the event under study ( i.e. , the free - energy landscape and free - energy 
 minima identification ) . the characterization of the two - dimensional 
 fes ( figure 2b ) also 
 therefore provides information on the possible pathways connecting 
 the monomeric state to the different dimeric structures . we divided 
 the fes into two regions separated by the free energy barrier observed 
 at a crossing angle of 0 in order to calculate mfeps for dimerization 
 ( i.e. , from high to low l values ) . we observed similar 
 features for the calculated pathways in the two regions ( figure 5a ) . at high l values , the mfeps linearly follow the free energy basins 
 at   18 and   + 18 
 discussed above , which represents the tilted orientations of tm helices 
 in their monomeric state . in contrast at intermediate l values , ( 1.5 nm < l < 3 nm ) corresponding 
 to predimeric states , the mfeps show a nonlinear pattern with an increase 
 of the magnitude of the crossing angle between the tm helices when 
 they first contact each other . this is followed by a decrease in the 
 magnitude of the crossing angle when the helices pack together more 
 closely ( figure 5a ) . 
 subsequently , the mfeps follow the curvature of the free energy wells 
 observed at low l values . ( a ) minimum free energy 
 pathways ( mfeps ) between monomeric and 
 dimeric states in the regions of the fes with  above or below 
 the zero crossing angle . mfeps are drawn as black lines between monomeric 
 states with high l ( > 5 nm ) and the free energy 
 minima 
 of the dimeric states , while mfeps between these latter points and 
 the dimeric states with lower l values are shown 
 as gray lines . free energy profiles of these pathways are in the supporting information ( figure s1f ) . ( b , c ) regions 
 of the fes explored during 5 s of standard cg - md simulation 
 starting from two randomly chosen initial monomeric states : ( b ) l = 3 nm and  = 28 ; and ( c ) l = 3 nm and  = 7. each frame of the 
 standard simulations is represented by a black dot . these findings prompted us to compare the dimerization pathways 
 predicted by cg - metad to those observed through unbiased cg simulations . we therefore performed three standard cg - md simulations , each of 5 
 s duration starting from monomeric states ( l  3 nm ) with different  values . these simulations 
 displayed dimerization of the tm helices with features similar to 
 those predicted by the mfeps . however , comparison of the cg - md simulations 
 showed a range of trajectories on the fes , highlighting the stochastic 
 behavior of individual dimerization events . dimerization of 
 the tm helices occurred early after 50 
 ns in the first simulation , following the mfep calculated in the region 
 of negative  values of the fes ( figure 5b ) . in the second simulation a wider exploration 
 of the fes at high l values the longer diffusion time 
 of the helices in the monomer state resulted in a well - balanced representation 
 of the two low energy populations , with both right - handed and left - handed 
 monomeric states sampled . in this simulation the two tm domains first 
 interacted with each other through the extracellular n - loops . these 
 predimeric states can readily convert from left- to right - handed forms 
 in agreement with the quasi - linear isoenergy lines observed along 
 the  axis in this region of the fes ( 1.5 nm subsequently , more stable interactions are formed via the residues 
 at the core of the helices . in these dimeric states ( l  1 nm ) the system still explores both right - handed and left - handed 
 forms , through an equilibrium between the lzip and 
 the nter populations . the final stage of the simulation 
 shows the stabilization of the dimer in right - handed forms corresponding 
 to the global minimum of the fes ( figure 5c ) . the third simulation ( supporting information , figure s8 ) showed a mixed profile 
 with respect to the first two simulations , namely a fast association 
 of the helices followed by the exploration of both right - handed and 
 left - handed forms to reach the final dimeric state . interestingly 
 these results suggest that the equilibrium between the lzip and the nter populations may be a preliminary pathway 
 to access the most stable nter+ state . the 
 agreement between cg - metad and standard cg - md in finding the 
 lowest energy pathways in the monomeric , predimeric , and dimeric states 
 demonstrates that the metadynamics bias did not corrupt the thermodynamics 
 of the system , and that the cvs l and  properly 
 describe the motion of the helices during the dimerization / dissociation 
 processes . elucidating 
 the structural transitions that follow the initial dimerization step 
 is fundamental to our mechanistic understanding of egfr activation . in particular , the equilibrium dynamics of the tm domain may reveal 
 how the helix interactions change when the active and inactive dimeric 
 states interconvert . thus , we evaluated through standard ( i.e. , unbiased ) 
 cg - md simulations whether the four higher energy dimeric states identified 
 in the multidimensional landscape could convert into the nter+ conformation , which is thought to correspond to the activated state 
 of the receptor . first , we performed 
 a standard cg - md simulation ( of 5 s duration ) starting from 
 the cter+ conformation , which corresponds to a possible 
 inactive state of the receptor . this explored the various low energy conformations of the right - handed 
 ( negative  value ) region of the fes ( figure 6a ) . analysis and structural examination revealed 
 that there is equilibrium between the cter+ and the nter+ populations , in agreement with the free energy calculations 
 that indicated a low energy barrier between the two dimers ( 2 
 kj / mol from nter+ to cter+ ; see figure 3b ) . in particular , 
 we found that these two populations can interconvert through a  pivot  
 motion whereby the contacts formed by the n - terminal small - x3-small motif alternate with those formed by c - terminal motif ( figure 6b ) . a similar mechanism 
 has been invoked to explain the correlation between the tm conformational 
 transitions and the activation of the egfr intracellular kinase domains . ( a ) region of the fes explored during 
 5 s of standard cg - md 
 simulation starting from a cter+ conformation , and 
 ( b ) time - evolution of the d cv . ( c ) region 
 of the fes explored during 7 s of standard cg - md simulation 
 starting from a lzip conformation , and ( d ) time - evolution 
 of the || cv . in ( b , d ) conformations observed during 
 the simulations are indicated by the colored labels , and the insets 
 depict the observed pivot and rotational motions . comparable standard cg - md simulations starting from the ( left - handed ) lzip dimer showed a transition to the more stable right - handed nter+ ensemble ( figure 6c ) . this transition occurred through a two - step mechanism 
 involving helix rotation as characterized by evaluation of || 
 ( figure 6d ) . in particular , 
 after 1.5 s one of the two helices rotated leading 
 to the intermediate ( asymmetric dimer ) nter population . subsequently rotation of the second helix leads to the cter+ conformation that rapidly evolves to the active nter+ ensemble ( figure 6d ) . alongside some experimental studies , our results suggest that the  rotation  motion of 
 the tm helices of the egfr may correspond to another possible activation 
 pathway . we also performed cg - md simulations starting from cter and from nter ( supporting information figure s9 ) . we found a rapid transition 
 ( after 0.25 s ) 
 from the cter population to nter+ , which then existed in equilibrium with cter+ . this is in line with the free energy landscape previously described 
 ( figures 3b ) . starting 
 from the nter dimer , an equilibrium with the lzip population was reached and followed by transition to nter+ . these observations confirm our suggestion that nter is an intermediate state between lzip and nter+ ( figure 3c ) . some structural asymmetry , observed in nter conformations ( figure s9c ) , may play a role in the kinetics of conversion to the nter+ active state . taken together these 
 simulations support the picture suggested by , e.g. , nmr studies of a number of different structures of the tm 
 dimer , interconversion between which may provide pathways for activation 
 of the egfr . for a number of tm 
 helix dimers , our cg - metad simulations have allowed us 
 to identify nter+ as the most stable dimeric state . we therefore simulated 
 the dimeric states in membranes differing in their bilayer thickness 
 ( supporting information figure s10 ) to 
 investigate the effect of bilayer thickness on the conversion pathways 
 from the higher energy states ( cter+ , nter , cter , lzip ) to this stable 
 conformation . overall , these simulations confirm the results obtained 
 in the 16:0 pc bilayer used in the majority of the simulations . however , 
 some interesting differences were found . in particular , in a thinner 
 ( 12:0 pc ) bilayer the amplitude and the frequency of the pivot motions 
 were increased , whereas these appeared to be slowed in a thicker ( 20:0 
 pc ) bilayer ( figure 7a ) . similar behavior was observed regarding the rotational motions , 
 even leading in some case to the absence of interconversion in the 
 course of simulations ( figure 7b ) . these findings support the view that the bilayer thickness 
 can modulate tm helix dimer stability , and that the physical properties of plasma membrane 
 ( nano)domains may be able to modulate receptor activation . indeed , in future studies it 
 would be of interest to explore the dependence of free energy landscape 
 on more complex lipid bilayer compositions , including interactions 
 with , e.g. , cholesterol and phosphatidylinositol 4,5-bisphosphate 
 such as have been observed in , e.g. , simulations with gpcrs . influence of the bilayer thickness ( as shown schematically 
 in the 
 center panels ) on the kinetics of ( a ) pivot , and ( b ) rotational mechanisms . 
 in ( a ) d is monitored during cg - md simulations 
 starting from the cter+ conformation in either 12:0 
 pc or 20:0 pc bilayers . in ( b ) || is monitored during 
 cg - md simulations starting from the lzip conformation 
 in either 12:0 pc or 20:0 pc bilayers . the structure of the tm domain 
 of the egfr has been determined by nmr in dpc micelles ( pdb i d : 2m0b ) . it shows interhelix contacts at the c - terminal small - x3-small motif similar to those found in the cter+ population identified in our simulations ( figure 8a ) . furthermore , we note that the structure 
 of the erbb1/b2 tm helix dimer ( pdb i d : 2ks1 ) obtained by nmr in dmpc / dhpc bicelles 
 resembles more closely our nter+ ensemble . another nmr structure of the egfr tm domain 
 plus the adjacent intracellular juxtamembrane ( jm ) region ( pdb i d : 2m20 ) also revealed interhelix 
 contacts at the n - terminal small - x3-small motif within 
 the tm domain . this dimer structure is 
 very similar to the nter+ ensemble identified as 
 the population of lowest free - energy in our cg - metad calculations 
 ( figure 8b ) . this agreement 
 prompted us to undertake a series of 5 s standard cg - md simulations 
 using the tm+jm nmr structure ( pdb i d : 2m20 ) as a starting conformation . in these 
 simulations , we found that the experimental conformation is stable , 
 displaying contacts between the tm helices at the n - terminal motifs 
 with additional hydrophobic leucine - zipper interactions between the 
 jm segments . encouragingly , the coarse - grained 
 model reproduced the global structure and especially the interhelical 
 contacts determined by nmr ( figure 8c ) . ( a ) superposition of the nmr structure of the tm domain 
 ( pdb i d : 2m0b ) and a representative 
 structure of the cter+ ensemble . ( b , c ) superposition 
 of the nmr structure of the tm+jm dimer ( pdb i d : 2m20 ) with representative 
 structures of ( b ) the nter+ population ; and ( c ) cg - md 
 simulation # 1 in a 16:0 pc bilayer ( see supporting information figure s11a ) . nmr structures are in cyan , and the 
 backbone of cg models are shown as gray traces . red and green spheres 
 respectively represent the nter and cter motifs of interhelical contact in the tm domain of cg models . in 
 the insets we also 
 simulated the behavior of the nmr tm+jm dimer using cg - md in a phospholipid 
 bilayer containing anionic lipids ( popc : popg 85:15 ) in order to mimic 
 the environment presented by a mammalian plasma membrane . after 2.7 s a conformational transition was observed 
 from the nter+ to the cter+ population 
 ( supporting information figure s11a ) . this 
 change occurred through a pivot motion of the tm domain and was associated 
 with the separation of the jm segments ( figure 9a ) . interestingly a recent simulation study 
 suggests that a similar motion is associated with an inhibitory effect 
 of binding of the herceptin monoclonal antibody to the her2 receptor . in another set of cg - md simulations of the tm+jm 
 dimer , this time performed in a zwitterionic ( pc ) bilayer and starting 
 from the nter+ conformation , we observed a rotational 
 motion of the helices that led to the lzip conformation , 
 before returning to the initial nter+ conformation 
 ( supporting information figure s11b ) . similarly 
 to what was observed for the tm domain , the transition between the 
 right - handed and the left - handed conformations is mediated by an asymmetric nter conformation ( supporting information figure s11c ) . structural characterization also shows that the jm 
 leucine - zipper interactions are disrupted during this process and 
 reformed at the end of the simulation ( figure 9b ) . extrapolating back to the full - length 
 receptor , reformation of the antiparallel leucine - zipper interactions 
 in the jm region would be expected to latch the kinase domains in 
 the asymmetric active conformation . thus , 
 our cg - md simulations reveal reversible transitions between conformations , 
 which may be identified as corresponding to the likely active and 
 inactive states of the receptor , providing structural insights into 
 the dynamics of egfr activation . ( a ) structural summary of the dynamics 
 observed during a cg - md 
 simulation starting from the nmr structure of the tm+jm dimer ( pdb 
 i d : 2m20 ) in 
 a bilayer containing anionic lipids ( popc : popg 85:15 ) , indicating 
 pivot motion and average values of d for each 
 population . ( b ) structural summary of the dynamics observed during 
 a simulation in a 16:0 pc lipid bilayer , indicating rotation events 
 of individual monomer and average value of || for each 
 population . residues from motifs involved in direct interhelix contacts 
 are drawn as opaque spheres , while motif residues not engaged in direct 
 contacts are shown as transparent spheres . red , green and blue spheres 
 represent nter , cter and lzip motifs of the tm domain , and leucine zipper residues 
 in the jm domain are in pink . ( c ) suggested mechanisms of egfr activation 
 via either pivot or rotational motions . kinase domains of inactive and active states 
 are shown schematically with symmetric and asymmetric states , respectively , 
 in cyan and blue . overall , our simulations 
 suggest that when the tm helices are in contact through the c - terminal 
 small - x3-small or lzip motifs , the jm 
 segments are distant and can not readily form antiparallel leucine - zipper 
 interactions ( supporting information figure 
 s11 ) . this further suggests that the lzip , cter and cter+ conformations all correspond 
 to inactive states of the receptor . we note that all these inactive 
 dimer states exhibit an increased distance between the two k618 residues 
 in the n - loop that have been previously identified as critical for 
 regulation of egfr activation by glycolipids such as gm3 . we also note that such variation of distance 
 between charged residues in this specific upper region of the tm domain 
 was found sufficient to regulate activation of the growth hormone 
 receptor . thus , our findings provide 
 a description of the possible role of the tm domain in relaying conformational 
 changes between the extracellular and the intracellular regions of 
 the receptor ( supporting information figure 
 s11 ) . mechanism , both leading to the right - handed nter+ active dimeric state ( figure 9c ) . large - scale motions 
 of tm domains play a key role in the activation 
 mechanisms of membrane receptors such as rtks . these functional motion 
 of tm helices can be difficult to characterize by either experimental 
 or theoretical approaches . here , we have described an advanced computational 
 protocol that combines two well - established techniques , namely coarse - grained 
 molecular dynamics and well - tempered 
 metadynamics ( cg - metad ) , allowing 
 us to overcome the time scale limits of current simulations . using 
 this approach , we have simulated helix / helix interactions of the egfr 
 tm domain to reconstruct the underlying free energy landscape , providing 
 a detailed description of the conformational motions of this canonical 
 rtk tm dimer at near - atomic resolution . in particular , the enhancement 
 of the phase space exploration using cg - metad allows us to describe 
 states and transitions between states that would be otherwise only 
 partially described by more standard approaches . running large ensembles 
 of multiple standard cg - md simulations with different initial configurations would be an alternative approach but requires a priori structural knowledge on starting configurations 
 that is achievable only by a thorough sampling such as is possible 
 by cg - metad . the agreement between the simulated packing modes and 
 the available structural data supports our thermodynamic characterization . our calculations predict the existence of five dimeric states of the 
 egfr tm domain , which are separated from one another by less than 
 10 kj / mol ( i.e. , 4 kt ) , and for which the relatively fast 
 kinetics of interchange may be modulated by the thickness of phospholipid 
 bilayers ( figures 6 and 7 ) . these features underlie considerable 
 conformational flexibility of the tm domain , supporting its likely 
 role in the activation mechanism of the egfr . we have described two 
 mechanisms leading to the population of lowest free - energy , nter+ , corresponding to the activated state of the whole 
 receptor . the first involves a  pivoting  motion of 
 the helices that converts the cter+ ( inactive ) population 
 in the nter+ ( active ) population , while the second 
 is a  rotational  mechanism that leads from a left - handed lzip ( inactive ) population to the right - handed nter+ ( active ) population . simulations in the presence of the juxtamembrane 
 ( jm ) domain confirm this scenario and demonstrate the likely critical 
 role of the tm domain in the allosteric regulation of the egfr ( figure 9 ) . thus , we are able 
 to fully describe functionally relevant motions of the tm domain  switch  
 which links the ecto and intracellular domains within the intact receptor 
 dimer . it is of course worthwhile to consider possible limitations 
 of 
 the cg model ( i.e. , martini ) employed in these studies . the simplifications of the molecular systems 
 include a mapping of four non - hydrogen atoms into one single cg particle , 
 as well as the short - range nature of the interactions between particles . 
 these features lead to smoothed free energy landscape , which intrinsically 
 modifies the kinetics of the system . however , 
 the enthalpy - driven parametrization of the martini force field potential 
 aids prediction of energetically relevant states , including intermediate 
 states . it has been demonstrated that 
 the cg force field reproduces physicochemical properties of solvated 
 phospholipid bilayers and describes protein  lipid and protein  protein interactions . 
 to date in particular , the in - house variant 
 used in the present study has previously shown quantitative agreement 
 with experimental measures on a number of tm systems . furthermore , comparisons of pmfs 
 for helix / helix interactions of glycophorin a derived from martini 
 v2.1 vs v2.2 suggest that the overall well depth changes by less than 
 kt ( unpublished observations ) . although the standard formulations 
 of martini present limitations that may affect the dynamics of the 
 system , such as the lack of solvent polarizability and entropic components , the correlation between our results and available experimental data 
 suggests that this cg force field is able to predict thermodynamically 
 relevant structural properties ( figure 8) . using the cg - metad approach we have addressed 
 the challenging problem of characterizing 
 accurately the free energy 
 landscape of protein / protein interactions within membranes the association 
 and dissociation of the tm helices of egfr in a phospholipid bilayer 
 environment are described by a large number of interconversion events 
 between the dimeric and monomeric states ( see movie s1 ) . through the accelerated sampling of the slowest 
 motions of the system , cg - metad simulations allow us to overcome the 
 time scale limits of current simulations . estimating effective rates 
 of the observed events will be important for future developments of 
 the cg - metad methodology . a protocol has been recently proposed to 
 extract transition rates from atomistic metadynamics simulations . however , in the present case the kinetic calculations are rather 
 complex since the sampling acceleration observed with the cg model 
 along the cv space needs to be rigorously evaluated . by enabling 
 us to go beyond the time scale and sampling of current 
 simulations , our cg - metad studies alongside standard equilibrium cg 
 and at simulations reveal new structural information which complements 
 those obtained previously . our results demonstrate that the combination 
 of cg and metad provides an exhaustive exploration of the phase space 
 leading to an accurate description of the free energy landscape underlying 
 the association of the egfr tm helices . this allows identification 
 of new dimeric states and characterization of multiple , reversible 
 interconversion events between these states of the tm domain dimer 
 of the egfr . this contributes to our understanding of the mechanisms 
 of egfr signaling by helping to bridge the gap between the extracellular 
 and intracellular regions . the intrinsic functional motions of the 
 tm domain will reciprocally modulate and be modulated by additional 
 molecular interactions in the structurally intact receptor in vivo . future investigations 
 will need to analyze the possible effects of more complex membrane 
 environments and to probe the similarities 
 and/or differences among homologous rtks . this study endorses cg - metad as a powerful method that opens new 
 avenues in free energy calculations of biomolecules and in the investigation of large - scale conformational motions underlying 
 the activation mechanism of rtks and related membrane proteins .
the epidermal growth 
 factor receptor ( egfr ) is a dimeric membrane 
 protein that regulates key aspects of cellular function . 
 activation 
 of the egfr is linked to changes in the conformation of the transmembrane 
 ( tm ) domain , brought about by changes in interactions of the tm helices 
 of the membrane lipid bilayer . using an advanced computational approach 
 that combines coarse - grained molecular dynamics and well - tempered 
 metadynamics ( cg - metad ) , we characterize the large - scale motions 
 of the tm helices , simulating multiple association and dissociation 
 events between the helices in membrane , thus leading to a free energy 
 landscape of the dimerization process . 
 the lowest energy state of 
 the tm domain is a right - handed dimer structure in which the tm helices 
 interact through the n - terminal small - x3-small sequence 
 motif . 
 in addition to this state , which is thought to correspond to 
 the active form of the receptor , we have identified further low - energy 
 states that allow us to integrate with a high level of detail a range 
 of previous experimental observations . 
 these conformations may lead 
 to the active state via two possible activation pathways , which involve 
 pivoting and rotational motions of the helices , respectively . 
 molecular 
 dynamics also reveals correlation between the conformational changes 
 of the tm domains and of the intracellular juxtamembrane domains , 
 paving the way for a comprehensive understanding of egfr signaling 
 at the cell membrane .