protein is intrinsically 
 dynamical , and conformational changes 
 are essential to its functions .
large - scale conformational rearrangements 
 are involved in signaling transduction , enzyme catalysis , and protein 
 folding .
molecular dynamics ( md ) simulations 
 are a powerful tool to model such motions and deepen our understanding 
 of the relationship between protein structure and function
. studies of protein conformational changes and 
 the intermediates that are formed along the transition pathway provide 
 valuable insights into the process by which the protein is activated 
 or inactivated . however , characterization of the transition pathway 
 and intermediates as well as establishing their role in conformational 
 rearrangements remains a major challenge . to tackle this problem , 

a plethora of methods have been developed 
 to enhance sampling of the conformational change , which is a rare 
 event .
it is , however , still difficult 
 to obtain adequate sampling of conformational changes which usually 
 occur on microsecond to millisecond time scales or beyond , especially 
 for a large biomolecular system . to overcome this limitation ,
one approach aims to speed up 
 conformational transitions by use of elevated temperatures or a biasing 
 potential such as parallel tempering or replica exchange md ( remd ) , umbrella sampling , metadynamics , accelerated md , temperature - accelerated md , adaptive biasing force , and so on .
remd and its variants are popular in sampling rare events 
 while its efficiency has been discussed in detail . for a bias potential , a finite number of collective 
 variables ( order parameters ) needs to be defined for description of 
 transitions .
the other 
 class of rare events method , path - based sampling , is of greatest relevance 
 to the present work .
these path - building techniques have been developed 
 to construct a path between two end points , including milestoning , the string method , path metadynamics , transition path 
 sampling , weighted ensemble , fast marching method , nudged elastic 
 band , targeted md ( tmd ) , steered md , and so on .
these approaches have 
 been applied for exploration of large - scale conformational changes 
 of biomolecules .
however , fundamental questions still remain to be 
 answered . how can these pathways be compared ?
we present a strategy combining 
 multiply targeted md ( mtmd ) , 
 unbiased md simulations , and bayesian clustering to address these 
 questions .
this approach is helpful not only for exploration of transition 
 pathways and metastable states along them but also for comparison 
 of distinct pathways . in our method , the phase space is first partitioned 
 by a series of seeding structures , which are generated by mtmd .
then , 
 a long unbiased md simulation ( 200 ns ) is initiated on each seeding 
 structure to achieve local equilibration and explore the subspace . 

the ensemble of md trajectories is analyzed by a clustering method 
 based on the bayesian model to determine metastable states and transitions 
 between them .
the strategy of running multiple md trajectories for 
 sampling rare events has been investigated previously .
the improvements in our approach are that , ( i ) with mtmd , we can 
 define transition pathways from an initial state to a final state 
 via different intermediate states ; ( ii ) the seeding structures are 
 not necessary to be in equilibrium ; ( iii ) long unbiased md simulations 
 are performed to explore likely transitions between neighboring seeding 
 structures ; ( iv ) specification of collective variables is not needed ; 
 and ( v ) bayesian clustering is carried out to reconstruct conformational 
 transitions from the ensemble of md trajectories . with these features , 
 it is likely to design different transition pathways between two end 
 points and compare them .
the initial structure , final structure , and 
 intermediate structures can be obtained from previous crystallographic 
 studies or md simulations . to test this approach ,
we apply it 
 to the activation / inactivation 
 process of the kinase domain of the epidermal growth factor receptor 
 ( egfr ) .
egfr , a receptor tyrosine kinase , has emerged as a target 
 of choice for treatment of various cancers .
crystallographic studies have revealed different conformational 
 states of egfr kinase domain ( supporting information 
 ( si ) figure s1 ) : the dfg - in active state ( a ) , the src - like 
 inactive state ( t ) , and the dfg - out inactive state ( i ) . for simplicity , 
 we will use a , t , and i in the following to represent these three 
 forms , respectively .
the major structural changes among them are movement 
 of the activation loop ( a - loop ) , dfg - flip ( conversion between dfg - in 
 and dfg - out ) , and displacement of c helix . in form
a , it adopts 
 dfg - in conformation ; the a - loop is in the extended state ; glu762 located 
 in c helix forms a salt bridge with lys745 , which is highly 
 conserved and important for kinase activity . form
