Efficient and Accurate Modeling of Conformational Transitions in

Subscriber access provided by UNIV OF DURHAM

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Efficient and Accurate Modelling of Conformational Transitions in Proteins: The Case of c-Src Kinase Edoardo Milanetti, Alexandra G. Trandafir, Josephine Alba, Domenico Raimondo, and Marco D'Abramo J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b07155 • Publication Date (Web): 05 Sep 2018 Downloaded from http://pubs.acs.org on September 5, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Efficient and Accurate Modelling of Conformational Transitions in Proteins: The Case of c-Src Kinase Edoardo Milanetti1,2, Alexandra G. Trandafir3, Josephine Alba3, Domenico Raimondo4, Marco D’Abramo3* 1

Department of Physics, Sapienza University of Rome, P. le A. Moro, 5, 00185 Rome, Italy

2

Center for Life Nano Science @ Sapienza, Istituto Italiano di Tecnologia, Viale Regina Elena

291, 00161 Rome, Italy 3

Department of Chemistry, Sapienza University of Rome, P.le A. Moro, 5, 00185 Rome, Italy  

4

Department of Molecular Medicine, Sapienza University of Rome, Viale Regina Elena 324,

00161 Rome, Italy

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 30

ABSTRACT The theoretical computational modelling of large conformational transitions occurring in biomolecules still represents a challenge. Here we present an accurate “in-silico” description of the activation and deactivation mechanism of human c-Src kinase, a fundamental process regulating several crucial cell functions. Our results clearly show that, applying an efficient and automated algorithm able to drive the molecular dynamics sampling along the pathway between the two c-Src conformational states – the active and the inactive state – it is possible to accurately describe, at reduced computational costs, the molecular mechanism underlying these large conformational rearrangements. This procedure, combining molecular dynamics simulations with the sampling along well defined principal motions connecting the two conformational state, allows to provide a description well beyond the present computational limits and it is easily applicable to different systems where the structures of both the initial and final state are known. 1. INTRODUCTION The Src family kinases (SFKs) are the main group of non-receptor tyrosine kinases. They are enzymes involved in different regulation processes with important roles in the signal transduction following engagement of many different classes of cellular receptors. They control cell growth, differentiation, migration, cell survival and their activity is regulated by the phosphorylation status of two tyrosine residues in key positions and by ligand interaction through kinase domains. These interactions stabilize the kinase in active or inactive form. A dysfunction in SFK kinase control mechanisms can cause a change in their activation which is


2

Page 3 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


directly related to the onset of many pathologies such neoplastic progression1,2, preeclampsia3, and immune disorders4. The SFKs are proteins of 52-62 kDa, sharing a conserved modular domain organization as we reported in Figure 1. The human Src gene encodes 536 amino acids (UniProtKB/Swiss-Prot accession code: P12931). From the N-terminus to C-terminus direction, c-Src protein is composed by an SH3 domain, an SH2 domain, an SH2-kinase linker portion, a protein-tyrosine kinase domain (SH1), and a C-terminal regulatory portion (Fig. 1). The c-Src kinase domain – which was investigated here – is a well conserved protein when we consider SFK family proteins or other tyrosine kinases; moreover, it is both indispensable and sufficient for enzymatic activity. The kinase domain (residues 270-523 according UniProtKB/Swiss-Prot annotation) is composed by the N-terminal region (N-lobe) and the C-terminal region (C-lobe). Active site includes 1) the glycine-rich loop region and C-helix from N-lobe 2) the catalytic and activation loops from the C-lobe. Enzymatic activity is regulated by the conformational variation and three-dimensional rearrangements of these structural elements, which in turn are modulated by phosphorylation and interaction with other regulatory proteins (Fig. 1). A crucial residue involved in c-Src auto phosphorylation is Tyr1655. Auto phosphorylation of catalytically inactive kinases can occur, demonstrating that phosphorylation must proceed in trans6. SH2 and SH3 domains (residues 151248 and 84-145, respectively) are crucial elements substrate recruitment and in the catalytic activity modulation7. In physiological condition c-Src enzyme is almost inactive in cells, while dramatic changes in its activity occurs when the fine balance between phosphorylation and dephosphorylation is altered. However, if the fine balance between phosphorylation and dephosphorylation is disrupted, changes can occur in c-Src activity with drastic results8,9. c-Src contains a catalytic loop, (residues 385 to 403) formed by Asp and Asn residues that are highly


3


Page 4 of 30

