Using Atomistic Simulations to Explore the Role of Methylation and ATP in Chemotaxis Signal Transduction.

A bacterial chemotaxis mechanism is activated when nutrients bind to surface receptors. The sequence of intra- and interprotein events in this signal cascade from the receptors to the eventual molecular motors has been clearly identified. However, the atomistic details remain elusive, as in general may be expected of intraprotein signal transduction pathways, especially when fibrillar proteins are involved. We performed atomistic calculations of the methyl accepting chemoprotein (MCP)-CheA-CheW multidomain complex from Escherichia coli, simulating the methylated and unmethylated conditions in the chemoreceptors and the ATP-bound and apo conditions of the CheA. Our results indicate that these atomistic simulations, especially with one of the two force fields we tried, capture several relevant features of the downstream effects, such as the methylation favoring an intermediate structure that is more toward a dipped state and increases the chance of ATP hydrolysis. The results thus suggest the sensitivity of the model to reflect the nutrient signal response, a nontrivial validation considering the complexity of the system, encouraging even more detailed studies on the thermodynamic quantification of the effects and the identification of the signaling networks.

Introduction

Bacteria swim in search of nutrients, also known as the chemotaxis movement, displaying what appears to be a biased random walk. The frequencies of the long strides (runs) and the sudden changes of direction (tumbles) that define this random walk are influenced by the nutrient concentration gradients. The runs and tumbles in this random walk are influenced by the anticlockwise or clockwise rotation of the flagella.[1,2] The surprisingly complex sequence of events from sensing to actuation of the flagella all occurs within a single cell, mediated through the physical and chemical modifications in the different groups of proteins that are involved.[2,3]

A schematic of the chemotaxis signaling system is shown in Figure : stimulus sensing, by the ligand (nutrient) binding domain and the trans-membrane (TM) helices;[4] input output control, conveying the signals from the periplasmic region to the cytoplasmic region via the HAMP domain;[5] kinase control, which occurs through a complex of methyl accepting proteins (MCPs),[6] CheA and CheW, that interact with CheB and CheR, of which the histidine autokinase (CheA) is responsible for the ATP hydrolysis for further downstream signaling;[7] and the motor control, involving the CheY and CheZ proteins, which are mainly responsible for taking signals from CheA to the flagella motor.[8,9]

Figure 1

System and model. (a) The core signal transduction system in chemotaxis involves several proteins that are shown schematically. (b) The model used for the molecular dynamics study, based on the structure from Cassidy et al.,[12] is also shown. The top terminal amino acids of the MCP proteins that are meant to be in the membrane were immobilized in our simulations.

The sequence of events in the signaling involves the binding of an attractant ligand (nutrient) to the receptor molecules, reducing the activity of the kinase, which decreases the amount of phosphorylated proteins, CheY-P and CheB-P.[10] CheY-P then moves in the cytoplasm and binds to a component of the flagellum known as the switch. This binding increases the probability that the flagellum will spin in the clockwise direction, which makes the bacteria tumble in random directions. When bacteria swim in the direction of increasing attractant concentration, more attractant is bound to the receptor, which decreases the number of CheY-P and results in less CheY-P binding with the switch, which makes the bacteria move in the same direction.[11]

One of the central and outstanding questions in bacterial chemotaxis concerns how the receptors can modulate the activity of the CheA kinase that is more than 200 Å away. This is extremely important from the perspectives of understanding the molecular details of chemotaxis activation,[13] engineering chemotaxis,[14] and also developing an understanding of allostery and signaling in proteins. The identification and annotation of the structural intermediates are important for developing functional insights. However, solving the large structural complex or studying it computationally has not been easy. The advances in the observations of the structural complexes[15] and the detailed molecular architecture[9,16] have been fairly recent and mainly through cryoelectron microscopy (cryoEM) studies. Recent studies combining cryoEM with molecular simulations[12] identified the presence of a dipped state of the P4 domain of the CheA protein and further validated it using generalized simulated annealing calculations.[17] The dipped state is expected to play a significant role in the chemotaxis mechanism, although the details are not clear.