i adopts dfg - out conformation , and the a - loop is also extended 
 while lys745 and glu762 ( the ion pair ) are separated because of outward 
 rotation of the c helix .
form t adopts the dfg - in conformation , 
 and the ion pair are apart from each other while the a - loop is in 
 the compact mode .
the mechanism that controls the transition from 
 inactive forms of kinases to active ones as well as intermediate states 
 has attracted much interest . based on analysis of abl crystal structures and molecular simulations , 
 it was proposed that the src - like inactive form ( t ) functioned as 
 an intermediate in the activation pathway to facilitate dfg - flip .
however , it is still not clear if dfg - flip can happen without passing 
 form t. if it is likely , is any intermediate state involved which 
 has not been determined experimentally ?
in addition , the roles of 
 a - loop and c helix in this process need to be elucidated .
deep 
 understanding of the activation / inactivation pathways and intermediate 
 states of kinases will be helpful for development of high specific 
 anticancer drugs .
in this work , we design and compare two distinct 
 transition pathways 
 of egfr kinase domain to investigate under what condition dfg - flip 
 can happen .
one is the conformational change between a and i via t , 
 which will be called ati pathway in the following .
the other one is 
 direct transition between a and i without any intermediate , which 
 will be called ai pathway .
then , in comparison with crystal 
 structures , several intermediate states are recognized based on clustering 
 results .
finally , rough two - dimensional ( 2d ) free energy surfaces 
 are plotted to illustrate transitions between different conformations . 

comparison of two predefined pathways yields new insights into the 
 conformational transition of egfr kinase domain .
moreover , our results 
 verify that this approach is valid and efficient in sampling of large - scale 
 conformational changes and determination of the energy - favorable pathway 
 from available ones .
two mutations in the structure of 4i20 ( l858r and v948r ) 
 were changed back to the wild type within maestro v9.3 .
some residues at the terminus were 
 deleted so that the three structures had the same 284 residues .
the targeted molecular 
 dynamics ( tmd ) and multiply targeted molecular dynamics ( mtmd ) implemented 
 in the amber 11 program suite were employed 
 to construct the transition pathways between egfr active and inactive 
 states .
it has the following form : tgtmdfrc 
 is the force constant for targeted md .
nattgtrms is the 
 number of atoms specified for root mean square deviation ( rmsd ) calculation . 

in contrast , the force constant in mtmd is changing 
 with step number by a single factor r : increments is the number 
 of times the target 
 value changes .
all three structures were subjected 
 to energy 
 minimization first using the steepest descent algorithm in 1000 steps 
 and then the conjugate gradient algorithm in 10000 steps . for tmd 
 simulations ,
the starting structure is 2itw ( a ) and the target structure is 3w32 ( t ) .
three tmd simulations 
 ran for 0.5 ns with different force constants ( 1 , 0.1 , and 0.01 kcal/(mol ) ) .
for mtmd , each simulation ran for 1 ns . for the i - t - a 
 transition , form i ( 4i20 ) was used as the starting structure ; form t ( 3w32 ) was first set as 
 the target structure .
the force constant was applied to all heavy 
 atoms and varied from 0.0001 to 2.0 kcal/(mol ) in 0.5 ns as the simulation steps increased to drive the conformational 
 changes to form t. then , form a ( 2itw ) was set as the target structure .
the 
 force constant was reset to 0.0001 kcal/(mol ) and increased to 2.0 kcal/(mol ) in continuous 
 0.5 ns to drive the further transition to form a ( 2itw ) . the same setting 

was also applied in the simulations of a - t - i and i - a - i .
the tmd and 
 mtmd simulations were performed in the implicit solvent model of generalized 
 born using the sander module .
snapshot 
 structures were extracted from produced trajectories at a time interval 
 of 1 ps . to generate an ensemble of representative structures from each mtmd 
 trajectory ,
snapshot structures from each mtmd trajectory were superposed onto 
 the reference structure using backbone atoms to remove overall rotation 
 and translation .
then rmsd values were calculated with respect to 
 the reference structure as the distance metrics .
three mtmd trajectories 
 were divided into six transition trajectories each lasting for 0.5 
 ns . for each transition , 500 snapshot structures were extracted and 
 classified with the k - means method implemented in 
 the statistical software r v3.0 . for ati pathway ,
