Mechanism of activation and autophosphorylation of a histidine kinase

Histidine kinases (HK) are one of the main prokaryotic signaling systems. Two structurally conserved catalytic domains inside the HK enable autokinase, phosphotransfer, and phosphatase activities. Here, we focus on a detailed mechanistic understanding of the functional cycle of the WalK HK by a multi-scale simulation approach, consisting of classical as well as hybrid QM/MM molecular dynamics simulation. Strikingly, a conformational transition induced solely in DHp leads to the correct activated conformation in CA crucial for autophosphorylation. This finding explains how variable sensor domains induce the transition from inactive to active state. The subsequent autophosphorylation inside DHp proceeds via a penta-coordinated transition state to a protonated phosphohistidine intermediate. This intermediate is consequently deprotonated by a suitable nearby base. The reaction energetics are controlled by the final proton acceptor and presence of a magnesium cation. The slow rates of the process result from the high energy barrier of the conformational transition between inactive and active states. The phosphorylation step exhibits a lower barrier and down-the-hill energetics. Thus, our work suggests a detailed mechanistic model for HK autophosphorylation.

Two component systems (TCS) are one of the major signal transduction mechanisms in bacteria.They regulate the response to a variety of environmental or cellular signals (temperature changes, changes in pH, ligands, etc.) 1 .The individual components are the sensor histidine kinase (HK) that detects the signal and the response regulator (RR) protein that coordinates the response, most commonly by acting as a transcription factor (see Fig. 1).The two proteins communicate via histidine to aspartate phosphoryl-group transfer.Based on domain architectures, evolutionary origin and activities there are numerous variations of TCS 2,3 .While TCS are employed by some eukaryotes, they are notably absent from the animal kingdom.That, paired with their importance to bacteria makes these enzymes promising targets for developing novel compounds that selectively inhibit the growth of bacteria or suppress virulence.For instance, waldiomycin, an angucycline antibiotic, inhibits the HK activity of WalK 4,5 in Staphylococcus aureus, a human pathogen responsible for a variety of acute and chronic diseases [6][7][8] .The molecular signal of this system is still unknown but emanates from the bacterial cell wall 9 .In general, the WalRK system has garnered significant experimental attention since it is conserved across Gram-positive bacteria of the order Firmicutes where it has been shown to be essential for viability in a variety of different species of bacteria.Due to their prevalence and the associated wealth of genomic data, TCS are also a common target of bioinformatics studies to, e.g., investigate TCS complex formation 10 , predict 11 or investigate 12 conformational transitions, or redesign protein signaling 13 .Other extensive studies and reviews highlight the range of TCSs and their activities [14][15][16][17][18] .
In a prototypical TCS (cf.Fig. 1), a stimulus (environmental or cellular signal) is detected by the periplasmic and/or cytoplasmic sensor domains at the N-terminus of the dimeric HK.This signal is propagated from N-to C-terminus along the multidomain HK via a series of conformational transitions ruled by unstable transient interactions 19 .While the structural properties of HKs differ, they all have at the C-terminus a conserved kinase core (~450 amino acids) consisting of the homodimeric dimerization histidine phosphotransfer (DHp) domain and the ATP-binding catalytic domain (CA).Following signal detection, the conserved core adopts an asymmetric conformation such that one of the two subunits of the homodimer is kinase active while the other is inactive (cf.Fig. 2).In the kinase inactive conformation, ATP can enter the CA domain and the binding site of the DHp domain is accessible to a RR for phosphoryl-group transfer.In the kinase active conformation, the RR cannot bind the DHp domain.Here, the gamma-phosphoryl group of the bound ATP of one CA is positioned in close proximity to a specific phosphorylatable histidine of DHp.Two different auto-phosphorylation mechanisms have been attributed to individual HKs: In cis-phosphorylation the ATP from the CA domain phosphorylates its own DHp domain within the homodimer, while in trans-phosphorylation the DHp domain on the other monomer within the homodimer is phosphorylated 20 .It appears that the difference in phosphorylation mechanism is merely structural, based on a left handed versus right handed orientation of the dimeric four-helix bundle that forms the DHp domain.As soon as the histidine is phosphorylated, transfer of this phosphoryl group to an aspartate of a bound RR for communication between the two proteins is possible.Bifunctional HK (e.g., EnvZ) also function as phosphatase for the RR and therefore catalyse the hydrolysis of the phosphoryl group 21 .The activation and inactivation mechanisms of the protein are reviewed in detail in ref. 22.
In this study, the activation mechanism and subsequent autophosphorylation reaction of the WalK HK [also known as YycG 15 ] are explored.WalK is an essential HK involved in the regulation of cell growth and division 23 in low GC-content Gram-positive bacteria such as Bacillus subtilis 24 and S. aureus 25 .Building on the tremendous insights into HK function already gained from structural investigations, we focus on several key-aspects of HK autophosphorylation that remain only partially or poorly understood.To complete our understanding of the autophosphorylation mechanisms we utilize a complex multi-scale computational study on supercomputers to complement previous experimental findings.
First, we aim to understand the conformational dynamics and free energy barrier of the transition from inactive to active conformation.While the initial and final structures of the transition have been well resolved in X-ray experiments, insight into conformational dynamicsexperimentally poorly accessiblecan be provided by molecular dynamics (MD) simulations.Pending on the task, we employ both coarse-grained and regular simulations.Another critical gap in our understanding of the signal transduction event is how signal detection by the N-terminal signal domains induce the conformational transition by the C-terminal catalytic core that leads to activation.We address this gap of knowledge by mimicking effective forces from an incoming signal on the DHp domain and show that these signaling forces trigger the conformational transition of the entire HK.Finally, we complete the picture of HK activation by investigating the mechanism and the energetics of the histidine phosphorylation reaction.This histidine phosphorylation requires an acidic residue near the target histidine [26][27][28][29] .We resolve the sequence of elementary events of the phosphorylation and the reaction energetics by extended-sampling QM/MM simulations on a microsecond scale 30 to distinguish the order of phosphoryl transfer and proton transfer and the influence of a critical magnesium cation 28,31,32 .Taken together, the combination of results for the conformational transition and the chemical phosphorylation step renders the course of energy of the beginning of the TCS phosphorylation cascade, thereby revealing the rate determining step and providing insights on the free energy that is released during the phosphoryl transfer.

