Modelling the Activation Pathways in Full-Length Src Kinase

Src kinases play fundamental roles in several crucial cell processes. Their activity is tightly regulated by conformational transitions between the active and the inactive forms, which are carried out by complex protein structural rearrangements. Here, we present an in-depth study of such structural transitions coupling extensive all-atoms molecular dynamic simulations coupled to an algorithm able to drive the system between defined conformational states. Our results, in line with the available experimental data, confirm the complexity of such a process indicating the main molecular determinants involved. Moreover, the role of an Src inhibitor—able to bind to the protein inactive state—is discussed and compared with available experimental data.


Introduction
c-Src is a tyrosine kinase protein that regulates many intracellular processes such as cell growth, differentiation, and proliferation phosphorylating intracellular proteins. c-Src is expressed in most the human cells, and a dysfunction [1][2][3] in its activation is related to the onset of neoplastic progression [4][5][6][7][8][9]. Such a protein is formed by a myristoylated SH4 domain that anchors the kinase on the membrane surface, the Unique domain, the structural SH3 and SH2 domains and the kinase domain (KD), which contains the active site [10][11][12] (see Figure 1). As other kinase proteins, c-Src lives at least in three conformations [7]: (i) a close-inactive form, stabilized by the phosphorylation of Tyr 527, located in the C-terminus region; (ii) an unphosphorylated primed form; and (iii) an open-active form, characterized by the phosphorylation of Tyr 416 in the kinase domain (KD). The SH2 and SH3 domains participate in the regulation of the conformational equilibrium between the structurally different open and closed forms by intermolecular interactions. The catalytic site [13,14] is located in the cleft between the C-lobe and the N-lobe of the KD. The phosphorylation of Tyr 527, belonging to the C-terminal tail of the KD, results in the KD-SH2 binding that stabilizes, together with additional interactions between the SH3 domain and linker located between the SH2 and the KD domain, the closed form of the protein [15,16]. The dephosphorylation of Tyr 527 triggers a conformational rearrangement leading to the open/active form, where the interaction between the C-terminal tail and the SH2 domain vanished. The successive autophosphorylation of Tyr 416 further stabilizes the active form [2,17]. Therefore, the conformational equilibrium between the active and inactive states is tightly regulated by different molecular interactions that are dependent on different input signals (Figure 1).
To this end, we characterize these interaction patterns and their role in the conformational pathways in the c-Src full construct at the atomic level by means of a computational study performed by means of Molecular Dynamics (MD) simulations and an enhanced sampling technique, the Essential Dynamic sampling (EDS) [18]. By such an approach, we were able to describe c-Src full-length activation/inactivation paths and to characterize the crucial movements of the different protein regions that lead to the conformational transitions.
Due to the fact that the deregulation of protein tyrosine kinase activities leads to several pathologies, such a protein family represents an important target from a medical perspective. Therefore, several inhibitors have been developed so far. Among others, imatinib (gleveec) [19][20][21], belonging to the type II class of inhibitors, was one of the most effective for a closely related kinase, the Abl tyrosine kinase. The differences in the inhibitory power of imatinib with respect to c-Src and Abl tyrosine kinases-which showed a high degree of similiarity in terms of both sequence and structure-stimulated structural investigations to explain this issue [22,23]. From a structural point of view, imatinib binds to the c-Src kinases in the closed conformation, as illustrated by crystallographic studies of the c-Src-imatinib bound structure [24]. In this respect, we study here the effect of imatinib binding to the closed form of the c-Src kinase to highlight its role in conformational transitions. We anticipated that, in line with the literature, such an inhibitor does not affect the conformational behavior and/or transitions, thus supporting previous studies [22,23,25,26] indicating that imatinib's main effect is to stabilize the relative population of the closed form, thus lowering the c-Src activity.

MD Simulations
From the human c-Src crystallographic structures of the closed (PDB ID:2src) and the open (PDB ID:1y57) states, the phosphate groups were removed from the phosphorylation sites to obtain the primed conformations [27][28][29]. The imatinib coordinates were taken from the DRUGBANK [30] database, and through the superimposition on the crystal structure of the chicken c-Src in complex with imatinib (PDB ID: 2oiq), the drug was docked into the closed form of the human c-Src kinase domain. The imatinib topology has been generated using the Antechamber package [31] implemented to be used with GAFF [32], with the RESP charge method [33]. In order to convert the Amber topology in the Gromacs format, the acpype (AnteChamber PYthon Parser interfacE) tool [34] has been used. The systems (c-Src closed state, c-Src open state, and c-Src closed state bound to imatinib) were solvated with the TIP3P water model [35] in a cubic box, imposing a minimum distance of 1.5 nm between the protein and the box. Sodium ions were added to neutralize the systems. The particle mesh Ewald method [36], with a cut-off of 1.2 nm, was adopted to treat the electrostatic interactions, and the Verlet cut-off scheme was used for the van der Waals interactions (with a cut-off of 1.2 nm). The MD simulations were performed in the NVT ensemble using the velocity rescale algorithm [37], imposing a τ = 0.1 and a temperature of 300 K. The Gromacs 2018 Software [38] was used to carry out the simulations using the Amber99sb force field. For the apo states, seven simulations of the open structures, each lasting ≈600 ns, and two simulations of the closed forms of ≈800 ns each were performed. The holo state was simulated for 500 ns for both the configurations.