conserved and very important for the Mg-ATP complex coordination. Residues 404 to 432 form the activation loop. This region plays a key role in kinase activity because it is involved in the phosphorylation of Tyr416, a protein modification required for full enzymatic activation. In particular, c-Src lives in “open” active and “closed” inactive conformations that are structurally well understood, a behavior similar to others Src family members10. In order to bind Mg-ATP and the protein substrate the enzyme has to assume the open form. Starting from this open conformation the protein starts converting to the closed one as the catalysis occurs. As the catalysis ends the enzyme is reconverted to the open form, phosphorylated protein and Mg-ADP are released and the next catalytic cycle can occur. The open form of the active enzyme binds Mg-ATP and the protein substrate; this is accompanied by the conversion to the closed form as catalysis occurs. After catalysis, phosphorylated protein and then Mg-ADP are released as the enzyme is reconverted to the open form prior to the next catalytic cycle. The two main phosphorylation sites on c-Src are represented by Tyr416 and Tyr527. Phosphorylation of Tyr527 is one of the mechanisms for the regulation of c-Src activity. When Tyr527 is phosphorylated the molecule adopts a folded, inaccessible bundle and the inactive configuration is mediated by an intra-molecular interaction between the phosphorylated. Tyr527 and the SH2 domain of c-Src11. This “inactive conformation” is stabilized by a network of interactions between several residues. In particular, Asp404 forms electrostatic interaction with Lys295 and in this way it disrupts Lys295-Glu310 interaction, essential for maintenance of the active state. As consequence of Asp404 – Lys295 interaction, Phe424 moves toward the α-C-helix of the c-Src N-lobe, forcing the helix to move away from the catalytic cleft12. Moreover, Tyr416 is oriented towards the catalytic center of the kinase domain when residues 413-418 of the activation loop form a short α-helix (A-loop helix). In this way Tyr416 not only is able to form hydrogen bonds with Arg385


4

Page 5 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and Asp386 but it is also buried into a hydrophobic cavity in order to be protected from autophosphorylation13. In addition, the active conformation of the α-C-helix is stabilized by the formation of a salt bridge between Glu310 of the α-C-helix and Arg40914,15. When Tyr527 becomes dephosphorylated, this bond is disrupted, the kinase is opened to an active state and the autophosphorylation of the residue Tyr416 occurs. Protein interactions also act to regulate c-Src by either directly activating c-Src, or by moving c-Src to sites of action16. As such, it results that the complex network of interactions determining the structural behavior of the c-Src suggests a conformational plasticity intrinsic to this protein family. On the other hand, ATP analogs such as tyrophostin and pyrimidine compounds, directly inhibit the tyrosine kinase activity of c-Src and/or related kinases17,18,19,20,21. However, the inhibition of c-Src by means of small molecules remains quite difficult due to the very likely to generate mechanisms of resistance such as point mutations in the drug binding-site21. Therefore, alternative routes to develop selective drugs represent an invaluable resource. Among others, the identification of allosteric binding-sites able to affect the c-Src conformational behavior represent an effective strategy able in principle to overcome the drug-resistance and selectivity issues22. The effectiveness of such a strategy strongly depends on the knowledge of the structural modifications occurring between the inactive and the active state of the kinase. In this way, the complex regulatory network - which dynamically evolves along the transition pathway - could be effectively handled to regulate the protein activity. From a computational point of view, the description of the pathway – although possible in principle – is strongly limited by the process timescale, preventing the effectiveness of computational techniques such as molecular dynamics23. To overcome such a limit, several methods able to enhance the sampling have been applied to study the kinase conformational transitions, such as coarse-grained representations and enhanced sampling techniques24,25,26.


5


Page 6 of 30

Recently, an unprecedented computational effort based on distributed computing has been used to describe the c-Src activation pathways by means of all-atoms unbiased molecular dynamics simulations27,28. This work has not only been able to map the c-Src activation conformational transition without the use of external biases, but it provides an invaluable reference to measure the effectiveness and the accuracy of different enhancing sampling techniques. In this framework, we present here the study of the c-Src conformational transition pathways by means of the Essential Dynamics Sampling, an enhanced sampling technique able to efficiently sample such slow and wide conformational rearrangements. The Essential Dynamics Sampling (EDS) reject movements (in the conformational space) not in the desired direction, driving the sampling towards a known structure (target). The advantages of the EDS are that i) no ad hoc terms are introduced in the Hamiltonian and ii) the system moves along its own essential eigenvectors as provided by unbiased molecular dynamics (MD). Moreover, the contraction procedure used in this work guarantees – by construction – that the system reach the final state, being only slightly constrained all along the pathway29. With respect to other theoretical methods potentially able to model these kind of processes (i.e. free-energy methods or coarse-grained models), it is worth noting that in the EDS procedure, i) the model is not biased and, thus it is easily transferable to other systems/processes; ii) the model does not require ad-hoc parametrization and iii) the description is at atomistic level. The main limitation of the EDS is the difficulty to estimate the relative weight of the intermediate conformations sampled along the path. In this work, the EDS technique has been applied to the both opening and closing process of cSrc to provide a detailed atomistic description of the corresponding pathways, unrevealing for


6