four transition 
 trajectories were included and each of them was clustered into 11 
 classes . for ai pathway , two transition trajectories were included 
 and each of them was clustered into 9 classes . from each class ,
the 
 structure closest to the centroid was chosen as the representative 
 structure for following md simulations .
the equations 
 of motion were solved with the leapfrog integration algorithm with 
 a time step of 2 fs .
the lengths of all bonds involving hydrogen atoms 
 were kept constrained with the shake algorithm .
the particle mesh ewald ( pme ) method was applied for treating 
 long - range electrostatic interactions .
a random 
 seed was generated based on the current date and time for every run 
 to assign initial velocities .
for each 
 seeding structure , it was explicitly solvated in a rectangular box 
 of tip3p water molecules with a minimal 
 distance of 10  from the protein to the box edges .
after the whole 
 system was set up , a series of energy minimizations and equilibrations 
 were performed .
first , the water molecules , hydrogen atoms , and ions 
 were subjected to 3000 steps of steepest descent minimization followed 
 by 12000 steps of conjugate gradient minimization while other heavy 
 atoms were constrained with the harmonic force of 2 kcalmol .
next , the whole 
 system was energy minimized with 10000 steps of l - bfgs algorithm without 
 any harmonic restraint .
then , coupled to a langevin thermostat , the 
 system was heated from 10 k up to 300 k by increments of 100 k in 
 20 ps and continued to run for 40 ps at 300 k at constant volume . 

finally , the system was equilibrated for 200 ps in npt ensemble with the langevin thermostat and isotropic position scaling , 
 at 300 k and 1 bar .
the production run for each seeding structure 
 was carried out for 200 ns in nvt ensemble with the 
 langevin thermostat at 300 k using the parallel cuda version of pmemd 
 on 2 gpus .
all simulations were performed on our linux cluster with 2 gpus and 
 12 cpus on each node .
bayesian clustering is an unsupervised 
 method that seeks a maximum 
 posterior probability classification by an expectation - maximization 
 ( em ) algorithm .
the advantage is that there is no need to specify 
 the cutoff radius or the number of clusters in advance .
the optimal 
 solution is determined by calculating the maximum posterior parameter 
 values and the most probable classification .
the structure of the egfr kinase domain 
 was described by 282 pairs of backbone dihedral angles .
the structures 
 for bayesian clustering were sampled from unbiased md trajectories . 
 for ati pathway ,
snapshot structures were extracted at a time interval 
 of 0.4 ns from each trajectory and a total of 23500 structures was 
 obtained . for ai pathway ,
snapshot structures were taken at an interval 
 of 0.2 ns from each trajectory and a total of 20000 structures was 
 obtained .
the clustering process started with a serial of seed numbers 
 which gave an initial guess of the number of clusters .
then a random 
 classification was generated and refined until a local maximum was 
 found .
this process was repeated 1000 times , and the top two classifications 
 with the highest probability were output .
the 2d free energy surfaces 
 based on these sampled structures were plotted ( supporting information ) .
the calculation 
 of free energy is given bywhere w0 is the 
 depth of the free energy surfaces ( fes ) , p is the 
 probability distribution , kb is the boltzmann 
 constant , and t is the temperature .
all sampled structures 
 from unbiased md simulations are employed for 2d fes calculations 
 and 1d convergence analysis . for 2d fess , the 2d distribution function p(r )
is computed by the histogram analysis 
 method hist2d implemented in the package of gplots with a bin area 
 of 0.3   0.3  .
one - dimensional free energy surfaces 
 are presented in the supporting information , which are calculated over the three metrics ( a - loop rmsd , dfg rmsd , 
 and ion - pair distance ) at an interval of 40 ns of simulation time . 

the 1d distribution function is calculated by the histogram analysis 
 method implemented in r v3.0 with a bin width of 0.1  .
compared to other approaches , tmd 
 is 
 a valid tool for construction of a transition pathway , which does 
 not require the specification of collective variables . in the standard 
 tmd method
one is that a weak 
 constraint may hold back the trajectory from reaching the target . 

the other one is the so - called  large - scale - first  , which means that large - scale motions occur before 
 small - scale ones .
thus , the order of events in conformational transition 
 is associated with the direction in which the simulation is performed . to avoid these questions ,