In principle, computational studies can shed light on the mechanistic details. However, the signal generation by methylation or its propagation to the kinase domain or the subsequent modulation of the ATP hydrolysis has not been studied computationally. The very large structural complex inspires coarse-grained simulations;[18] however, such simulations are not suited for atomistic investigations that involve methylation or ATP. In fact, even the ATP-bound MCP–kinase domain structural complexes have not yet been simulated.[12] Signaling and allostery in globular proteins have received some attention in computational studies.[19,20] However, signaling in fibrillar proteins is a relatively unexplored territory, which raises questions even about the suitability of the force fields to study the signaling mechanisms in fibrillar proteins. Historically several such questions arose on the suitability of the force fields,[21−23] when longer simulations became possible,[24] or when new types of structures such as those with intrinsic disorder were investigated[23] and even when considering ATP hydrolysis.[25] Thus, there are several stages of work that are required before interesting mechanisms can be unveiled by simulations.

In this work, we study the MCP–kinase complex of the chemotaxis in search of evidence in the simulations for biochemical observations about the signaling. We perform and analyze all-atom molecular dynamics on the different systems with the specific conditions simulating methylation/no-methylation and ATP-bound/no ATP. Given the difficulty of capturing and quantifying the effects of methylation over distances of 180 Å and the associated long time scales, the scope of the present work is to take the first step in establishing early evidence for the relation between the methylation and potential kinase activity in the MCP–CheA–CheW complex of Escherichia coli.

Results and Discussion

Defining the Model System for Studying Kinase Activity

The kinase control region in which we are interested consists of a group methyl accepting chemotaxis proteins (MCPs), coupling protein CheW, and histidine kinase CheA. CheA protein is a multidomain protein consisting of P1 phosphoryl transfer domain, P2 substrate binding domain, P3 dimerization domain, P4 kinase domain, and P5 regulatory domain. The P4 domain has a binding pocket for ATP which transfers the phosphate group to the P2 domain, which in turn transfers it to the CheY protein.[26,27] CheR methylates the MCPs when the nutrients are bound to the receptor, thereby contributing to the increase of the kinase activity.[28] Similarly, demethylation of the MCPs is performed by the CheB proteins.[29,30]

As a step toward developing detailed insights using computational models, we explored the sensitivity of the atomistic simulations to capture the relevant pathways between the methylation and the downstream hydrolysis of ATP. The simulations were aimed at modeling the changes internal to the MCP complex, and transport proteins such as CheR, CheB, and CheY were not included in the model. The models we simulated were based on the simplest system that contained all the relevant details at the atomic level. The model consisted of 20 protein chains, including MCPs, CheW, and the P3, P4, and P5 domains of CheA. The cryoEM structure of the complex from ref (12) was used as a starting structure, and four different scenarios were created—with and without methylation of MCPs and with and without ATP in the P4 domain of CheA. Each of these four model systems contained about 775 000 atoms each, and the simulations on all these system choices were repeated with GROMOS and CHARMM force fields (Materials and Methods Section). Overall, 4.8 μs of all-atom explicit solvent simulations was performed on this 775 000 atom system (4 systems ×150 ns per system per copy × (5 copies with GROMOS + 3 copies with CHARMM)). From each of the copies, the first 50 ns of equilibration was not used in the subsequent analysis, making it an analysis on 3.2 μs of data. Using these simulations, several parameters, such as the P4−P3−P4 interdomain angle (Figure ), the P4−MCP distance, the microenvironment of the hydrogen bonds, and the overall coordination of the ATPγ-phosphate, were computed and analyzed.

Figure 2

Interdomain angle. To understand the flexibility of the domains, the interdomain angle was computed using the centers of mass of the P3 domain and the two P4 domains.

Validating the Simulations with Known Interdomain Contacts

The knowledge of a few static structures may not be sufficient to interpret protein dynamics, functional or nonfunctional. Since the MCP–CheA–CheW protein complex is large, one way of interpreting the dynamic contacts is through mutational studies. A summary of these experimentally identified functional contacts was compiled in Supplementary Table 2 of Cassidy et al.[17] From the different simulations we performed, we compared the contacts at the beginning or end of the simulation. The comparisons are shown in Supporting Information Tables 1 and 2. As can be seen, most of the contacts that were found to be important were explored during the course of our simulations, with a typical cutoff of 8.5 Å for the Cα–Cα distance. Cassidy et al.[12] also show a comparison of these experimentally derived contacts with the contacts they found using simulated annealing. The performance of our calculations was comparable to that from the simulations of Cassidy et al.[12] We consider the occurrence of the several contacts from the mutational studies and from the simulated annealing calculations noted in the literature as a validation of our simulation, and we further analyze the trajectories. The results from the GROMOS simulations are shown in the main article, and the CHARMM force field is shown in the Supporting Information; both are compared in Section At this stage, however, neither the experiments indicate the thermodynamic propensities of the contact formation nor do the simulations in our work or those of Cassidy et al.[12] estimate the thermodynamic free energies of contact formation in the very large systems we chose for this work.