Page 7 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the first time the differences at the base of such complex conformational transitions. The EDS trajectories have been analyzed in detail and compared with the reference trajectories as obtained by means of massive distributed computing27. The results clearly show that this kind of processes can be described by the proposed procedure with very limited computational efforts and with an accuracy in the pathway description comparable to unbiased molecular dynamics simulations. 2. COMPUTATIONAL METHOD Molecular dynamic simulations were performed on c-Src using Gromacs software package version 5.0.5.30. The initial structures were the c-Src in the closed (inactive state, PDB ID: 2src) and open conformation (active state, PDB ID: 1y57), respectively. The kinase domains only were considered: residues from Ala175 to Leu449 for 1y57, and Ala259 to Leu533 for 2src. The two conformational state of the kinase domains were (separately) solvated in rectangular boxes large enough to contain the enzyme and 1 nm of the solvent on all sides, resulting in 84514 atoms and 72049 atoms for the open and closed form, respectively. Sodium ions were added to provide a neutral simulation cell. During the simulations the temperature were maintained at 300K using the velocity rescaling algorithm31. The Amber99sb force field was used32. For the solvent, the water SPC model33 was used at standard density (1 kg/dm3). The Verlet cut-off scheme was used34; long range electrostatic interactions were treated by means of the particle-mesh Ewald method35,36. For all simulations the protein secondary structures were monitored to ensure that no unfolding occurred. All the analysis were performed using the Gromacs tools included in the package, unless otherwise specified. The first 20 ns of the two trajectories were considered as equilibration and they have been removed from the analysis. In line with our previous work37, the Principal Component Analysis (PCA) was performed using the covariance matrix of the fluctuations of the C-alpha motions as obtained by the combined trajectories of both the


7


Page 8 of 30