Methods
A detailed description of the simulation methods is provided supplementary information.

Molecular dynamics (MD) simulations
The inactive and active states of histidine kinase WalK, an essential HK from Lactobacillus plantarum, were obtained from the protein data bank (PDB ID: 4U7N and 4U7O respectively) 33 .The non-terminal missing residues were modelled using MODELLER version 10.0 34 .The phosphoryl acceptor His391 residue was prepared in the HID state, so that the deprotonated Nϵ atom is a suitable phosphorylation target, and there is a strong hydrogen Fig. 1 | A typical two-component system (TCS) featuring domains for signal recognition, transmission, and catalysis.A stimulus is first detected at the periplasmic sensor domain (red) and this signal is transmitted along the transmembrane helices (blue) and the linker domains (green) before reaching the catalytic core at the C-terminus.The catalytic core comprises the dimerization and histidine phosphotransfer (DHp, purple) and catalytic ATP-binding (CA, yellow) domains.Signal detection results in a phosphoryl transfer reaction from ATP in the CA domain to a conserved histidine in the DHp domain.This phosphoryl group is then transferred to an aspartic acid in the response regulator (orange) protein, resulting in an appropriate cellular response.A labeled cartoon model can be found in Fig. 4. bridge between the protonated Nδ-H and the Glu392 sidechain, which also supports a proton transfer between pHis391 and Glu392.The terminal Glutamate is negatively charged at physiological pH.All MD simulations of WalK HK were conducted at 300 K using GROMACS 2020.4 [35][36][37] patched with PLUMED v2.7 38 , with AMBER99SB-ILDN force-field 39 , and with a time step of 2 fs.Each HK state was solvated with the TIP3P water model 40 .Electrostatic interactions were calculated using the particle-mesh Ewald (PME) method 41 with a direct cutoff of 1.0 nm, and a cutoff of 1.0 nm was used for the van der Waals interactions.All bonds involving hydrogen were constrained by the LINCS algorithm 42 .Each HK-water system was first minimized using the steepest descent algorithm for 1000 steps.Then, NVT equilibration was performed for 1 ns using the velocity-rescaling thermostat 43 , allowing the system to reach temperature convergence.Once the temperature was stabilized, NPT equilibration was carried out for 1 ns with the Parrinello-Rahman barostat 44 , effectively achieving pressure convergence.The 2D metadynamics simulations used two different collective variables (CVs): the difference of root-mean-squared deviations from the inactive and active structures (PDB ID: 4U7N and 4U7O, respectively) termed "ΔRMSD", and a torsional angle defined by the centers of mass (COM) of four different parts of the protein: COM1: chain A residues 454-610, COM2: chain A residues 389-453, COM3: chain B residues 454-610, COM4: chain B residues 389-453.The torsional angle COM1-COM2-COM3-COM4 is then termed "COMTOR".

Inactive to active state conformational transition
To study the large-scale conformational changes involved in the activation pathway of the WalK HK, we used steered MD simulation, and then metadynamics simulations to obtain the free energy profile.We chose these methods since the endpoint conformations of the WalK autophosphorylation reaction have been experimentally determined.In the steered MD we used the difference in RMSD from the two reference states as the collective variable to obtain the intermediate configurations along the transition pathway.The RMSD bias was applied to Cα atoms of the DHp and the active CA domain as detailed in Materials and Methods.Finally, a converged free energy landscape of the conformational transition was obtained by performing a multiple-walker metadynamics simulation with two collective variables, ΔRMSD and additionally, an inter-domain torsional angle COMTOR, to facilitate convergence.
Hybrid QM/MM metadynamics simulations All QM/MM simulations were performed using a combination of locally modified versions [45][46][47] of Gromacs 35,36 and DFTB+ 48 interfaced with Plumed 49 .The initial structural model is based on the crystal structure PDB ID 4U7O 33 ; the details of the initial modelling can be referred to in the SI.The MM region was treated with the AMBER99SB-ILDN force field for the protein 39 and TIP3P water model 40 .The QM region, defined in the Results section, was described with the DFTB3 method 50 using the 3OB parameter set 51 .Special DFTB repulsive potentials were used for the interactions P-N 52 and P-O 53 , specifically calibrated to properly describe the phosphoryl transfer and hydrolysis reactions.The electrostatic interactions between the QM and MM regions were treated by means of electronic embedding, and all of the QM-MM interactions were evaluated with PME, thus no long-range cut-off was applied on the electrostatics.A constant temperature of 300 K was maintained by means of Bussi's thermostat 43 , and the simulations ran at constant volume.PMFs were generated using the multiple walker well-tempered metadynamics protocol 54,55 involving 96 individual QM/MM simulations ("walkers").The convergence of the PMFs was monitored, and the metadynamics sampling was extended to 1 μs in each of the two cases.