we employ mtmd , a tmd variant , to generate 
 the transition pathway for sampling of seeding structures .
compared 
 with tmd , mtmd allows not only multiple reference structures being 
 used but also a constraint evolving in time .
the mtmd simulations 
 were performed with the force constant varying from very weak ( 0.0001 
 kcal/(mol ) ) to weak ( 2.0 kcal/(mol ) ) to guarantee that the final structure would be within heavy - atom 
 rmsd 1.0  of the target structure and large - scale motions would 
 not happen at the beginning .
we also ran mtmd simulations in both 
 activation and inactivation directions to neutralize the possible 
 bias in one - direction simulations .
two of them are for ati pathway : 
 one is from state i to a via t and the other one is from a to i via 
 t. the third trajectory is for ai pathway : from state i to a and then 
 back to i. in the mtmd simulations the starting structure is gradually 
 changed into the target structure . and dfg - flip occurs after 0.4 ns 
 in each transition ( si figure s2 ) .
in contrast , 
 three tmd trajectories with different constraints are also shown in 
 figure 1 . in tmd simulations
this simple comparison indicates 
 that mtmd produces a trajectory of conformational changes more smooth 
 and consecutive than tmd does .
detailed comparison of mtmd and tmd 
 is beyond the scope of this work and will be left for future studies .
three mtmd simulations are the conformational transition from state 
 i to a via t ( i - t - a ) , the transition from state a to i via t ( a - t - i ) , 
 and the transition between states i and a ( i - a - i ) , respectively .
conformational 
 transition from state a to t is also simulated by tmd with different 
 force constants ( 1.0 , 0.1 , and 0.01 kcal/(mol ) ) .
three crystal structures , 2itw , 3w32 , and 4i20 are the dfg - in active form ( a ) , src - like 
 inactive form ( t ) , and dfg - out inactive form ( i ) , respectively . the 
 rmsd values are calculated with respect to the target structure ( 2itw , 3w32 , or 4i20 ) .
the sampled structures from mtmd simulations were 
 clustered using k - means method based on backbone 
 rmsds of both egfr kinase 
 domain and a - loop , considering that a - loop displays the highest flexibility 
 in crystal structures .
the structure closest to the centroid in one 
 cluster was selected as the representative structure . for ati pathway , 

44 snapshot structures were chosen and 47 seeding structures ( plus 
 3 crystal structures ) were prepared for md simulations . for ai pathway , 
 18 snapshot structures were selected and 20 seeding structures ( plus 
 2 crystal structures ) were prepared for md simulations . the structural fluctuation 
 from the ensemble of md trajectories in two distinct pathways
is characterized 
 by four metrics : backbone rmsd of egfr kinase domain ( backbone rmsd ) , 
 backbone rmsd of the a - loop ( residues 854878 , a - loop rmsd ) , 
 heavy - atom rmsd of the dfg motif ( residues 855857 , dfg rmsd ) , 
 and distance between the center of mass of lys745 and glu762 ( the 
 ion - pair distance ) .
backbone rmsd describes the overall conformational 
 changes of the egfr kinase domain ; a - loop rmsd evaluates the movement 
 of the a - loop between the extended state and the compact mode ; dfg 
 rmsd indicates the conformation of the dfg motif ( dfg - in or dfg - out ) ; 
 the ion - pair distance ( lys745-glu762 ) reflects displacement of c 
 helix . as a reference , the ion - pair distance and rmsd values of three 
 crystal structures with respect to the crystal structure of 2itw ( form a ) are listed 
 in table 1 .
distribution of the measurements 
 sampled from md trajectories is displayed in figure 2 , in which the four metrics have different features of distribution 
 in two pathways , suggesting that sampling is affected by predefined 
 pathways .
as expected , the a - loop rmsd varies more widely in ati pathway 
 than in ai pathway , though the overall structure ( backbone rmsd ) does 
 not change much in both pathways .
the dfg rmsd and ion - pair distance 
 fluctuate in the same range with distinct characterization of distribution 
 in two pathways , suggesting the influence of designed pathways on 
 their movements .
distribution of rmsds and ion - pair distance : ( a ) data 
 from 47 md 
 trajectories in ati pathway ; ( b ) data from 20 md trajectories in ai 
 pathway .
