Dynamic energy correlation analysis of E. coli aspartokinase III and alteration of allosteric regulation by manipulating energy transduction pathways

Abstract Conformational change associated with allosteric regulation in a protein is ultimately driven by energy transformation. However, little is known about the latter process. In this work, we combined steered molecular dynamics simulations and sequence conservation analysis to investigate the conformational changes and energy transformation in the allosteric enzyme aspartokinase III (AK III) from Escherichia coli. Correlation analysis of energy change at residue level indicated significant transformation between electrostatic energy and dihedral angle energy during the allosteric regulation. Key amino acid residues located in the corresponding energy transduction pathways were identified by dynamic energy correlation analysis. To verify their functions, residues with a high energy correlation in the pathways were altered and their effects on allosteric regulation of AKIII were determined. This study sheds new insights into energy transformation during allosteric regulation of AK III and proposes a strategy to identify key residues that are involved in intramolecular energy transduction and thus in driving the allosteric process.


INTRODUCTION
As a typical dynamical model for investigating the relationship between protein structure and function, allosteric proteins have been studied for decades and are still of great interest [1][2][3]. Although static structures are known for for proteins (CURP) was proposed for the analysis of a neuronal protein PSD-95 [7]. An energy dissipation model has been proposed to investigate Escherichia coli (E. coli) aspartokinase III (AK III) mutations for better understanding the desensitization of product feedback inhibition via allostery [8].
As a key enzyme in amino acid biosynthesis, aspartokinase catalyzes the phosphorylation of the amino acid aspartate. This reaction is the first step in the biosynthesis of three essential amino acids: methionine, lysine, and threonine. In E. coli, there are three independently regulated isozymes of aspartokinase, each of which is specific to one of the three downstream biochemical pathways [9]. This allows independent regulations of synthesis of methionine, lysine, and threonine respectively. Among them, the enzyme AKIII coded by the gene lysC is l-lysinesensitive. It consists of a dimer of protein encoded by lysC.
In the presence of lysine, it can adopt an inactive tetrameric conformation [10].
For better understanding the mechanisms of allosteric regulation, E. coli AKIII is employed as a model system in this study. We aimed to reveal the structural and energetic features of the key residues in the dynamic process of allosteric regulation of this enzyme. To this end, steered molecular dynamics simulations were performed to figure out residues involved in the allosteric regulation and this was followed by energy transduction analysis. In addition, sequence conservation analysis of the aspartokinase family was carried out to calculate the conservation score of each amino acid residue and to identify conserved sites. Key amino acid residues that participate in the allostery process were identified to be those holding a high energy correlation score with both of the lysine binding site and the aspartic acid binding site. To verify the predictions, AK III mutants with mutation of key residues were experimentally studied.

System preparation
The x-ray diffraction structures of E. coli AK III were retrieved from the Protein Data Bank (PDB). In this work, the crystal structure of the T-state (PDB code 2J0X) was employed for study. After deleting the lysine located at the regulatory site, the T-state structure was neutralized by adding sodium and chlorine ions with an ionic concentration of 0.5 mol/L and solvated in a rectangular box of TIP3P water molecules with a minimum solute-wall distance of 20 Å. The solvated system was energy-minimized by 5000 steps employing the software of NAMD.

PRACTICAL APPLICATION
In this work, we combined steered molecular dynamics simulations and sequence conservation analysis to investigate the conformational changes and energy transformation in the allosteric enzyme aspartokinase III (AK III) from Escherichia coli. Energy correlation analysis indicates dynamic transformation and compensation between electrostatic energy and dihedral angle energy during allosteric regulation of AKIII. Key residues involved in the allostery are identified to exhibit a high energy correlation score with both of the lysine and aspartic acid binding sites. Energy transduction between the binding sites is triggered by ligand binding/leaving, leading to amplitude structural motions in the protein. The energy transduction pathways were redesigned to modulate the allosteric regulation and validations were conducted by site-directed mutagenesis. This work highlights the importance of structural dynamics and energy landscape of conformations and provides new insights for kinase engineering and design.