Results
The conformational transition is endergonic but creates new favorable inter-domain contacts The steered MD showed that the catalytic ATP-binding domain gradually rotated by 57 °, and the center-of-mass (COM) distance between the DHp and CA domains decreased from approximately 2.55 to 1.92 nm.This brings the phosphorylatable His391 in the DHp domain closer to the ATP-lid (Leu545-Leu577) in the CA domain to form the asymmetrical active state.It should be noted that ATP molecules were not present in these simulations and thus ATP interactions with the HK cannot be quantified and characterized by these simulations.The making and breaking of the interdomain contacts will be reported in more detail below.
The free energy surface resulting from the 2D-MTD simulations is shown in Fig. 3 as a function of ΔRMSD with COMTOR integrated out, while the original 2D free energy surface may be referred to in Fig. S1.Convergence of the Free Energy profile (FEP) for this simulation is shown in Fig. S12.The conformational transition barrier (state C to State D) (rate determining step) remains unchanged.The initial, inactive conformation 'I' is followed by a closely similar conformation 'B' with the same energy.Then, a barrier of 10 kcal/mol has to be overcome at ΔRMSD ≈ 0 nm.That barrier is largely due to the rotation of the CA domain, leading to the intermediate conformation 'C' with an introduction of asymmetry, and lies 6 kcal/mol higher than the initial, inactive state.Note that one of the protein monomers is now beginning the transition into the active conformation, while the other monomer will remain in the inactive conformation throughout our investigation.After that appears the ratelimiting barrier of 16 kcal/mol, which is due to the translation of the CA domain towards the DHp domain, leading to another intermediate conformation 'D'.The final product 'A' lying +10 kcal/mol in free energy higher than the reactant I is reached after passing another, minor barrier.The transition from the state D to A consists mostly of the movement of the loop region Arg558-Thr573 from the interface between the CA and DHp domains, flanking outwards, accompanied by a further approach of the ATP-lid to the phosphorylatable His391 situated in DHp domain.It is in this step that the ATP-lid finally assumes the conformation needed to bind ATP, and consequently, carry out the autophosphorylation reaction.Interestingly, the loop region Arg558-Thr573 is more flexible in the A state than in all of the states B, C and D, as can be inferred from the corresponding Ramachandran plots in Figs.S5-S8.
To characterize the conformational changes during the activation of WalK, we analyzed the formation and breaking of contacts between the different parts of the CA and DHp domains.First, let us turn our attention to any inter-domain salt-brige interactions between the CA and DHp domains.Interestingly, the number of salt bridges increases along the entire process: While the inactive states I and B possess two and one salt bridges, respectively, three salt bridges are observed in each of the intermediate states C and D, before as many as six salt-bridge interactions are identified in the active state A. In addition, there is a new salt bridge between Asp563 and the residue Arg453 on the other WalK monomer in the active state, providing an extra stabilization to the previously mentioned loop region in the active site, thus creating a suitable ATP binding site.All the salt-bridge interactions of the active conformer are shown in Fig. S2.
To describe the interaction between the different parts of CA and DHp domains more generally, we have estimated the formation and breaking of general contacts during the transitions B → C, C → D and D → A. In order to do so, we ran additional unrestrained MD simulations starting from geometries taken from the metadynamics simulation, belonging to each of the states B, C, D and A, and estimated the frequencies of contacts between the individual amino acid residues.These data were used to generate the corresponding differential contact frequency maps, which are shown in Fig. 4.
Four regions in the difference contact frequency maps, where contact formation is favorable, have been identified.The number of CA-DHp contacts between residues Thr395-Gly410 of the DHp and residues Lys475-Pro480 in the CA α-helix increased significantly in the B → C step.So did the contacts between DHp residues Leu422-Asp430 and CA residues Asn460-Lys476, with some reordering observed in the D → A step.For example, the interaction between Val424-Lys476 broke in order for Val424-Arg469, Val424-Met472 and Ser402-Lys476 interactions to form.A further region of interest is between DHp residues Glu392-Ser402 and the CA loop region Phe559-Ile579.Contact pair formation in this region became increasingly more favorable as the two domains were brought closer together in steps C → D → A. Finally, the contact between the α-helical regions 391-399 (αA) and 504-517 (α2) comes into being in the step D → A, and is only found in the active state A.
Finally, any contacts between hydrophobic and aromatic residues were analyzed specifically in the interface region of DHp and CA domains.This was performed by restricting the procedure for creating differential contact maps, to the residues of Val, Phe, Ile, Met, Leu and Tyr.The resulting contact maps shown in Fig. S4 reveal that there are no significant changes of hydrophobic/aromatic contacts between these regions.