the dfg - in active conformation ( 2itw ) is used as the reference structure . 

then rmsd values of backbone , a - loop , and dfg 
 are calculated .
the density is computed by the histogram analysis 
 implemented in r v3.0 with a bin width of 0.1  .
the crystal structure of 2itw ( a ) is used as the 
 reference structure for rmsd calculations .
clustering based on the bayesian model was previously employed for 
 analysis of md trajectories .
an analysis 
 of clustering simulation data indicated that the clustering results 
 to some extent depended on the number of available sampled structures .
considering their results and computation time , 
 we extracted 23500 structures for ati pathway ( 500 structures from 
 each md trajectory ) and 20000 structures for ai pathway ( 1000 structures 
 from each md trajectory ) for bayesian clustering .
the egfr structure was described 
 by 282 pairs of backbone dihedral angles . for ati pathway ,
23500 structures 
 were grouped into 94 classes ( class 093 ) in the optimal classification . 
 for ai pathway , 85 classes ( class 084 )
the four metrics aforementioned are employed 
 for comparison of sampled structures to the crystal structures ( for 
 the full data of the top 20 populated classes , see si tables s1 and s2 ) . by them
, we not only quantify the overall 
 agreement of structures in classes with crystal structures but also 
 identify intermediate structures classified by bayesian clustering , 
 which has not been reported yet . for ati pathway ,
structures in classes 
 0 , 2 , and 11 have very good agreement with all structural properties 
 of forms t , i , and a , respectively , with average rmsd values below 
 3.0  and similar ion - pair distances ( for selected data , see 
 table 2 ; for full data , see si table s1 ) . except the known structures , conformations in 
 some other classes
for example , structures in class 7 exhibit backbone 
 and a - loop rmsds as well as an ion - pair distance close to form t , 
 but adopt the dfg - out conformation .
similar properties are also found 
 in classes 9 and 16 . for ai pathway ( for selected data ,
see table 2 ; for full data , see si table s2 ) , structures in class 0 have good agreement with the structural 
 properties of form i. structures in class 3 exhibit properties similar 
 to form a , typically with a good agreement to the dfg - in conformation , 
 backbone rmsd , and ion - pair distance as well as a better correspondence 
 to the a - loop rmsd .
structures in class 1 have similar 
 backbone rmsd and a - loop rmsd to form i while dfg rmsd and the distance 
 of the ion pair are close to form t. structures in class 9 adopt the 
 dfg - out conformation , and rmsds of the backbone and a - loop are closer 
 to form i than to form a while the ion pair is likely to be in contact , 
 which makes them look like the imatinib - induced abl structure .
although this structure has been solved for 
 some kinases , it has not been observed in egfr structures experimentally 
 determined .
further 
 examination of fluctuation of the structural properties 
 shows that egfr undergoes conformational transitions in some classes , 
 though it is pretty stable in most of the classes ( for the top 20 
 populated classes , see si figure s4 ) . for 
 example , in class 1 of ati pathway the conformational change between 
 dfg - out and dfg - in is clearly observable ( figure 3a ) .
it is coupled with movement of the a - loop and alteration 
 of the ion - pair distance .
similarly , structures in class 13 of ai 
 pathway also exhibit conformational change of the dfg motif while 
 the a - loop does not move far away from the extended state ( figure 3b ) .
these kinds of changes in adopted structural 
 properties occur more frequently in less - populated classes than in 
 populated ones , suggesting transitions between intermediate states .
evolution 
 of structural properties in the class from bayesian clustering : 
 ( a ) class 1 in ati pathway ; ( b ) class 13 in ai pathway .
four metrics 
 are displayed : backbone rmsd , a - loop rmsd , dfg rmsd , and distance 
 of the ion pair ( lys - glu ) .
the rmsd values are calculated with respect 
 to the crystal structure of 2itw ( form a ) . to visualize the transition pathways and intermediate states , we 
 tried to project the dihedral angles space into a two - dimensional 
 fes .
the principle component analysis was carried out , but the first 
 two eigenvalues 
 could only account for 13% of the variability of the data , which made 
 them not a good choice .
three metrics ( dfg rmsd , a - loop rmsd , and distance of the ion pair ) 
 were employed for plotting of rough 2d fess to investigate the relationship 
 between dfg - flip and movement of the a - loop as well as displacement 
 of the c helix ( figure 4 ) .