Molecular dynamics simulations
In this work, steered molecular dynamics (SMD) simulations were performed by auto-searching the pathway of energy minimization. With a periodic boundary condition in the NPT ensemble, Langevin dynamics applied at 310 K with the damping coefficient of 5.0 ps-1 and constant pressure of 1 atm. No constraint was applied to the protein during the molecular dynamics simulations. A time step of 2 fs was used and the coordinates of the simulated complexes were saved every 1.0 ps. The simulations lasted for 10 ns and were performed employing the software Amber with the CHARMM27 force field. The Collective Variable-based Calculations (Colvars) module embedded in the NAMD software was used to conduct the SMD. A harmonic restraint was set with the force constant of 10.0 kcal/mol and the ligand was pulled out of the binding pocket until its center from that of the protein was extended by 10 Å. Parameters were as follow: force-Constant 10.0; centers 9.0310; targetCenters 19.0310; tar-getNumSteps 5000000. Lysine was completely pulled out from AKIII at the end of simulation. In total, 2000 snapshots were obtained. More details of dynamic simulation are submitted as Supplementary Material. Protein side-chain networks (PScN) were then constructed by considering amino acid residues as nodes; edges between the nodes were added on the basis of noncovalent interactions and side-chain interaction between them [8]. A residue-based cut-off for the non-covalent interactions was 20 Å.

Dynamic correlation analysis
Dynamic correlation analysis was employed to identify residues that mediate the energy transduction within AKIII. This approach is based on the correlation coefficients between each pair of residues which were obtained from the trajectory of steered molecular dynamics simulations. The correlation coefficients are calculated according to the following equation [8]: where C(Rs,Rt) is the correlation coefficient between the residues s and t. Rs and Rt are the coordinates of the residues s and t; 〈Rs〉 N and 〈Rt〉 N are the average coordinates of residues s and t calculated based on the whole trajectory and N is the total number of snapshots recorded in the trajectory. The correlation coefficients of distance, energy and energy change in small time scale between each pair of residues were calculated using Matlab. Correlation analysis of energy was carried out by using Matlab function corr2 (A, B), which computes the correlation coefficient between arrays A and B. The so-called C Value, which indicates the energy correlation score of residue X with both residue A and B, was calculated using Equation 2:

Conservation analysis
Conservation analysis of AKIII (PDB ID: 2j0x) was carried out by using the ConSurf webserver (http://consurf. tau.ac.il/2016/) with HMMER as the Homolog search algorithm and the HMMER E-value was set to be 0.0001.There are 902 HMMER hits of sequences. Amino acid sequences which are similar to AK III in the PDB were collected and multiply aligned using CSI-BLAST and MAFFT, respectively. The evolutionary conservation of each amino acid position in the alignment was calculated using the Rate4Site algorithm. The algorithm takes explicitly into account the phylogenetic relations between the aligned proteins and the stochastic nature of the evolutionary process. Multiple sequence alignment was built using MAFFT. The homologues were collected from SWISS-PROT. Homolog search algorithm is BLAST. PSI-BLAST E-value was 0.001. Maximal %ID between sequences is 95 Minimal %ID for homologs is 35.

Enzyme activity assay and site-directed mutagenesis
Activity assay and site-directed mutagenesis of AK III were carried out according to Chen et al. [10]. Specifically, AK III encoded by the gene lysC was obtained from E. coli K12 MG1655 genome and inserted into pET22b at the NdeI and XhoI restrict enzyme sites. For over-expression and purification of the wild-type AK III and relevant mutants, the wild-type lysC gene was cloned to pET-22b(+) (Novagen, Darmstadt, Germany) with the introduction of an additional His-tag at the C-terminal to generate the plasmid pET22-lysC. Strain cultivations were performed in shake flasks at 30 • C and 250 rpm. The activity was measured with a method adapted from Blank and Wright (1954). The preparation for the enzyme reaction assay was composed of 200 mM Tris-HCl (pH 7.5), 10 mM MgSO 4 ⋅6H 2 O, 10 mM aspartate, 10 mM ATP, 160 mM NH 2 OH⋅HCl (neutralized with KOH), and appropriate amounts of enzyme. Total reaction volume was 500 µL. After incubation at 37 • C for 30 min, the reaction was stopped by mixing with 1 mL of a 5% (wt/vol) FeCl 3 solution, and the absorbance at 540 nm was monitored.
Site-directed point mutations were generated by PCR reaction using QuickChange site-directed mutagenesis kit (Agilent Technologies, Germany). PCR reactions were carried out with a pair of synthetic complementary oligonucleotides containing the desired mutation as primers and a plasmid with the inserted gene as a template. A mutant plasmid containing staggered nicks was generated by extension of the oligonucleotide primers during temperature cycles. All PCR reactions were performed using Pfu or Phusion polymerase obtained from Fermentas (Germany) with a DNA Thermal Cycler DOPPIO (VWR, Germany).

Steered molecular dynamics simulations of E. coli AK III
Steered molecular dynamics simulations apply forces to a molecule in order to manipulate its structure by pulling F I G U R E 1 (A) Residues with a RMSD value higher than 10 Å referred from the steered molecular dynamics simulations; (B) Residues with total energy change higher than 8 kJ/mol (ID3, 9,  it along desired degrees of freedom. This approach can be used to study the allosteric regulation in a large timescale [7]. In this study steered molecular dynamics simulation was employed to study the process of pulling lysine out of the ligand binding site of AKIII. Backbone root meansquare deviations (RMSD) were calculated based on the trajectory with the starting configuration as the reference and all coordinate frames from the trajectories were first superimposed on the reference conformation. Structural analysis based on RMSD was carried out. Residues with a RMSD value higher than 10 Å include 54, 55, 59, 84, 88, 95, 98, 99, 101, 102, 103, and 104 ( Figure 1A).
Energy variation analysis was also carried out for each residue. Residues with a total energy change of more than 8 kJ/mol are shown in Figure 1B and they are supposed to trigger energy redistribution for driving the allosteric regulation of AKIII. It is noteworthy that most of the residues did not exhibit a high RMSD value, indicating that the structure/position-based parameters such as RMSD are not able to explore the energy response mechanism that undergoes an allosteric process.

Energy transformation during the allosteric process
For a given protein, the energy function that defines its so-called energy landscape is a combination of the conformational energy and non-bond energy [13]. The conformational energy is related to the covalent interactions, namely, the energy of bonds, bond angles and dihedral angle energy. The non-bond energy includes the van der Waals', and electrostatic energy. In the case of AKIII, the correlation coefficient of conformational energy and nonbond energy of all residues was calculated to be -0.45 dur-ing the steered dynamic process. Correlation coefficient which is less than -0.3, indicates that transformation and compensation occur among these two energy items during the allostery procedure.
Subsequently, correlation analysis was carried out to define the major driving energy. Correlation analysis of all the energy items listed in Table 1. From Table 1, it can be seen that the correlation coefficient between nonbond energy and electrostatic energy is as high as 0.999, which indicated that the electrostatic energy is the main source for non-bond energy redistribution [14]. Then, the correlation coefficient of conformational energy and dihedral angle energy was calculated to be 0.627, indicating that dihedral angle energy makes great contributions to the conformational energy change. And the correlation coefficient of electrostatic energy and dihedral angle energy is -0.414, implying that the electrostatic interaction energy and dihedral angle energy are transformed to each other and may drive the allosteric modulation. A recent study also proposed that allostery in the well-studied PDZ domain is driven by changes in electrostatic effects rather than solely changes in dynamics [15].
Due to the compensating effects of different energy items, namely electrostatic interaction energy and dihedral angle energy, the total energy of the protein is perturbed to a much lower extent. Electrostatic interaction energy becomes the key determinant in distinguishing the structural ensembles between the ligand bound and unbound states, and in elucidating the allosteric modulation. The electrostatic interaction redistribution is a common feature among "dynamic-driven" allosteric proteins [16]. However, electrostatic interactions are the main source for non-bound energetic redistribution, while van de Waals interactions and H-bonded networks are both involved in energetic redistribution to some extent [17]. Residues with TA B L E 1 Correlation coefficients among different items of energy obtained from correlation analysis of energy transformation and redistribution during steered molecule dynamic simulation  [18]. These results indicated that the mutation can change the energy transformation and then affect the allostery.

Energy change and conservation of the binding sites
The energy change of residues in the starting and final states of the steered MD calculations was calculated. The changes of lysine, aspartate and ADP binding sites during the steered molecule dynamic simulation are listed in Tables 2-4. For the lysine binding site, there are significant changes of the electrostatic energy and dihedral angle energy of residue D340, which were -4.656 and 5.204 kJ/mol, respectively. Residue 338 experienced a large electrostatic energy change of -5.134 kJ/mol. According to the crystal structure of AKIII, D340 and E346 can form two salt bridges with lysine. Mutation of S338L, E346R and V339A were reported to deregulate allosteric process by affecting the binding of lysine [10]. Obviously, redistribution of internal energies happened within the protein, especially the electrostatic energy.
For the aspartate binding site, there are significant electrostatic energy changes in residues 30, 199 and 201, which were -8.241, -4.943, -4.695 kJ/mol, respectively. S201 binds with both aspartate and ADP. However, there is much less electrostatic energy change in the ADP binding site, except the residue 201, which is also shared with the aspartate binding site. Therefore, the aspartate binding site undergoes a greater energy change than the ADP binding site.
As for the conservation of the residues located at the binding sites, the lysine binding site is less conserved than the aspartate/ADP binding sites according to their average conservation scores. For the lysine binding site, its average conservation score was -0.849, while for the binding sites of aspartate and ADP, they were -1.062, and -0.949, respectively. Several studies have been reported about altering the inhibitor from lysine to norlecuine et al. by mutation of the ligand binding site [10,11]. Engineered proteins were TA B L E 2 Energy change of lysine binding residues (kJ/mol)

Electrostatic energy correlation analysis
The release of lysine from the binding pocket causes energy redistribution, which is associated with the conforma-tional change. As a converse process, this is also true for the binding process of the ligand into the pocket. To figure out how this perturbation is propagated within the protein, correlation coefficients of the electrostatic energy were calculated for each residue. Then, all amino acid residues were clustered using the software Cytoscape by considering amino acid residues as nodes and electrostatic energy correlation scores as edges.

F I G U R E 3
Key residues referred from intersects of correlation coefficient of electrostatic energy less than -0.5 (blue circles) and more than 0.5 (red circles) Residues 341, 340, 292, 75, 245 and 296 were clustered with electrostatic energy correlation less than -0.5, which obviously indicated an electrostatic energy transduction among them. And residues 443, 245, 340, 165, 306, 106, 201, 95 were clustered with an electrostatic energy correlation coefficient which is more than 0.5, which indicated that they may act as an energy pool for driving the conformational change. Of those residues, F407 was reported to be in the protein dynamical modules based on energy dissipation [8]. V347M has previously been reported to be able to deregulate AKIII through the approach of random mutation [17]. As given in Figure 3, residues 106, 201, 306, 245, and 340 have high electrostatic energy correlation scores. Among them, residue 201 is located in the aspartate binding site. These results indicated that correlation analysis of electrostatic energy can provide a highly sensitive probe toward identifying the key residues in the energy redistribution [15].

Energy redistribution pathways
To reveal how the energy is transferred in proteins through the interaction network, 12 pathways were proposed from the key residues of lysine binding site (Residue A: S338, D340, E346) to aspartate binding site (Residue B: S39, G119, D202, S201) were studied based on energy correlation coefficients of each residues. Only residues that obtain high energy correlation coefficients with both of the beginning residue and the end residue will obtain high C scores.
Residues with a C score of higher than 0.3 are considered as key residues in the energy transduction pathway. Residues in the pathway between E346 to S201 showed the highest correlation among all the 12 pathways revealed by SMD. This indicates that it is the main pathway of energy transfer. The pathway from E346 to S39 includes 26 residues with C score more than 0.3. The pathways from D340 to S39 and S201 also showed a high correlation. Due to the fact that S39 and S201 are the binding residues shared by both of the aspartate and the ADP binding sites, it can be assumed that energy redistribution pathways which connect the lysine and aspartate binding sites . And the 345-222 route is also a key signal transfer route according to the previous work done by Ma et al. [8].
Shared amino acid residues in the routes of D340 to S201, E346 to S201 and E346 to S39 includes (ID 106,158,160,165,245,260,292,296,306,341 and 443), which formed two big clusters. Key residues which shared in all of the four pathways (D340, E346 to S201, S39) were residues 106, 165, 245, 292, 296, 306 and 443. The average conserved scores of the residues (ID 106, 165, 245, 292, 296, 306 and 443) which shared in many pathways have a high value of -0.0471. This is moderately conserved, while those residues which only exist in a single pathway exhibited less conserved. For example, 19 residues in the pathway of E346 to S201 with an average value of 0.316 (Figure 4). This indicates that the shared residues in the energy transduction pathway are relatively more conserved.
Besides, the residues 106, 245 and 306 also show a high value of correlation coefficient for the electrostatic energy which is more than 0.5. Residue 306, which belongs to the ACT1 domain, contributes the major energy transduction in the pathway based on the electrostatic energy redistribution. These results are in agreement with the dynamical mechanism for allosteric transition and inhibition by lysine in AKIII proposed by Kotaka et al. [9].

Experimental verification by site-directed mutagenesis
Most of the mutations reported in the literature (Table 5) which are able to deregulate the feedback inhibition of AK III are around the lysine binding site [19]. Few mutations were reported that are within the energy transduction pathways. In this study, key residues in the energy transduction pathways were selected according to the analysis of energy correlation. Several residues that are not located within the active site motifs were selected for mutation F I G U R E 4 (A) Key residues in the energy transfer route of 346-201 with high correlation coefficients of electrostatic energy more than 0.3; (B) Key residues in the energy transfer route of 340-201 with high correlation coefficients of electrostatic energy more than 0.3 ( Figure 5). Thus, the experimental measurements that follow in this work would not arise the problem of changing substrate binding. Allosteric regulation has been affect by point mutations in the energy change pathway.

TA B L E 5 Reported mutations of AKIII
Five key residues have been mutated, namely 106, 165, 245, 292, and 296, which were shared in the all four pathways of D340, E346 to S201, and S39. The relative activities were calculated based on their initial activity without lysine. Mutation of K292D deregulated the inhibition effect, while double mutations of K292D/P296G have almost the same inhibition effect as with the wild type. Residue 245 was highly conserved with a conservation score as -0.8590. The mutation of F245G greatly decreased the enzyme activity, indicating that 245 is a key residue for enzyme activity even it is far away from the binding sites. The mutation of 165 from Leu to the less bulk Gly enhanced the inhibition effect by decreasing Dihedral angle energy. The mutation of T106G deregulated the inhibition effect greatly, which verified the energy transduction pathway we predicted. The mutation of K292D changed the F I G U R E 5 Residues selected for mutation in this work (Residues 106, 165, 292, and 296)  positively charged residue into a negative one, indicated interference of electrostatic interaction deregulated inhibitor effect. These results indicated that electrostatic interaction is a key factor for energy transduction of the allosteric system for AKIII ( Figure 6).

CONCLUDING REMARKS
With AKIII from E. coli as a model system, the dynamic energy transformation between electrostatic energy and dihedral angle energy during allosteric regulation was explored based on steered molecular dynamics simulations. A strategy for identification of key residues involved in energy transduction during the allosteric process was proposed based on the dynamic energy correlation analysis. Mutations were conducted and their effects on allosteric regulation were determined. The results demonstrated that electrostatic energy is the key driving force for allosteric regulation. The intramolecular energy transduction that drives the allosteric process is crucial for a more F I G U R E 6 Relative activities of the wild type and mutants of AKIII with lysine as the effector. The relative activities were calculated based on their activity without lysine holistic and better understanding allosteric mechanism. Energetic features of the key residues in dynamic allosteric regulation process can apply for better understanding of the allosteric regulation of AKIII and for designing and optimizing of allosteric systems with a better performance.

C O N F L I C T O F I N T E R E S T D ATA AVA I L I B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.