Helical bending and catalytic ATP-binding domain rotation occurs concurrently
We next aimed to explain the relationship between the motions of the DHp and CA domains.Above, we had explored a scenario in which HK activation conformation transition occurs via a concerted mechanism as the biasing ΔRMSD potential was applied to both domains.In this approach, the helical bending of the DHp and the rotation of the CA domain necessary for activation occurs simultaneously.In a full-length protein, transition forces stem from N-terminal signal binding domains that are directly tethered to the DHp but not the C-terminal CA domain.We thus investigated an alternative approach where the domains move independently step-wise.The RMSD biasing potential in the steered MD simulation was only applied to the Cα atoms of the DHp domain (cf.Methods).This means that the CA domain was left unperturbed and the structural transition was only induced on DHp by the local perturbation.
Significantly, we observed that the DHp helical bending lead to the rearrangement of the entire protein's 3D structure such that the CA domain rotates concurrently.The CA domain rotates by 50 °, which is just 7 °less than what had occurred when the RMSD bias was applied onto both domains.Additionally, the center of mass distance between the CA and DHp domains decreased from 2.55 nm to 1.95 nm.By the end of the simulation, the HK model had a backbone RMSD of 0.45 nm relative to WalK active state (PDB ID 4U7O) which indicates a large degree of structural similarity.Furthermore, important CA-DHp contacts including Arg568-Glu428, Arg434-Asp502 and Phe559-Asn388 were formed.
Overall, there is a larger number of CA-DHp non-bonded interactions in the final state relative to the starting structure, see Fig. 5.
Thus we observed that the CA domain rotates in response to local force perturbations at the helical bundle of the DHp domain, which essentially mimics the forces exerted by ligand binding to N-terminal sensing domains in the full-length protein.These forces lead to a conformational transition of the DHp domain, which likely acts as a switch to trigger a transition from inactive to active state of the entire HK catalytic core in spite of the earlier observed large transition barrier.

Phosphoryl transfer from ATP to His391
As soon as the active state is established, an autophosphorylation reaction takes place, consisting of a transfer of the γ-phosphate group from ATP to the Nϵ atom of His391 in the DHp domain, followed by a deprotonation of His391.This complex reaction is only possible with His391 in the HID protonation state, which was considered in the current simulations also: His391 can be phosphorylated on the unprotonated Nϵ atom, while the protonated Nδ is stabilized by a strong hydrogen bridge to the Glu392 side chain.Another prerequisite for the autophosphorylation is the presence of an Mg 2+ cation, which in this work was always found coordinating with six oxygen atoms after the QM/MM equilibration (one oxygen each in the side chain of Asn541, in the γ-phosphate, β-phosphate and α-phosphate of ATP as well as in two water molecules).Its mechanism cannot be studied by experimental means in any feasible way, while it poses a serious challenge for a computational investigation, requiring an approach that is both, sufficiently accurate and efficient.
The methodical choice taken in this work is a QM/MM multiplewalker metadynamics simulation employing the semi-empirical densityfunctional approach DFTB as the quantum chemical method.This easily parallelizable protocol makes it possible to reach microsecond sampling, while the accuracy approaches 1 kcal/mol due to a reparametrization of the P-N repulsive potential of DFTB.The QM region consisted of the side chains of His391 and Asn541, the ATP molecule with the coordinated Mg 2+ ion, five nearby water molecules and a suitable proton acceptor, see Fig. 6 (a detailed discussion for the influence of the Mg 2+ ion can be found in the SI).Two different simulations were performed, differing in the identity of the proton acceptor: the side chain of Glu392 in system 1, or an OH -ion as proton acceptor in system 2. The metadynamics simulations involved two collective variables (CV) to describe the progress of the chemical reactions and express the potentials of the mean force (PMF): The O-P-N antisymmetric stretch [i.e., difference of the distances Pγ(ATP)-Oβ(ATP) and Pγ(ATP)-Nϵ(His391)] describes the transfer of the phosphoryl group, while the N-H-O antisymmetric stretch, [i.e., the difference of distances Nδ(His391)-Nδ(His391) and Hδ(His391)-O(proton acceptor)] describes the transfer of the proton to the acceptor, see Fig. 6B.For illustration, a negative O-P-N means that the γ-phosphate group has transferred from ATP to His391, and a positive N-H-O denotes a completed proton transfer to the acceptor.
During the metadynamics simulation, the coordination sphere of Mg 2+ always involved five or six ligands.In the final product, Mg 2+ is coordinated with oxygens of three different water molecules, those of the β-phosphate and α-phosphate groups of ATP and, interestingly, an oxygen of the phosphohistidine.To further investigate the role of Mg 2+ in the reaction, we additionally performed a restrained 2D QM/MM metadynamics simulation as follows: Five ligands of the coordination sphere of the Mg 2+ were kept fixed (in accordance to product structure) and an additional water molecule (sixth ligand) was allowed to freely move in and out of the coordination sphere of Mg 2+ and we monitored the free energy.The obtained PMF shown in Fig. S9 reveals a tiny energy difference between the coordination numbers of five and six, as well as a very low barrier of 2 kcal/mol to the un-and re-binding of the sixth ligand.This could mean that pHis is (almost) free to be released from the coordination sphere, and that would help the processes following the phosphoryl transfer: The protein domains could easily move away from each other if the coordination bond to be broken is not quite as strong.