conformational states. In this way, the eigenvectors describe the motion directions connecting the close and the open form of the c-Src, and they can be used in the EDS procedure to decrease the distance from the reference conformation (i.e. contraction procedure). In such a procedure, at each time frame a classic MD step is performed and the distance between the current and the reference conformation in the subspace is calculated. The step is accepted if the distance decreases. If this does not occur, the atomic coordinates are projected radially on the subspace hypersphere centered on the target conformation. The hypersphere is defined by the set of the essential eigenvectors of the C-α atoms given by the PCA analysis performed on the two concatenated (free) MD trajectories. The displacement in the essential subspace during one time step of molecular dynamics can be represented by: ∆ = ∆ + ∆ where ∆ is the displacement produced by the MD step and ∆ is the correction for the application of the constraint. The constraint in the essential subspace can be defined as: (( ; = 0 To obtain a unique solution to this equation, the requirement that the total perturbation |∆ | is minimized adding one Lagrangian multiplier (the mathematical details can be found in ref.29). Following our previous work37, we performed a preliminary set of EDS simulations varying the number of eigenvectors used to define the subspace. In this way, we assess the minimum number of eigenvectors able to describe the transitions but still leading final structures very similar to those extracted from the free MD simulations and used as target. We found that the use of the first 100 eigenvectors leads to structures which are still very similar to the corresponding targets


8

Page 9 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(see Table S1 in S.I.), providing at the same time a very limited bias on the trajectory. To select proper reference structures representative of the two conformational states, a cluster analysis of the unbiased MD simulations has been performed (the Gromacs tool g_cluster with the default cut-off value of 0.1 nm on the RMSD was used), providing 9 representative structures for both the open and the closed form. Every possible structure pairs have been then used as starting and arrival state in the EDS procedure, providing 162 (81x2) possible paths connecting the two protein conformational states. A scheme of our approach is reported in S.I. (Fig. S1). These conformational pathways have been compared to those obtained by Shukla et al.27. In that study, a massive distributed computing (m.d.c.) approach was used to reproduce the transition between the closed and open form of c-Src. To this end the authors employed a very large number of small trajectories which required a huge computational effort. In that work, the authors collected about 500 ms of all-atoms unbiased MD trajectories and ~ 20.000 structures have been considered as representative of the whole set of MD simulations (freely available from https://searchworks.stanford.edu/view/cm993jk8755).

This

representative

trajectory

was

converted in “dcd” format to import and analyze it by means of the “bio3d” package of R software. To compare our approach with the m.d.c. results, we calculated for each c-Src structure, as sampled by both m.d.c and EDS approach, two structural observables (defined in Shukla et al.27 and able to describe the conformational pathways - were calculated: i) the Root Mean Square Deviation (RMSD) of the activation loop (A-loop) and ii) the distance difference between two pairs of residues (i.e. Glu310 and Lys295). The RMSD was calculated using the “rmsd” function of “bio3d” package of R software.


9


Page 10 of 30

The native contacts analysis has been performed considering the distance between all the Calpha atoms of the protein. For any frame obtained by the EDS, the c-alpha-c-alpha distance is measured and the residues pair is considered in contact if their distance is less than 6 Å. Two residues are not considered in contact if they have a sequence distance of less than 4, even if their distance in structure is less than 6 Å. The first frame of each EDS trajectory (which represent the protein in the starting conformation) is considered as the reference structure for the analysis. For each frame we defined the fraction of native contacts (Q) as the ratio between the number of contacts in common with the reference structure and the number of contacts for such a frame. 3. RESULTS & DISCUSSION To obtain an accurate description of the c-Src kinase domain in the two reference states (i.e. open and closed), we firstly performed a 200 ns long classical molecular dynamics simulation for both the active and the inactive form of c-Src. As expected, the secondary structure elements as well the overall fold do not change along the trajectories and both the proteins remains in the conformational basin common to the corresponding the starting structure. Following our previous work37, we apply the EDS technique using the set of the first 100 eigenvectors and eigenvalues obtained by the diagonalization of the atomic fluctuations covariance matrix of the trajectory obtained merging the active and the inactive molecular dynamics simulations (hereafter called the concatenated trajectory).The resulting first two eigenvectors obtained by the Essential Dynamics (ED) analysis on the concatenated trajectory represent the 92% of the variance of the whole protein motion, which includes the open and the closed conformational basins. As expected, the first eigenvector - characterizing the motion direction linking the two


10

Page 11 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


conformational states -describes almost the total motion of the system (~ 87%). By means of a clustering analysis 18 representative structures (9 for the “active” and 9 for the “inactive” form) have been selected from the two unbiased MD and used as starting and target structures for the EDS procedure. That is, the EDS runs were performed to connect each of the 9 representative structures of the open state with each one of the 9 structures of the closed states (and vice versa). In this work, the EDS algorithm have been used in the contraction scheme, i.e. to decrease along the simulation - the distance in the essential subspace between the selected reference structure and the target structure. By such an approach, 81 trajectories from close to open conformations and 81 trajectories from open to close conformations have been obtained.

Activation pathway identified through EDS technique The analysis of the EDS trajectories performed along the activation pathway, reveal that the main changes are carried – as expected - by the A-loop and the C-helix (see Fig. 2). In particular, in the opening process, the rearrangements of such regions occur in a concerted fashion: the A-loop runs outwards while the C-helix moves inwards, in completely agreement with previous computational studies13,38. Other conformational changes involve residues buried within the protein domain, Leu325, Met314, Phe405 and His384. They form a skeleton of four nonconsecutive hydrophobic residues that constitute a regulatory region, called R-spine, linking the two lobes of the catalytic domain. It is critical for the catalytic activity of the kinase, because it is usually broken in the deactivation process39,40. Our analysis confirms that this interaction pathway is rapidly broken during the transition from the active to the inactive form. Our results also show that His384, Arg385 and Asp386 (called HRD motif) maintain the same orientation during the opening process. In the


11


Page 12 of 30

open conformation, the HRD motif does not interact with the A-loop nor with the C-helix. Only after that the C-helix have reached the open conformation and the A-loop moved outwards, Arg385 starts to interact with the Tyr416 belonging to the A-loop. The electrostatic interactions between Arg385, Tyr416, Arg409, Glu310, Lys295 and Asp404 are known to be involved in the conformational changes occurring during the activation and the deactivation processes39. Previous studies have shown that in the closed conformation there is a peculiar electrostatic interaction network formed by the following couples: Glu310 with Arg409, Arg385 with Tyr416 and Lys295 with Asp404. In the open conformation these interactions do not take place. Our simulations, confirming such a behavior (Fig. 3 A and B), also outline that in the activation process the first interaction to be broken is between Arg385 with Tyr416, anchoring the A-loop to the HRD motif (Tyr416 moves away so that the A-loop starts to move outwards and, almost simultaneously, the C-helix starts moving inwards). As a result of these motions, the residues Arg409 (belonging to the A-loop), and Glu310 (in the C-helix), no longer interact (Fig. 3 A). Deactivation pathway identified through EDS technique In the deactivation process we found that the HRD motif interacts with the Tyr416. As the A-loop moves inwards, Tyr416 (belonging to the A-loop) starts to interact with R385 as observed in a previous study of Src kinase conformational activation41. During the deactivation process, our simulations show that the interaction between Tyr416 and Arg409 early disappears; the A-loop starts to fold, and the C-helix is almost in the closed conformation. However, the interaction between Lys295 and Glu310 is still in place. At the same time, the interaction between Lys295 and Glu310 vanishes. The C-helix is in the closed conformation and starts to twist and tilt, driving the Glu310 outwards where it finally interacts


12

Page 13 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with Arg409. At the end of the deactivation process, the interactions between Arg385 and Tyr416 take place, stabilizing the structure in the closed conformation (Fig. 3 B). The closed conformation is stabilized by the interaction between the A-loop with the HRD motif, the C-helix interacting with the N-lobe. After the Lys295 and Glu310 interaction is broken, the A-loop and C-helix show a concerted movement until the interactions stabilizing the closed conformation take place. Activation pathway: Comparison between EDS and m.d.c approach As highlighted in the Introduction section, the main goal of the present work has been the validation of the sampling as obtained by the EDS technique. To this end, the EDS trajectories were compared to the recent work of Shukla et al.27 where the conformational transition pathways of the c-Src activation process were described by massive distributed computing (m.d.c.). In that work, making use of huge computational resources, the authors performed ~ 24000 MD simulations totally lasting 50 microseconds. The results were summarized in a single trajectory (using an adaptive sampling algorithm based on the Markov State Model approach), which we used to compare our EDS data. In the work of Shukla et al.27 the authors characterized two intermediate states of c-Src kinase domain by means of two simple structural observables: i) the RMSD of the A-loop with respect to the inactive state and ii) the distance between two residues pairs (Glu310-Arg409 and Lys295Glu310). The comparison between the behavior of these two observables in the two different approaches is reported in Fig. 4. We can clearly notice that the EDS essentially provides the same intermediate structures of the unbiased simulations with respect to these observables (RMSD of the A-loop


13


Page 14 of 30

and the difference between Glu310-Arg409 and Lys295-Glu310 distances). In fact, starting from similar (i.e. sampled within the same conformational basin) closed structures both the methods describe the activation process by an initial change in the A-loop conformation followed by the formation of the Lys295-Glu310 H-bond and the rupture of the Glu310-Arg409 one. The differences in the point densities is due to the fact the EDS algorithm does not provide the correct weight of the sampling as the m.d.c does. Very importantly, such results strongly suggest that at orders of magnitude faster - the EDS algorithm actually provides the same intermediates as those sampled by unbiased massive distributed computing approach (at least in the subspace described by the A-loop RMSD and by the difference between Glu310-Arg409 and Lys295Glu310 distances). To further compare the two sampling methodologies, the EDS and the m.d.c., we projected the two sets (i.e. EDS and m.d.c.) of structures sampled along the activation pathway on the subspace spanned by the two principal eigenvectors (Fig. 4). As the eigenvalues correspond to the mean square eigenvector coordinate fluctuation, and, therefore, contain the contribution of each principal component to the total fluctuation, the projections have been scaled according to the ratio between the two eigenvalues (0.057). The projections along the first principal component as obtained by the EDS technique reported in Fig. 5 closely resembles those provided by the m.d.c. approach. The discrepancy along the second principal component, between the EDS and m.d.c. inactive state (still very small), is due to different condition used to define the productive phase of the corresponding simulations. In fact, Shukla et al. have reconstructed a unique trajectory to describe the conformational transitions of Src kinase combining several short trajectories – lasting 20 ns plus 1 ns considered as equilibration time – into a single model27. On the other hand, the 9 reference structures used as