Intermediate Structure Is Stabilized by the ATP

CryoEM studies on the core signaling proteins of bacterial chemotaxis reported two conformations of the P4 domain, having the two P4 units in the CheA protein further apart (undipped) or bringing them closer (dipped).[12] These dipped and undipped forms of the CheA protein are expected to play an important role in the kinase activity[17] by bringing the P4 domain closer to the P1 domain to which phosphoryl group has to be transferred.

We first analyzed the angle between the two P4 domains from our simulations without ATP, similar to the conditions used in the cryoEM as well as the earlier MD simulations.[17] As seen in the histogram of the angle between these two domains (Figure ), the methylated and unmethylated structures both explored the dipped and undipped configurations with comparable frequencies. Incidentally, these profiles resemble the profile observed in Supplementary Figure 5 of ref (17).

Figure 3

Domain dipping effect. Distribution of the angle between P4 domains of chain C and chain E in the CheA protein. For each of these systems with/without methylation and with/without ATP, five copies were simulated. The calculations are with the GROMOS54a8 force field; a similar analysis with CHARMM36 is shown in Supporting Information Figure 1.

We then repeated the analysis with the ATP-bound structures. The results (Figure ) show a clear difference in the flexibility of the P4 domain when ATP is bound. Interestingly, the methylated and ATP-bound structure appears to stabilize an intermediate between the dipped and the undipped configurations. Although the presence of an intermediate structure in the spectrum between the dipped and undipped structures was noted,[17] because the ATP-bound structures have not been studied, the possible functional role of this intermediate structure has not been highlighted. The mechanistic details of how the ATP in the methylated proteins affects the structure or gets hydrolyzed are still unknown. To make an attempt at an understanding of the possible functional role of the intermediate, below we analyzed these structures further.

ATP Binding Increases the Distance between the P4 Domain and MCPs

In all our studies with ATP, ATP was docked in only one of the two CheA-P4 domains (chain C), while the other (chain E) was left in the apo state, giving us an opportunity to compare the role of ATP within the same system. Figure a and b shows the distribution of the distances between the MCPs and the P4 domain (chain C as well as chain E) with and without methylation, respectively. We see that the distance between the P4 domains of chain C and chain O of MCPs is greater with methylation, and the distance between the P4 domains of chain E of CheA and chain I of MCPs is less with methylation and without methylation. These observations suggest that the ATP in the P4 domain of chain C of CheA contributes to an increased distance between the MCPs and the P4 domains.

Figure 4

MCP–P4 distance. Distance between the nearest MCP chain and the P4 domains of the CheA protein. Chain C has a bound ATP, while chain E does not. The calculations are for five copies of the simulation (150 ns each) with the GROMOS54a8 force field; a similar analysis with CHARMM36 is shown in Supporting Information Figure 2.

Methylated Protein Creates Favorable Conditions for ATP Hydrolysis

To understand the microenvironment of the γ-phosphate, we analyzed the number of protein atoms that are within 2.5 Å from the γ-phosphate of ATP. As shown in Figures and 6, we found that with methylation the γ-phosphate has a higher coordination by the protein, leading to a higher possibility of ATP hydrolysis. The three important events in the ATP hydrolysis are the hydrolysis, the separation from the Mg2+ ion, and the downstream signaling involving a transfer of the phosphate group.[26,27,31,32] The sequence or concerted nature of these different events is not yet clear. According to the mechanisms of ATP hydrolysis,[31,32] the affinity of the γ-phosphate of ATP to Mg2+ decreases as its interactions with the amino acids of the P2 domain increase.[32−34]

Figure 5

Local structure near the γ-phosphate. Different representations from a methylated structure where the γ-phosphate is (a) poorly and (b) better coordinated by the protein. Within this simulation, no significant changes in the coordination of the Mg2+ were observed. However, it is likely that a better coordination of the γ-phosphate by the protein is related to its separation from the Mg2+, as is required for the downstream signal transfer.

Figure 6

Microenvironment of the γ-phosphate. To note the differences in the microenvironments, (a) the number of protein atoms within 2.5 Å of the γ-phosphate of ATP were measured. The higher coordination in the methylated case suggests the possibility of easier ATP hydrolysis. (b) Details of the specific contacts in the methylated and unmethylated scenarios are also shown. The calculations are with the GROMOS54a8 force field; a similar analysis with CHARMM36 in shown in Supporting Information Figure 3.

Intermediate Structure Coordinates the γ-Phosphate Better