Principal Component Analysis (PCA) and Essential Dynamics Sampling (EDS)
Principal component analysis (PCA) [39] was performed in order to extract the more relevant motions from the trajectory. In particular, the diagonalization of the covariance matrix of the C-α fluctuations, provides a set of eigenvectors that are ordered according to the corresponding eigenvalues. The eigenvectors with the largest eigenvalues describe the largest-amplitude collective motions of the protein. PCA was performed on the combined trajectories of both the open and close conformational states. In this way, the eigenvectors describe the main motion directions connecting the close and the open forms of the protein, and they can be used in the Essential Dynamics Sampling (EDS) procedure [18] to model the transition pathways between the two conformational states, i.e., the closed and the open states. The EDS technique is based on the idea of driving the sampling towards a target structure rejecting modes that are not in the desired direction. In this procedure, two unbiased simulations of the starting and final states, describing the two equilibrium basins are required together with a certain number of short simulations that drive the system to the final state. For each MD step, if the distance between the intermediate and the reference configuration decreases in the selected subspace, the step is accepted; otherwise, a correction is applied on the atomic coordinates projected on the essential subspace. The mathematical details of such a procedure are reported in previous works [18,39]. Following this protocol, the possible transition paths that connect the closed and open protein conformations as well as the structural changes involved in the activation/deactivation mechanism are identified.

Conformational Equilibria
The analysis of the MD simulations of the closed state indicates that the sampled configuration is structurally similar to the starting crystallographic structure (Figure 2).
On the other hand, the trajectories of the open (active) state show a pronounced conformational variability (Figure 2), mainly due to mutual orientation of the three protein domains-SH3, SH2, and KD-occurring along the MD trajectories. The fluctuation analysis indicates that they are more pronounced in the SH2 and KD and in the linker region connecting these two domains (data not shown) as also confirmed by hydrogen bond analysis (see Supplementary Materials). This large conformational variability is even more evident using principal component analysis (PCA). In fact, combining open and closed trajectories and applying the PCA , the protein projection on the essential plane defined by the first two eigenvectors clearly shows that the space explored by the open states is much wider than the closed state MDs (Figure 3). The eigenvalue associated with the first eigenvector-which clearly discriminates between the two states-is very high, describing the 87% of the total variance, thus indicating that the two states sample very different region of the common essential conformational space.  Interestingly, the same analysis applied to the SH3 and SH2 domains shows limited differences, indicating that the most relevant protein motions are due to the mutual domain rearrangements and/or to the KD domains (Figure 4a). Indeed, the ED projections of the KD domains reveal that they are confined in conformational regions near those explored by the corresponding starting structures (Figure 4b).

Activation and Inactivation Pathways
Following our previous approach, where the conformational pathways connecting the open and the closed states of Src protein kinase domains have been described [27,28], here, we apply the same procedure to the full-length c-Src construct. The EDS trajectories describing the activation and the inactivation pathways are represented in Figure 5. Interestingly, the two pathways are quite similar in the essential subspace described by the first two eigenvectors. To better describe these pathways, we monitored some structural observables such as the angle between the three protein domains θ K23 (defined by their center of mass) and the distance d DR between the SH2 and KD domains, computed considering two specific residues: the Asp117 located in the SH2 domain and the Arg318 in the KD. These residues were selected because they form an hydrogen bond that stabilizes the closed conformation and such an interaction is disrupted in the open state (see Figure 6).    These results point out that, during the activation process, the mutual orientation between the three domains starts to change without affecting their inter-domain contacts. Once the domains begin to separate, these interactions disappear, allowing the protein to reach the open/active state. Similarly, in the inactivation process, the formation of favorable contacts between the domains is already observed when the θ K23 values are outside of the range of the closed form. In summary, a certain degree of conformational variability in the mutual domain arrangements is observed even in the presence of specific local contacts between the protein domains.