14

Page 15 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


starting point and the 9 structures representing the target (arrival) points of the EDS trajectories were extracted from the corresponding conformational basins as provided by two single and independent “free” 200 ns MD runs (note that these 200 ns long MD simulations were only used to describe the two forms, whereas the paths connecting the two forms have been modelled by means of EDS). In such two MD simulations, the first 20 ns were considered as an equilibration phase and, hence discarded from the analysis. Including such 20 ns of sampling (Fig. S2), it is clear that the differences between our approach and the m.d.c. one is due to such a difference in the starting points of the pathway. In summary, the presented data clearly provide evidences that the EDS trajectories describe conformational changes on the activation pathway very similar to those obtained by means of the m.d.c. approach. There results strongly suggest that our methodology can provide a reliable sampling along the opening and closing process occurring in the c-Src kinase. Activation and deactivation pathway Due to the EDS computational efficiency, the EDS procedure has also been applied to model the deactivation process. The major structural observables changing along the process has been highlighted in the previous section. Nevertheless, additional information arises when comparing the opening and the closing pathways. First, the projection of the closing pathways on the essential subspace indicates no relevant differences with respect to the c-Src activation within such a subspace (Fig. 6). On the contrary, the plot of the distances between the Glu310, Arg409 and Lys295 versus the RMSD of the A-loop as occurring along the open-to-close pathway (Fig. 7) shows some interesting differences with respect to the activation process (Fig. 4). In fact, whereas in the opening process the H-bond network changes follow the A-loop “unfolding”


15


Page 16 of 30