fess for ati pathway are shown in a and b ; fess for ai pathway 
 are shown in c and d. three metrics are employed for fes plotting : 
 dfg rmsd , a - loop rmsd , and distance of the ion pair .
the rmsd values 
 are calculated with the crystal structure of 2itw ( form a ) as the 
 reference .
energy basins corresponding to forms a , t , and i are denoted . 
 intermediate states are labeled with prefix is .
energy wells on fess are compared to the classes from bayesian 
 clustering . for ati pathway , three populated ensembles labeled as 
 t ,
i , and a ( figure 4a , b ) correspond to classes 
 0 , 2 , and 11 ( table 2 and si figure s3 ) , displaying similar features with the crystal 
 structures .
is1 , is2 , and is3 ( figure 4a ) denote 
 three intermediate states , corresponding to classes 7 , 31 , and 14 
 ( table 2 ) , respectively .
states a , t , and is3 
 adopt the dfg - in conformation while i , is1 , and is2 adopt the dfg - out 
 conformation . for ai pathway , energy basins labeled as i and a ( figure 4c , d ) correspond to classes 0 and 3 , having properties 
 similar to those of the crystal structures .
is4 , is5 , and is6 denote 
 three different forms ( figure 4d ) , corresponding 
 to classes 1 , 9 , and 7 , respectively ( table 2 ) .
states i , is5 , and is6 adopt the dfg - out conformation while a 
 and is4 adopt the dfg - in conformation . from fess ,
the transition 
 between egfr active and inactive forms 
 is apparently visible . in ati pathway ,
dfg - flip occurs with little 
 energy cost ( figure 4a , b ) .
in contrast , dfg - flip 
 happens less frequently in ai pathway , which is clearer in figure 4d than in figure 4c .
our 
 results suggest that ati pathway is more energetically favorable than 
 ai pathway , though dfg - flip can happen in both pathways .
previous 
 studies suggest that displacement of the c helix or a - loop movement is the major energy 
 barrier mediating the transition . here
, we further examine the role 
 of a - loop and c helix in dfg - flip . in both b and d of figure 4 , dfg - flip occurs when the ion pair is broken apart . 

when the ion pair is in contact , no transition of the dfg motif is 
 observed , possibly because of sterical hindrance .
this suggests that 
 dfg - flip is coupled with disruption of the ion pair . on the other 
 hand , in figure 4c
the a - loop fluctuates around 
 the extended state , while , in figure 4a , the 
 a - loop floats away from the extended state and reaches the compact 
 mode .
dfg - flip happens more frequently in figure 4a than in figure 4c , suggesting that 
 the transition is strongly influenced by movement of the a - loop .
thus , 
 our results yield a view that breaking the ion pair is critically 
 important for dfg - flip and movement of a - loop significantly facilitates 
 the conformational change .
interestingly , except 
 the predefined pathway ( a - t - i ) , another transition 
 pathway appears between form a and i ( a - is3-is2-i in figure 4a ) . in this pathway ,
we deduce that disruption of the ion 
 pair plus a small movement of a - loop provides enough space for dfg - flip . 

if this is the case , considering that the a - loop displays high flexibility 
 in kinase crystal structures and distinct a - loop conformations of 
 egfr have been recognized , it is very 
 likely that there exist multiple transition pathways between the active 
 and inactive forms as well as various intermediate states . in 
 addition
, we notice that state is5 appears in ai pathway ( figure 4d ) but is not sampled in ati pathway ( figure 4b ) .
examination of clustering results also shows 
 that none of the classes in ati pathway has structural properties 
 similar to those of is5 .
this may suggest that some conformation is 
 not preferred for unbound egfr . in ati pathway ,
the inactive forms 
 ( t and i ) and the conversion between them are energetically favorable 
 ( figure 4b ) , and sampling is trapped in the 
 energy wells .
this agrees with a previous report that the isolated 
 egfr kinase domain in monomeric form has low catalytic activity . in ai pathway , however , transition of the dfg 
 motif is impeded because of constraint of the a - loop ( figure 4d ) , which may make state is5 energetically favorable . 