Mechanism of His Phosphorylation, and Nature of the Transition State
To analyze the appearance of the transition state, another QM/MM metadynamics simulation was performed.This 2D metadynamics simulation used a pair of CVs designed to describe the phosphoryl transfer: the distances Pγ(ATP)-Nϵ(His391) and Pγ(ATP)-Oβ(ATP).
The resulting PMF is shown in Fig. 7A.In the PMF can be observed both the associative (formation of a five-fold hyper-valent phosphate transition state, tight TS) and the dissociative (formation of metaphosphate, loose TS) mechanisms 56,57 .It appears, though, that the associative mechanism corresponds to the minimum energy pathway, with the tight TS yielding a (lower) energy barrier of 8 kcal/mol above the reactant.A very similar barrier height will be found in the simulation of the entire chemical step below, indicating that this is in fact the real mechanism.
The structure of the tight TS, which exhibits a five-fold hypervalent state of the phosphorus atom, as obtained from the simulation is shown in Fig. 7B.Interestingly, the bonding on the phosphorus atom is asymmetric, with the P-O distance of 2.17 Å being markedly shorter than the P-N distance of 2.53 Å as shown in Fig. 7B.This might be caused by the asymmetric composition of the coordination sphere of the Mg 2+ ion, with the oxygens of the negatively charged phosphate on one side and the oxygens of the neutral Asn541 and waters on the other.
It is interesting to note that a certain asymmetry was found for the phosphoryl transfer reaction in the DesK:DesR HK:RR complex from B. subtilis 58 .Just like the autophosphorylation site of WalK, the phosphoryl transfer site of DesK:DesR contains an Mg 2+ ion, and the authors argued that the location of Mg 2+ closer to the RR-Asp being phosphorylated suppresses the backward phosphoryl transfer, making the reaction unidirectional.In the authophosphorylation of WalK, the Mg 2+ ion is clearly located  close to the ATP (later ADP) molecule all the time because it is by far the most attractive interaction partner, but there is a different sort of asymmetry in the TS: the P-O distance is much shorter than the P-N distance.
Trajtenberg et al. further argued that the distance between the phoshophoryl donor and acceptor modulates the thermodynamics of the reaction, via the different capability of the Mg 2+ ion to stabilize either molecule or both 59 : At longer distance, Mg 2+ can only stabilize either the phosphoreactant or the product, leading to (larger) free energy difference between them.At shorter distance, Mg 2+ can interact with and stabilize both, equalizing the reaction free energy somewhat.In the current case of WalK, however, this relationship appears less important because the reaction appears down-the-hill (unlike the phosphoryl transfer from the HK to the RR investigated in Ref. 59) due to the energy gain from the deprotonation of phosphorylated His391.

The chemical step is exergonic and base-dependent
We modelled the reaction considering the conserved Glu392, which is in a close contact with His391, as the proton acceptor.A good convergence of the resulting free energy surface (FES) shown in Fig. 8 (left) is indicated by the analysis presented in Fig. S10; no change in the resulting FES is visible after 700 ns of QM/MM simulation time.Passing over a barrier of 8 kcal/mol, the phosphoryl transfer leads to a protonated phosphohistidine intermediate, lying 20 kcal/mol below the reactant.Then, a nearly barrierless proton transfer from His391 to Glu392 leads to a final product that lies 13 kcal/mol above the protonated intermediate.This indicates that a stronger base is needed as the final proton acceptor to make the reaction sequence exergonic.
To this end, an OH -ion was placed near the Hδ atom of His391, which is the proton to be transferred.A new QM/MM metadynamics simulation was performed, see Fig. 8 (right) for the resulting free energy surface.Also, the analysis in Fig. S11 shows a good convergence of the FES after 700 ns of the QM/MM sampling.The final product lies less than 1 kcal/mol below the intermediate, and is a global minimum of free energy now.The energy barrier to the reaction sequence of 8 kcal/mol is identical to the previous case where considering a glutamate as the proton acceptor.The reason for this is that the higher barrier applies to the phosphoryl group transfer from ATP to histidine, which is exactly the same process in both cases.The subsequent proton transfer passes over a much lower barrier of 4 kcal/mol.