process (Fig. 5), in the closing process the H-bond network starts to change quite quickly (Fig. 7). In other words, our data suggest that, differently from the activation pathway, the open-to-close mechanism might involve the rupture of the H-bond between Lys295 and Glu310 in the very first part of the process. Differences between the two pathways are also confirmed by the C-helix tilt angle. Considering that the C-helix movement is one of the major conformational changes involved in the processes, in Fig. 8 we report the probability to have a tilt angle variation greater than 60° along the paths. That is, we consider that the A-helix conformational change occurred when the tilt angle (variation) reaches this value, taking as zero its initial value (for sake of clarity graphical representations of c-Src structures at various tilt angle values are given in Fig. S3). As shown in the Fig. S3, the C-helix transition occurs in a wide interval, indicating that different EDS trajectories employs a different number of EDS steps to move the C-helix in the final inactive state. On the contrary, in almost all the EDS opening trajectories, the C-helix suddenly reaches its final open conformation. In Fig. S4, we also reported the fraction of conserved native contacts (Q) trajectories along the conformational transitions for 6 representative EDS trajectories (three for the activation and three for the deactivation process). Interestingly, both the processes retain a remarkable fraction (between 70% and 80%) of the native contacts of the corresponding starting forms, indicating that a large part of local contacts is not influenced by such large conformational rearrangement. As expected, the protein regions mainly determining such native contact decrease are the C-helix (residues 298-309) and the A-loop (res 411-419) ones (see Fig. S4).


16

Page 17 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


These results clearly point out different mechanisms between the activation and the deactivation processes: in the former, the H-bond network is maintained (Fig. 4) and the C-helix transition can occur on a wide interval of frames (Fig. 8, blue line), whereas in the deactivation mechanism, the H-bond (Fig. 7) and C-helix changes (Fig. 8, green line) are likely to occur in the very first part of the process. The described mechanism of activation/inactivation has the potential to improve the drug selectivity (which is very limited in the case of conserved binding site), for example designing new molecules targeting a low-conserved protein region which might affect one specific protein conformational transition pathway). Therefore, the proposed approach, able to describe slow processes occurring in large biomolecular systems with extreme computational efficiency, seems very well suited for the design of new inhibitors, acting on different protein regions as well as to shed new lights on complex conformational equilibria. In addition, it can give insights into the conformational rearrangements occurring in different processes such as small ligand binding or transition channel opening. 4. CONCLUSIONS In this study we modeled the conformational transitions occurring along the pathways connecting the active and inactive states of c-Src kinase. The complex molecular rearrangements occurring along the activation pathways are in line with previous observations based on unbiased all-atoms simulations obtained by means of massive distributed computing. The active to inactive transition is mainly described by the C-helix and A-loop rearrangements as well as the dynamical change of the interactions in important c-Src motifs (R-spine and HRD motif). Our results clearly show that our approach provides a similar accuracy with respect to


17


Page 18 of 30

unbiased molecular dynamics simulations at orders of magnitude faster and, very importantly, without arbitrariness on the choice of the parameters needed to describe the process. Due to the EDS computational efficiency and the agreement of our data with that produced by Shukla et al.

27

for the c-Src activation pathway, the deactivation pathway has been also

described. Interestingly, the associated mechanism is slightly different with respect to the activation: we have found that - in contrast to the opening pathway - the C-helix rearrangement occurs in the very first few EDS steps during the deactivation and the characteristic H-bond network of the closed structure involving the residues Lys295, Glu310 and Arg409. Finally, we think that the possibility to describe the molecular pathways connecting the active and inactive states at reduced computational cost might improve the discovery of new drugs able to interact with different protein regions and traps the protein in a specific conformational state as well as to better understand the effect of point mutations on complex conformational equilibria.

ASSOCIATED CONTENT Supporting Information. Table S1, Figure S1, S2, S3 and S4 are supplied in the Supporting Information.

AUTHOR INFORMATION Corresponding Author *Email: [email protected]


18

Page 19 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Aknowledgement The authors acknowledge Sapienza University of Rome for financial support and CINECA and NVIDIA for the computational support.


19


Page 20 of 30