our results are consistent with the observation that the imatinib - induced 
 abl structure ( corresponding to is5 in figure 4d and class 9 in si figure s3b ) is rarely 
 found in unbound kinases but appears in kinase  ligand complexes . 

possibly because this form is not energetically preferred and becomes 
 stable under certain conditions such as ligand binding or constraint 
 of the a - loop .
previously , conformational transitions of src 
 and abl have been 
 investigated with biased md simulations and the free energy may be 
 overestimated .
compared to their data , 
 the free energy calculated in this work may be underestimated . interestingly , 

our results are close to the free energy computed for conformational 
 changes of wide type and mutant egfrs recently reported .
the difference in estimated free energy may 
 result from distinct kinases , simulation methods or choice of order 
 parameters . to improve performance of our approach , more seeding structures 
 may be needed .
studies 
 have demonstrated that local protein unfolding ( cracking ) is common 
 in conformational changes , and c helix is disordered in egfr 
 conformations .
based on previous reports , displacement of c helix may be 
 caused by outward rotation or local unfolding .
outward rotation of 
 c helix among kinase crystal structures has been discussed 
 in detail . here
we further examine variability 
 of secondary structure of c helix by calculating the time when 
 residues of c helix are in the helix status ( si figure s5 ) . in all md trajectories three residues ( residue 
 nos .
disorder 
 transitions . as indicated in the previous section , disruption of the 
 ion pair induced by displacement of the c helix
our results suggest that partial unfolding of he c 
 helix may accompany the transition between egfr active and inactive 
 forms and facilitate dfg - flip .
we present a strategy combining md simulations and bayesian clustering 
 to investigate conformational transitions between egfr active and 
 inactive forms . compared to previous reports ,
our results verify that 
 this approach is valid and efficient in sampling rare events .
moreover , 
 with this approach , distinct transition pathways can be compared and 
 the energy - favorable one can be distinguished .
it is usually difficult 
 to observe the conformational changes using constant - temperature simulations 
 even on microsecond time scales . running
multiple md simulations initiated 
 on different conformational states may effectively cross the energy 
 barrier and overcome the conformational sampling problem . with the 
 help of a high parallel computing platform , this strategy is easily 
 implemented .
bayesian clustering is then employed to recognize intermediate 
 states and transitions along the pathway from the ensemble of md trajectories . 

although it is still a challenge to project the transition pathway 
 from a high - dimensional space to a two - dimensional space , selection 
 of appropriate order parameters is helpful for investigating the effect 
 of local fluctuations on global conformational changes .
since 
 the catalytic domain is highly conserved in all kinases and 
 the three forms of egfr ( a , t , and i ) have also been found in other 
 distantly related kinases , the conformational changes between them 
 may be a general mechanism of kinases activation and inactivation . on the other hand , it is believed that each 
 inactive conformation of kinase has some special features which distinguish 
 it from other states . and kinase inactive forms have been the targets 
 for development of high specific anticancer drugs in the pharmaceutical 
 industry .
therefore , determination of unreported inactive conformations 
 and transitions between them provides new targets for cancer treatment .
 the 
 epidermal growth factor receptor ( egfr ) is aberrantly activated 
 in various cancer cells and an important target for cancer treatment . 
 
 deep understanding of egfr conformational changes between the active 
 and inactive states is of pharmaceutical interest . 
 here we present 
 a strategy combining multiply targeted molecular dynamics simulations , 
 unbiased molecular dynamics simulations , and bayesian clustering to 
 investigate transition pathways during the activation / inactivation 
 process of egfr kinase domain . 
 two distinct pathways between the active 
 and inactive forms are designed , explored , and compared . 
 based on 
 bayesian clustering and rough two - dimensional free energy surfaces , 
 the energy - favorable pathway is recognized , though dfg - flip happens 
 in both pathways . 
 in addition , another pathway with different intermediate 
 states appears in our simulations . 
 comparison of distinct pathways 
 also indicates that disruption of the lys745-glu762 interaction is 
 critically important in dfg - flip while movement of the a - loop significantly 
 facilitates the conformational change . 
 our simulations yield new insights 
 into egfr conformational transitions . moreover 
 , our results verify 
 that this approach is valid and efficient in sampling of protein conformational 
 changes and comparison of distinct pathways . 
Introduction
Methods
Results
Conclusion