While the interdomain angles were studied, the intermediate structure was noted specifically in the methylated systems. To develop an insight into its potential, we analyzed this intermediate structure for how the ATP is bound. As shown in Figure , the intermediate structures with a P4–P3–P4 angle in the range of around 80–100° have a higher coordination of the γ-phosphate of the ATP. This higher coordination suggests a higher propensity for the hydrolysis, hinting at a potential functional role for the intermediate structure.

Figure 7

P4 domain angle vs coordination number of the γ-phosphate. The coordination number was calculated as the number of atoms that are within 2.5 Å from the γ-phosphate. A change of this threshold to 3.0 Å did not change the results. A small random scatter has been added for clarity of understanding the density of points with a given coordination number.

Positional Correlation between Functional Amino Acids (with ATP)

The observation that the ATP is in a favorable environment for hydrolysis when the protein is methylated was an encouraging result for this complex system. Following up on this, we used the method of checking for dynamical cross-correlations from the MD simulations to infer if a pathway for the transmission of this methylation signal may be identified. Because the protein complex is large with 20 chains that add up to 5010 amino acids, we identified 274 functionally interesting amino acids and reorganized the correlation map to highlight correlations among these functional groups. These functionally critical amino acids are re-enumerated for convenience in Figure : methylation sites on each of the two trimers of the dimers;[2] glycine hinges in MCPs;[35,36] MCPs lower end hair pin loops;[37] connecting loops between the P3, P4, and P5 domains of CheA;[38] the ATP lid in the P4 domain of CheA;[39] and the amino acids near the ATP,[31,40] all of which have been identified to have a critical role in the kinase-on activity.

Figure 8

Correlation between functional groups. The index in these figures corresponds to the following functional groups (overall 274 Cα atoms) we identified: linker P4–P3 (1–10), loop near ATP (ATP lid) (11–21), protein near ATP (12–35), linker P4–P5 (35–45), lower end loops in MCPs (46–120), Gly residues in MCPs (121– 274), and methyl sites in MCPs (226–274). The inverse mapping of these re-enumerations to the amino acid numbers in the PDB is provided at https://github.com/Himanshu535/MD_Simulaion_data-. The correlation graphs were calculated with five copies of simulations (150 ns each) with the GROMOS54a8 force field; a similar analysis with CHARMM36 in shown in Supporting Information Figure 4.

We probed the correlations among the linkers, the methylation sites, and the ATP binding region, with the ultimate goal of probing the signal transmission pathways. In Figure , one sees strong positive or negative intraregion correlations for all selected regions. Specifically, looking at the proximity of the ATP binding pocket (11−35), it has negative or no correlation with P4−P3 and P4−P5 linkers when there is ATP in the P4 domain, while it has positive correlations with both of these linkers when the P4 domain has no ATP. Interestingly, in the case of protein with methylation and ATP, Figure (d), the correlations of P4 domains with the upper part of MCPs (121−274) are weakened. This is possibly because the dipping and P4 dimer interactions in the protein with methylation and ATP are related to the weak P4−MCP interactions. However, beyond these correlations, no causal relations were explored.

It was also experimentally observed[41] that with methylation the methyl accepting proteins of the MCPs (which are called trimer of dimers 1 (TOD1), trimer of dimers 2 (TOD2)) go apart from each other, and without methylation they come toward each other. Although we observe negative correlations between the TOD1 and TOD2 regions in the presence of the ATP, it is harder to compare with these experiments because in our simulations the top region of the MCPs is immobilized.

Contrasting the Two Force Fields

We performed our calculations with four different scenarios, with each of the two force fields. The results presented in the article are from the GROMOS simulations, while those from the CHARMM force field are shown in the Supporting Information. While the influence of the ATP on the P4–MCP distance was comparable to that with GROMOS, in most other cases such as the angle between the P4 domains (Supporting Information Figure 1), there was no noticeable difference in the methylated and unmethylated cases. This is an important place to note that the calculations we present are from unbiased molecular dynamics of very complex phenomena. Thus, unless extended studies are performed using accelerated sampling methods, it is not possible to comment on the thermodynamic propensities of either what appears to be favorable in GROMOS simulations or what appears to be indistinguishable in CHARMM simulations. However, the early evidence of the simulation results from the GROMOS calculations showing several interesting observations is encouraging and supports the pursuing of detailed thermodynamic studies on this extremely important signal transduction mechanism.

Materials and Methods

Structures for Simulation