REFERENCES 1. Warmuth, M.; Damoiseaux, R.; Liu, Y.; Fabbro, D.; Gray, N. SRC family kinases: potential targets for the treatment of human cancer and leukemia. Curr Pharm Des. 2003, 9(25), 2043-2059. 2. Guarino, M. Src signaling in cancer invasion. J Cell Physiol. 2010, 223(1), 14-26. 3. Irtegun, S.; Akcora-Yıldız, D.; Pektanc, G.; Karabulut, C. Deregulation of c-Src tyrosine kinase and its downstream targets in pre-eclamptic placenta. J Obstet Gynaecol Res. 2017, 43(8), 1278-1284. 4. Benati, D.; Baldari, C. T. SRC family kinases as potential therapeutic targets for malignancies and immunological disorders. Curr Med Chem. 2008, 15(12), 1154-1165. 5. Irtegun, S.; Wood, R. J.; Ormsby, A. R.; Mulhern, T. D.; Hatters, D. M. Tyrosine 416 is phosphorylated in the closed, repressed conformation of c-Src. PLoS One 2013, 8(7), e71035. 6. Patwardhan, P.; Miller; W. T. Processive phosphorylation: mechanism and biological importance. Cell Signal. 2007, 19(11), 2218-2226. 7. Thomas, S. M.; Brugge, J. S. Cellular functions regulated by Src family kinases. Annu Rev Cell Dev Biol. 1997, 13, 513-609. 8. Wheeler, D. L.; Iida, M.; Dunn, E. F. The role of Src in solid tumors. Oncologist. 2009, 14(7), 667-678. 9. Dehm, S. M.; Bonham, K. SRC gene expression in human cancer: The role of transcriptional activation. Biochem Cell Biol. 2004, 82(2), 263-274. 10. Roskoski, R., Jr. Src protein–tyrosine kinase structure and regulation. Biochem Biophys Res Commun. 2004, 324(4), 1155-1164. 11. Ahmad, S.; Banville, D.; Zhao, Z.; Fischer, E. H.; Shen, S. H. A widely expressed human protein-tyrosine phosphatase containing src homology 2 domains. Proc Natl Acad Sci U S A. 1993, 90(6), 2197-2201. 12. Xu, W.; Doshi, A.; Lei, M.; Eck, M. J.; Harrison, S. C. Crystal structures of c-Src reveal features of its autoinhibitory mechanism. Mol Cell. 1999, 3(5), 629-638.


20

Page 21 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


13. Xu, W.; Harrison, S. C.; Eck, M. J. Three-dimensional structure of the tyrosine kinase c-Src. Nature. 1997, 385(6617), 595-602. 14. Martin, G. S. The road to Src. Oncogene. 2004, 23(48), 7910-7917. 15. Yamaguchi, H.; Hendrickson, W. A. Structural basis for activation of human lymphocyte kinase Lck upon tyrosine phosphorylation. Nature. 1996, 384(6608), 484-489. 16. Kaplan, K. B.; Bibbins, K. B.; Swedlow, J. R.; Arnaud, M.; Morgan, D. O.; Varmus, H. E. Association of the amino-terminal half of c-Src with focal adhesions alters their properties and is regulated by phosphorylation of tyrosine 527. EMBO J. 1994, 13(20), 4745-4756. 17. Altmann, E.; Widler, L.; Missbach, M. N(7)-substituted-5-aryl-pyrrolo[2,3-d]pyrimidines represent a versatile class of potent inhibitors of the tyrosine kinase c-Src. Mini Rev Med Chem. 2002, 2(3), 201-208. 18. Druker, B. J.; Lydon, N. B. Lessons learned from the development of an abl tyrosine kinase inhibitor for chronic myelogenous leukemia. J Clin Invest. 2000, 105(1), 3-7. 19. Wakeling, A. E.; Guy, S. P.; Woodburn, J. R.; Ashton, S. E.; Curry, B. J.; Barker, A. J.; Gibson, K. H. ZD1839 (Iressa): An orally active inhibitor of epidermal growth factor signaling with potential for cancer therapy. Cancer Res. 2002, 62(20), 5749-5754. 20. Lynch, T. J.; Bell, D. W.; Sordella, R.; Gurubhagavatula, S.; Okimoto, R. A.; Brannigan, B. W.; Harris, P. L.; Haserlat, S. M.; Supko, J. G.; Haluska, F. G.; et al. Activating mutations in the epidermal growth factor receptor underlying responsiveness of non-small-cell lung cancer to gefitinib. N Engl J Med. 2004, 350(21), 2129-2139. 21. Roskoski, R. Src protein-tyrosine kinase structure, mechanism, and small molecule inhibitors. Pharmacol Res. 2015, 94, 9-25. 22. Zhang, S.; Yu, D. Targeting Src family kinases in anti-cancer therapies: Turning promise into triumph. Trends Pharmacol Sci. 2012, 33(3), 122-128. 23. D'Abramo, M.; Besker, N.; Chillemi, G.; Grottesi; A. Modeling conformational transitions in kinases by molecular dynamics simulations: Achievements, difficulties, and open challenges. Front Genet. 2014, 5, 128.


21


Page 22 of 30