Discussion
HK proteins and their role in TCS signal transduction cascades in bacteria have been widely studied for decades, and structural details of their activation mechanism have emerged more recently 27,28,33,60 .While the structural details of the endpoint conformations, i.e., phosphorylation inactive and phosphorylation active, are known, knowledge on the dynamics of the transition, the coupling between the domains and the relationship between the conformational changes and the subsequent phosphoryl transfer reaction is still limited.In a recent work, Olivieri et al. argued that the activation pathway likely occurs via a "walking" mechanism in which the CA gradually rotates relative to the DHp domain through the simultaneous breaking and formation of inter-domain contacts 12 .The alternative "release and rebind" pathway had a significantly higher energy barrier and hence was unfavorable.A study on CpxA investigated the influence of diphosphorylation of the individual monomers 61 and suggested effects on the equilibrium constant for autophosphorylation vs. formation of ATP.
In the present study, we aimed to build upon past studies and explore the coupled motions of the DHp and CA domains following signal detection in WalK HK.We proposed two distinct activation pathways: (1) helical bending of the DHp domain and CA rotation is concerted or (2) occurs via a step-wise mechanism.Our work supports the hypothesis that the large-scale conformational changes of the subdomains are tightly coupled and therefore the activation of the kinase core is a concerted process.
From many enhanced sampling techniques 62 , we chose metadynamics as the endpoint crystal structures are available and we had an intuition of the transition pathway.
The free energy change during the concerted activation pathway determined from the metadynamics simulation is +10 kcal/mol, and the free energy barrier is 16 kcal/mol.These values likely represent upper bounds as ligand binding might bring about negative contributions as seen in the case of adenylate kinase 63 .The study by Olivieri et al. found quite similar energetics of the conformational transition 12 , in spite of the differences in their and our computational setups: theirs included a bound ATP molecule and used a different method (umbrella sampling) and CVs (the separate RMSDs with respect to the inactive and active conformations), indicating that the resulting energetics may be robust, independent of these computational details.Also, Trajtenberg et al. recently proposed a highenergy phosphotransferase state for another HK, and the energy landscape observed here for the WalK HK shows that such a state may be a common feature present in HK proteins 59 .
In spite of the increase in inter-domain contact pairs, the final active conformation had a larger free energy than the inactive state.Many of the amino acids in the contact pairs are hydrophilic due to the presence of a polar or charged side chain.These interactions are weaker than hydrophobic interactions due to their tendency to also interact with the water molecules surrounding the protein.
Using steered MD, we investigated the alternative step-wise activation mechanism.When the force was applied exclusively to the DHp, the helical bending of DHp induced the rotation of the CA domain by 50 °.Although the rotation angle of the CA domain was slightly less than that from the concerted steered MD, His391 in the DHp was brought closer to the ATPbinding pocket, ready for the phosphoryl-transfer reaction.In addition, essential CA-DHp contacts were formed (Fig. 5).Therefore, the large-scale motions of the DHp and CA domain necessary for activation clearly occurred simultaneously.
In contrast, Marsico et al. demonstrated that while HAMP-DHp activation in the CpxA HK induces CA domain reorientation, in their study only approximately one third of the full transition is achieved 64 .A key difference between CpxA and WalK HK lies in the directionality of the autophosphorylation reaction: WalK HK undergoes cis-phosphorylation whereas CpxA HK undergoes trans-phosphorylation.For HK proteins that phosphorylate in cis, each monomer phosphorylates itself.Conversely, trans-phosphorylation occurs when each monomer in the dimer phosphorylates the other subunit.It is believed that whether phosphorylation occurs in cis or trans is dependent on the handedness of the hairpin loop that connects the individual DHp helices 20 .Moreover, our WalK HK model did not include HAMP as the crystal structure was not readily available.Therefore, the implications of the presence of a linker domain such as HAMP cannot be determined.Even given these caveats, our results hint at a more general signal transduction mechanism of physiological relevance for HK proteins, with specifics of the HK changing its details.
Regarding the chemical step of histidine phosphorylation, we performed two different semi-empirical QM/MM metadynamics simulations of the process combining the phosphoryl transfer from ATP to His391 and the proton transfer from His391 either to the initially deprotonated Glu392 or to an OH − anion placed near His391.The resulting free energy surfaces converged after simulations were extended to 1 μs.We observe a free energy barrier to the chemical reaction of 8 kcal/mol, independent of the identity of the proton acceptor, leading to a stable intermediate represented by protonated phosphorylated His391.
Adding this value to the free energy of the active conformation of HK (+10 kcal/mol), which represents the reactant state of the phosphoryl transfer reaction, leads to the overall free energy barrier of 18 kcal/mol.The experimentally reported catalytic rate of 0.027 min − 129 corresponds, using simple transition state theory expression for rate ¼ kT=h Á exp ÀE A =kT Â Ã , to an activation energy (barrier) of 22 kcal/mol.Such an agreement can be considered favorable, given that any uncertainties of the individual simulations (of the conformational transition and of the chemical step) in fact add up, as well as given that an approximate rate expression is used with an experimental catalytic rate.Since we here aim largely at qualitative interpretations of the observed energy relationships, the discussions of the thermodynamic and kinetic parameters remain valid even if their uncertainties are considered equal to the above "deviation" of our computational barrier from the converted experimental value, which is 4 kcal/mol.
Notably, the barrier opposing the conformational transition obtained from our classical simulations, of 16 kcal/mol, represents the major part of the overall barrier.Therefore, the conformational transition is the ratelimiting step of the whole catalytic activity, as long as that proceeds from the inactive protein structure in equilibrium.Note that this corresponds to the situation in an in vitro experiment, whereas the genuine in vivo process is being triggered by an energy input from the sensory domain reacting to a stimulus, and this also effectively reduces the barrier.
When terminating the reaction with the protonated Glu392, the product lies above the protonated pHis intermediate in free energy.This suggests that a stronger base has to be present as a real final proton acceptor, and this situation was considered in the second QM/MM metadynamics simulation, which had an OH − anion located near His391.With this setup, the deprotonated phosphorylated His391 is a stable product, the species lying the lowest in free energy, 20 kcal/mol below the initial reactant, and even 10 kcal/mol below the inactive conformation of the protein, see Fig. 9. Therefore, for the phosphorylation to occur more readily, there has to be a strong base present as a final proton acceptor somewhere in the system.It is not necessary for it to be directly near the histidine, as long as a proton transfer pathway between the histidine and the proton acceptor is available.
No matter what the pathway is like, the reaction energy will likely be similar; note that proton transfer may occur along rather long "water wires" exhibiting low energy barriers 65 .
We like to emphasize that the two simulations performed with different proton acceptors serve different purposes: The simulation with Glu392 as the acceptor shows the phosphorylation followed by proton transfer to Glu, which is likely the first step of the potentially complex deprotonation of the histidine.The other simulation with OH − representing a general strong base aims to estimate the reaction energy of the entire process, involving a real final proton acceptor.We do not claim that an OH - ion is genuinely present in close proximity to His391; most likely, it is not.The former simulation rules out the deprotonated unphosphorylated His391 (R') as a viable intermediate, as it was suggested earlier to be an "activated" phosphoacceptor 28 , because it lies too high in energy in the realistic pathway proceeding via a protonated Glu392.The latter simulation in turn reveals that the chemical step of the autophosphorylation is in fact a down-the-hill process.
Clausen et al. showed that the rate of autophosphorylation in WalK increases with increasing pH, indicating the need for a strong base present as a proton acceptor 29 .In agreement with that, our simulation performed with and without an OH -present provide a means to quantify the effect of a generic strong base available in the system, on the energetics and kinetics of the reaction.Prior computational studies of WalK and CpxA involved computationally considerably more costly DFT methods, which limited these studies to sub-nanosecond time scales 12,64 .The main finding on WalK was the "tight coupling" of the chemical step with the preceding conformational transition, whereby the protonation of His391 is prevented, which would otherwise hinder phosphorylation.The largely increased sampling in our current study unveils quantitative free energies surfaces, making it possible to explore the reaction mechanism in more detail.What we have learnt about the autophosphorylation, consisting of the phosphoryl-group transfer from ATP to His391 followed (at some point) by the deprotonation of pHis, leaves us with two possible scenarios of the whole process: One possibility is that pHis deprotonates immediately after the phosphoryl transfer has taken place.The first proton acceptor, Glu392 acts like a proton relay before the proton eventually is transfered to a sufficiently strong base (here represented by an OH − anion), which need not be located directly next to pHis.In vivo, that strong base might be, for instance, a suitable titratable molecule present in the solution, or an area of local basic environment in the cytosol.The other conceivable scenario is that the base does not act in this step, and pHis does not deprotonate.Then, the next phosphorylation reaction, which is the phosphoryl transfer from pHis to the conserved Asp residue on the RR WalR, would have to proceed from the protonated pHis.Future work will need to answer whether this process would run spontaneously and if so, whether the energy barrier would be low enough to allow for reasonably favorable kinetics of the process.
Looking at the big picture of the phosphorylation cascade in TCS, the autophosphorylation of His391 will be followed by the transfer of the phosphoryl group from His391 to the conserved aspartate residue on the WalR RR.The final product of the reaction sequence in this study, the WalK protein in an active state containing a deprotonated phosphorylated His391, represents a quite stable minimum on the free energy surface, yet still phosphohistidine is a relatively unstable, high-energetic species.This means that the energy of ATP hydrolysis has not been fully released, and is still available to drive consecutive processes, of which the first is a phosphoryl group transfer to the conserved Asp in RR.
It appears likely that the dynamics and energetics of WalK autophosphorylation found here are valid not only for a single HK, rather, broadly representative of the HK autophosphorylation mechanism.We base this on the fact that direct coupling analysis of vast sequence alignments of HK proteins identified highly correlated residue pairings between DHp and CA domains, that later were identified to be in close proximity in individual structural examples of either inactive or active conformation of the HK protein 11 .Thus, the vast majority in sequence alignments included in this study have to have similar active and inactive conformations, in which these contacts can be realized.By extension, we also anticipate similar dynamics in the transition between active and inactive conformations.

Conclusion
The initial processes in the signalling cascade of the WalK HK have been explored by means of a multi-scale simulation approach.Free and steered classical MD simulations were used to describe the conformational transition leading to the activated state.The interchanging patterns of favorable contacts between amino acid residues were analysed extensively, with emphasis on the contacts between the domains CA and DHp.Importantly, combination of two different simulation protocols made it possible to discover a correlation of the motions of the domains DHp and CA domains: Driving a conformation transition solely in the DHp domain induced a concurrent rotation of the CA domain also, and the process led to the experimentally known activated conformation.This thus appears to be one of a series of conformation coupling occurrences in a HK phosphorylation cascade.
The structure of the activated state of WalK HK served as a basis for an investigation of the autophosphorylation reaction, in which the γ-phosphate group of an ATP molecule bound to the CA domain is transferred to the His391 residue of the DHp domain.The applied semiempirical QM/MM MD multiple-walker metadynamics made it possible to achieve microsecond sampling and draw a reliable picture of the mechanism and energetics of the process.The reaction was shown to proceed via a pentacoordinated transition state to a protonated phosphohistidine intermediate, which is consequently deprotonated in favor of a suitable nearby base.The role of the basicity of the final proton acceptor was also described quantitatively.
The processes under study represent the first two steps of a long phosphorylation sequence in TCS, with the considered final product being still merely a high-energy intermediate in the sequence.The combined process is very slow, as known from a large body of previous experimental research, and in accordance with that, the obtained combined potential of the mean force indicated an energy barrier of 18 kcal/mol as seen in Fig. 9.The conformational transition to the active structure constitutes the ratedetermining step that contributes most of the energy barrier, and is endergonic.The phosphorylation step, on the other hand, exhibits downthe-hill energetics, with the exact shape being dependent on the nature of the final proton acceptor.When considering OH − in that role, the overall process combining the conformational transition and the phoshphoryl transfer reaction exhibits a reaction free energy of −10 kcal/mol.That draws a picture with of an exergonic process accompanied by a high energy barrier, in agreement with and extending the current state of knowledge of the reaction.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 2 |
Fig. 2 | Two different processes in an HK are taking place on different time scales.Purple -DHp domain; yellowcatalytic domain; black dotphosphate group.This work aims to investigate the processes I → II (classical atomistic simulation) as well as the phosphorylation in state II (QM/MM simulation).The chemical reaction consists of two distinct steps: II-A → II-B, phosphoryl transfer itself, and II-B → II-C, deprotonation of the phosphorylated histidine.

Fig. 3 |
Fig.3| Transition between inactive and active structures via intermediates.Top: The free energy profile of the conformational transition obtained from the 2D metadynamics simulation as a function of ΔRMSD, with COMTOR integrated out.Bottom: Representative structures of the conformational states found on the free energy profile, obtained as centroids from cluster analyses of additional unrestrained MD simulations started from geometries found in the metadynamics simulation.Iinactive structure; B, C, D-intermediates; Aactive structures.The loop region Arg558--Thr573 is highlighted in red.Green arrows are pointing at the large-scale movements observed.

Fig. 4 |
Fig. 4 | Concerted activation pathway of histidine kinase.Top: Inter-domain differential contact frequency maps for the interface between the DHp and CA domains.The residues in the DHp domain are on the horizontal axis, and the CA domain is on the vertical, with the ATP-lid highlighted in green.The three different plots are separate maps for each of the three consecutive conformational transitions occurring during the activation process: B → C, C → D and D → A, respectively.Contacts are defined as non-bonded interactions between heavy atoms within a cutoff distance of 0.45 nm.The color-coded values range between-1 (for contacts present entirely in the initial state and missing entirely in the final state) and + 1 (for contacts missing entirely in the initial state and present entirely in the final state).Bottom: Representative structures of the states I and A. The two DHp domains of the homodimer form the central hour-helical bundle, the two catalytic domains are connected by a linker region on the left and right of the helical bundle.Highlighted are the α-helical segments of one monomer: α1 (residues 460-476), α2 (502-520) and α3 (575-588) in blue; αA (383-413) and αB (415-445) in green.Contacts that only occur in the active state A are encircled in the contact maps, and highlighted in orange in the structure of state A.

Fig. 5 |
Fig.5 |Step-wise activation pathway of histidine kinase.Top: The initial and final conformations (violet) from steered MD superimposed with the crystal active WalK histidine kinase (PDB ID 4U7O; grey).Bottom: Contact maps of CA-DHp interactions of the initial (left) and final (right) states.Each point on the contact map represents the presence (green square) or absence (no marking) of contacts between residues of the CA and DHp domains of chain B. Contacts are defined as non-bonded interactions between heavy atoms within a cut-off distance of 0.45 nm.

Fig. 6 |
Fig. 6 | Model of WalK used as the initial structure for QM/MM metadynamics simulations of autophosphorylation.A The active structure adopted from PDB ID 4U7O with the non-reacting ATP-binding domain truncated; location of the reaction center highlighted in pink.B The QM region covering the reaction center; the antisymmetric stretch CVs presented in blue and pink.C Coordination sphere of the magnesium cation in the reactant structure involves three ATP oxygen atoms (one from each of the α-, βand γ-phosphate groups, going from top-right to bottom-left), the oxygen of the side chain of Asn541, and two water oxygens.

Fig. 7 |
Fig. 7 | Results from 2D QM/MM metadynamics simulation using the distances P-N and P-O as CVs.Left: The 2D potentials of the mean force; free energy color coded in kcal/mol, and contour lined separated by 2 kcal/mol.Bottom: A representative transition state structure from that simulation; highlighted are P-N and P-O distances (thick solid lines) and the six coordination bonds to the Mg 2+ ion (thin dashed lines).

Fig. 8 |
Fig. 8 | Results from the 2D QM/MM metadynamics simulations of autophosphorylation of WalK, using the antisymmetric stretches N-H-O and O-P-N as CVs.There are two different simulations: one involving the side chain of Glu392 as the proton acceptor, and the other with an OH − ion playing that role.Top: Potentials of the mean force for the phosphorylation reaction.The free energy is color-coded, and the spacing of contour lines is 3 kcal/mol.Bottom: Representative structures of the states found on the PMF: R-reactant, R'-intermediate (His391 is deprotonated and unphosphorylated), I-intermediate (protonated phosphorylated His391), P-final product (deprotonated phosphohistidine).

Fig. 9 |
Fig. 9 | Overall free energy profile.The free energy profile of the entire autophosphorylation sequence including the conformational transition from the inactive to the active structure, and the phosphorylation of the His391 residue.