Src-Imatinib
To characterize the conformational behaviour of c-Src bound to imatinib, a MD simulation of 500 ns of such a system was performed (Figure 9). The rationale here is to investigate how the imatinib binding affects the conformational behavior of specific regions of c-Src. In particular, it is known that imatinib does not only affect the binding site, but it also changes the dynamics of distant protein regions, suggesting an allosteric regulation [40]. To this end, we compared the atomic fluctuations of the apo and holo forms of the c-Src as provided by the corresponding MD simulations (see Figure 10). Very interestingly, a remarkable increase of the fluctuations is observed in the Activation loop (A-loop) well matching the hydrogen-deuterium exchange experiments recently reported [40]. In line with that experimental paper, the C-helix as well as the regions connecting the three domains are all influenced by the presence of imatinib, as it is shown by the corresponding RMSF increases. It is worth mentioning that, by means of free MD simulations, we observed here also a significant growth of the fluctuations in protein regions far from the inhibitor, thus corroborating the idea that imatinib acts on the structural dynamics of the whole protein construct, stabilizing the inactive form [24]. The other possible inhibition mechanism, i.e., that such an inhibitor affects the conformational pathways, is ruled out by our EDS data (see below). When the collective c-Src conformational behavior is investigated by means of principal component analysis, we observe the effect of the inhibitor shifting the protein in a different region of the conformational space ( Figure 11, panel a). Interestingly, imatinib leaves the conformational dynamics of the SH2 and the SH3 domains almost unaltered, but it affects the kinase domain and the linker region motions ( Figure 11, panels b, c, and d). As highlighted by the fluctuations analysis, the inhibitor increases the protein flexibility as indicated by the wider regions of the essential subspace sampled by the c-Src bound form with respect to the unbound one for all of the protein regions investigated (Figure 11).

Activation/Inactivation Pathways of c-Src:Imatinib
To describe the effect of the inhibitor on the activation/inactivation pathways, 322 EDS trajectories for both the activation and the deactivation processes were calculated ( Figure 12). Interestingly, the EDS trajectories of the bound form show almost identical pathways for the activation and the inactivation processes when compared to the corresponding pathways of the unbound forms, at least in the essential subspace.
When these pathways are compared through the structural descriptors θ K23 and d DR used in the characterization of the pathways in the unbound form (see above), a similar behavior is observed. This suggests that the presence of the inhibitor does not influence the evolution of these descriptors: the activation paths reported in Figure 13 show a sudden increase in the θ K23 , the angle between the three domains and, then the successive disappearance of the interaction between the Asp117 located in the SH2 domain and the Arg318 in the KD as described by the increase in the d DR values. Similarly, the inactivation pathways in the presence of the imatinib are described by the concurrent variation in θ K23 and d DR up to values compatible with the inactivate form, followed by the limited decrease in the θ K23 value ( Figure 13). In other words, the effect of the inhibitor, although clear by the alteration in both the fluctuations and the principal component analysis, is almost absent in the conformational changes describing the transition pathways. These results strongly support the picturecorroborated by experimental evidences [24]-that such an inhibitor affects the turnover of kinase, shifting the pool of conformations towards the inactive states.
Videos showing the conformational transitions as obtained by means of EDS approach in both the holo-and apo-states are available in the Supplementary Materials.

Conclusions
Here, we reported the description of the conformational behavior of full-length c-Src kinase in both the unbound and imatinib-bound forms as obtained by means of molecular dynamics simulations. As expected, the conformational variability of the open/active forms is more pronounced compared to the unbound ones. This effect is mainly due to both the mutual domain orientations and to the kinase domain dynamics as revealed by principal component analysis, which also show how the SH2 and the SH3 domains are quite unaffected by the protein conformation. Very interestingly, the effect of the inhibitor in the c-Src closed state is the movement of the essential conformational space in a different region as well as an increase in the fluctuations of some particular regions of the proteins well matching the experimental data based on hydrogen-deuterium exchange experiments. It is worth noting that the presence of the inhibitor increases the fluctuations of regions that are quite far from the binding pocket, thus confirming an allosteric effect. These unbiased molecular dynamics trajectories were then used to apply the essential dynamics sampling procedure to describe the pathways connecting two structurally different protein forms. In line with experimental evidences [24], the corresponding pathways seem quite unaffected by the presence of the inhibitor at both the global and local levels In summary, we have shown how it is possible-by means of extended molecular dynamics simulations and an enhanced sampling algorithm-to describe the conformational equilibria of an important protein member of the tyrosine kinase family and to accurately model the conformational rearrangements occurring during the activation and inactivation processes also in the presence of an inhibitor. The agreement with the available experimental evidences strongly indicates the reliability of such an approach, suggesting its applicability to a wide range of complex molecular systems.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/biophysica1020018/s1, Figure S1: HB between the residues 179T-527Y (black box) and 160R-365D (red box) stabilizing c-Src in close conformation, Figure S2: HB between the residues 11R-221E belonging, respectively, to the SH3 and KD domains for c-Src in open conformation, Table  S1: Hydrogen bonds (HB) in c-Src domains for the closed Apo conformation, Table S2 Funding: This research received no external funding Data Availability Statement: The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.