24. Huber, T.; Torda, A. E.; Gunsteren, W. F. Local elevation: A method for improving the searchig properties of molecular dynamics simulation. J Comput Aided Mol Des. 1994, 8(6), 695-708. 25. Grubmüller, H. Predicting slow structural transitions in macromolecular systems: Conformational flooding. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics. 1995, 52(3), 2893-2906. 26. Tozzini, V. Coarse-grained models for proteins. Curr Opin Struct Biol. 2005, 15(2), 144-150. 27. Shukla, D.; Meng, Y.; Roux, B.; Pande, V. S. Activation pathway of Src kinase reveals intermediate states as targets for drug design. Nat Commun. 2014, 5, 3397. 28. Meng, Y.; Shukla, D.; Pande, V. S.; Roux, B. Transition path theory analysis of c-Src kinase activation. Proc Natl Acad Sci U S A. 2016, 113(33), 9193-9198. 29. Amadei, A.; Linssen, A. B.; de Groot, B. L.; van Aalten, D. M. F.; Berendsen, H. J. C. An efficient method for sampling the essential subspace of proteins. J Biomol Struct Dyn. 1996, 13(4), 615-625. 30. Abraham, M. J.; van der Spoel, D.; Lindahl, E.; Hess, B.; and the GROMACS development team; GROMACS User Manual version 5.0.5; www.gromacs.org, 2015. 31. Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J Chem Phys. 2007, 126(1), 014101. 32. Lindorff‐Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side‐ chain torsion potentials for the Amber ff99SB protein force field. Proteins. 2010, 78(8), 1950-1958. 33. Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans; J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; p. 331. 34. Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. Comput Phys Commun; 1995, 91, 43–56. 35. Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N•log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089–10092.


22

Page 23 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


36. Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh ewald potential. J. Chem. Phys. 1995, 103, 8577–8592. 37. Bešker, N.; Amadei, A.; D'Abramo, M. Molecular mechanisms of activation in CDK2. J Biomol Struct Dyn. 2014, 32(12), 1929-1935. 38. Roskoski, R., Jr. Src protein–tyrosine kinase structure and regulation. Biochem Biophys Res Commun. 2004, 324(4), 1155-1164. 39. Meng, Y.; Roux, B. Locking the active conformation of c-Src kinase through the phosphorylation of the activation loop. J Mol Biol. 2014, 426, 423-435. 40. Ozkirimli, E.; Post, C. B. Src kinase activation: A switched electrostatic network. Protein Sci. 2006, 15(5), 10511062. 41. Yang, S.; Roux, B. Src kinase conformational activation: Thermodynamics, pathways, and mechanisms. PLoS Comput Biol. 2008, 4(3), e1000047.


23

linker


1 C-helix 2 3 4 5 6 7 8 9 10 11 A-loop 12 13 C-helix 14 15 16 17 A-loop 18 19 20 21 22 23 24 25 26 27 Fig. 28 1 29 Comparison of active (PDB id: 1FMK) inactive (PDB id: 2SRC) crystal structures of30the c-Src kinase. Coloring from N to C terminus: SH3 domain in green, SH2 31 domain in blue, linker connecting SH2 domain and N-lobe of the kinase domain in 32 purple, kinase domain N-lobe and C-lobe in light brown. αC helix and A-loop are 33 colored in deep brown and red respectively. 34 35 36 37 38 ACS Paragon Plus Environment 39 40 41

Page 24 of 30

linker

Page 25 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


A-loop moves outward

Fig. 2 Superimposition of active (in blue) and inactive (in gold) c-Src conformations along the activation pathway obtained by using EDS approach. Conformational changes of the αC helix and A-loop (resdiues 402-424) reveal outward movements of the A-loop while the αC-helix moves inwards.



A 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

B

Fig3 A,B. Changes in the position and interaction of specific residues (Arg385, Tyr416, Arg409, Glu310, Lys295 and Asp404) accompanying the large conformational changes during the activation and deactivation processes (A: open form; B: closed form).


Page 26 of 30

Page 27 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


Fig. 4 Two-dimensional conformational landscape of c-src tyrosine kinase with respect to r.m.s.d. of A-loop residues (x-axis) and difference of the distance between Glu310-Arg409 and Lys259-Glu310 (y-axis) residue pairs. Changes in the structural observables provided by the EDS and m.d.c. approaches are reported as blue and red points respectively.

Fig.5 Projection along the first and second principal components of 9 EDS opening trajectories (blue points) connecting the projections given by MD of open and closed states (black points in left and right region respectively). Projections of the closed and intermediate structures as provided by the m.d.c approach are reported as gray and red points.



Page 28 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 Fig. 7 24 Fig. 6 Changes in the structural observables as provided by the 25 Projection along the first and second principal components of 9 26 EDS deactivation trajectories (green points) connecting the EDS during the deactivation process. In x-axis is reported 27 projections given by MD of open and closed states (black points in the RMSD of A-loop, while in the y-axis is reported the 28 left and right region respectively). Projections of the closed and difference between two distances: the first and the second 29 intermediate structures as provided by the m.d.c approach are distance are related to pairs Glu310-Arg409 and Lys29530 Glu310 respectively. 31 reported as gray and red points 32 33 34 35 36 37 38 39 40 41 42 43 44 ACS Paragon Plus Environment 45 46 47

Page 29 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41


Fig. 8 The probability to observe the C-helix tilt transition, as measured by a corresponding tilt angle greater than 60°, versus the number of the EDS frames for opening (blu line and points) and for the closing (green line and points) processes


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

Transition probability


frames ACS Paragon Plus Environment

Page 30 of 30

Efficient and Accurate Modeling of Conformational Transitions in

Recommend Documents