We chose to work with the core kinase control region involved in the receptor signaling in E. coli, the MCPs, histidine kinase CheA, and adaptor protein CheW from RCSB (PDB: 3JA6).[12] To bind ATP to the P4 domain of CheA, we performed homology modeling using MODELER-9.21[42] using the PDB structure of the CheA domain from Thermotoga maritima (PDB: 1I5D) as a template. The complex structures used in our simulations are available at https://github.com/Himanshu535/MD_Simulaion_data-.

Earlier studies on the chemotaxis of T. maritima identified Q274, Q498 and E281, E505 as the key methylation targets that influence the chemotaxis.[43] In our simulations using the structure from E. coli, we methylated the equivalent amino acids, all 48 glutamic acids and glutamine in the MCPs (eight methylated sites in each dimer). The protein structures were solvated with simple point charge (SPC) water molecules.[44] In each of the systems used in our simulations, the MCP protein complex had around 5010 amino acids (around 48 000 atoms) and was approximately 16.8 × 11.4 × 27.2 nm in size. After solvation with around 240 000 water molecules, the system size was approximately 20 × 14 × 30 nm. The exact details for each system are given in Supporting Information Table 3.

For each of the four different systems that was designed, two different choices of the force fields GROMOS54a8[45] and CHARMM36[46] were used. For the methylation of glutamic acids and glutamine in MCPs, we used Vienna-PTM 2.0,[47] which is a web server for exploring protein post-translational modifications (PTMs) and provides the output consistent with our choice of GROMOS54a8 force fields. While using the CHARMM36 simulations, because the force fields for methylated amino acids were not easily found, we followed the protocol of replacing Q with E, making the QEQE amino acid occurrence in the wild-type to 4Q to mimic the effects of methylation.[17]

Molecular Dynamics

We performed unbiased molecular dynamics (MD) simulations using the GROMACS package (version 5.1.4).[48] After the systems were solvated, ions were added to neutralize the systems. Energy minimization (EM) using the steepest descent minimization algorithm[49] was first performed on these systems. A velocity-rescaling thermostat[50] with a coupling time of 0.1 ps was applied to maintain a constant temperature of 300 K. A Parrinello–Rahman barostat[51] with a coupling time of 2.0 ps and a reference pressure of 1 bar was used in all simulations. Periodic boundary conditions were applied to the system during the simulation, and the electrostatics was handled using particle mesh Ewald summation. The bonds were constrained using LINKS. The terminal amino acids of the MCPs which were supposed to be in the membrane were immobilized with a harmonic restraint of 1000 kJ/mol/nm2. Other than this, no other restraints were applied on the domains.

Four different systems were used for MD studies—with and without methylation and with and without ATP in one of the P4 domains in the CheA protein. With GROMOS54a8 simulations, each of the systems was simulated as five independent copies of 150 ns each. The first 50 ns from each copy was considered to be the equilibration time and was discarded from our analyses. The trajectories from the five copies were combined and used the 500 ns of each of the four different systems in our analyses. With CHARMM36 simulations, a similar procedure was followed with three copies for all four systems with 150 ns simulation and 50 ns equilibration. However, only eight copies of each system were simulated owing to our limited computational resources and the early observations that the GROMOS force field was having a better correlation with the experimental observations.

Correlations

The positional cross-correlation C between the ith and jth amino acid residues of the protein is calculated from the MD trajectory usingwhere r and r are the position coordinates of the ith and jth α-carbon (Cα) of the amino acids and ⟨ r⟩ and ⟨r⟩ are the corresponding mean positions. The Gromacs suite implementation (gmx-covar) was used to perform these calculations. The correlation coefficient (C) will have values in the range of −1 (perfectly anticorrelated) to 1 (perfectly correlated).

The protein complex simulated has thousands of amino acids, and hence visualizing and interpreting the pairwise correlations are difficult. The correlation maps that are shown were generated by regrouping the different groups of amino acids known for their influence on function. The table mapping the 274 amino acids shown in the figures to the amino acid identification in the PDB 3JA6 is given in the following: https://github.com/Himanshu535/MD_Simulaion_data-/.

Conclusions

Although the sequence of steps that happen in the signal transduction from nutrient binding to the change of flagellar rotation direction are known, the protein complex is too large to be studied by structural or computational methods. Despite very extensive molecular dynamic studies on several proteins, signal transduction, especially in fibrillar proteins, is very poorly studied. In this work, we present early evidence for several observations that suggest that it may be possible to capture the effects of methylation on the ATP hydrolysis in the kinase domain. The present work encourages further detailed studies, especially with the GROMOS force field.