Charging of transfer-RNA with cognate amino acid is accomplished by the aminoacyl-tRNA synthetases, and proceeds through an aminoacyl adenylate intermediate. The lysyl-tRNA synthetase has evolved an active site that specifically binds lysine and ATP. Previous molecular dynamics simulations of the heat-inducible Escherichia coli lysyl-tRNA synthetase, LysU, have revealed differences in the binding of ATP and aspects of asymmetry between the nominally equivalent active sites of this dimeric enzyme. The possibility that this asymmetry results in different binding affinities for the ligands is addressed here by a parallel computational and biochemical study.
Biochemical experiments employing isothermal calorimetry, steady-state fluorescence and circular dichroism are used to determine the order and stoichiometries of the lysine and nucleotide binding events, and the associated thermodynamic parameters. An ordered mechanism of substrate addition is found, with lysine having to bind prior to the nucleotide in a magnesium dependent process. Two lysines are found to bind per dimer, and trigger a large conformational change. Subsequent nucleotide binding causes little structural rearrangement and crucially only occurs at a single catalytic site, in accord with the simulations.
Molecular dynamics based free energy calculations of the ATP binding process are used to determine the binding affinities of each site. Significant differences in ATP binding affinities are observed, with only one active site capable of realizing the experimental binding free energy. Half-of-the-sites models in which the nucleotide is only present at one active site achieve their full binding potential irrespective of the subunit choice. This strongly suggests the involvement of an anti-cooperative mechanism. Pathways for relaying information between the two active sites are proposed.
The asymmetry uncovered here appears to be a common feature of oligomeric aminoacyl-tRNA synthetases, and may play an important functional role. We suggest a manner in which catalytic efficiency could be improved by LysU operating in an alternating sites mechanism.
Aminoacyl-tRNA synthetases (aaRSs) catalyse the specific attachment of a given amino acid to its cognate tRNA. As the link between protein structure and nucleic acid information, they are amongst the most ancient of protein families. Direct aminoacylation of tRNA catalysed by a discriminatory aaRS is the most common method of forming aminoacyl-tRNA, although it has recently been revealed that indirect pathways via alternatively aminoacylated tRNA intermediates are other routes to aminoacyl-tRNA . Both of these routes, and all aaRS, share a common first two steps wherein an amino acid is activated by reaction with ATP to form an enzyme-bound aminoacyl adenylate, followed by transfer of the amino acid to its cognate tRNA:
AA + ATP + aaRS aaRS•AA-AMP + PPi
aaRS•AA-AMP aaRS + AA-tRNA + AMP
The reaction is driven by the hydrolysis of ATP. The accuracy of acylation of tRNA is vital to the overall fidelity of protein synthesis, hence an understanding of the mechanism must provide insights into the ways in which aaRSs in particular have evolved to ensure such reliability.
Numerous crystal structures of substrate bound and free states of aminoacyl-tRNA synthetases have been reported [2,3]. These structural studies have enabled aminoacyl-tRNA synthesis to be elucidated at the atomic level. The reaction does not involve covalent or acid-base catalysis, but relies on the binding energy provided by the enzyme to stabilise the transition state . Amino acid specificity is achieved through idiosyncratic interactions of the side chain with the active site, which is shaped to accommodate the amino acid, often through an induced fit mechanism . In cases where discrimination is difficult due to structural similarity, such as between valine and isoleucine, the enzyme has evolved a proofreading mechanism to hydrolyse the misacylated adenylate or tRNA, thus preserving the fidelity of translation .
The twenty canonical aminoacyl-tRNA synthetases have been divided into two classes based on mutually exclusive sequence motifs that result in distinct active site topologies . The class I active site is characterised by a Rossman fold, whilst class II synthetases contain a multi-stranded anti-parallel β-fold, the oligonucleotide binding (OB) fold, surrounded by α-helices . One consequence of the different active site folds is that while ATP is bound in an extended conformation in the Rossman fold, the class II enzymes bind ATP in a bent conformation. Another difference is that the class I enzymes primarily aminoacylate the 2'OH of the tRNA terminal ribose by approaching the minor groove, whereas enzymes belonging to class II approach the major groove and thus charge the 3'OH. Within the class II family, a number of subfamilies have been identified based on more extensive sequence and structural homology [7,8]. Enzymes of both classes are commonly denoted by their three letter amino acid designation, followed by RS.
This study focuses on a lysyl-tRNA synthetase (LysRS) from Escherichia coli. The lysyl-tRNA synthetase, together with the aspartyl- and asparaginyl-tRNA synthetases, constitutes a member of class IIb . Such class membership is generally universal, although an exception has been observed in most archaea and some bacteria, where lysyl-tRNA synthetase contains the characteristic motifs of a class I enzyme . Unlike most prokaryotic aminoacyl tRNA synthetases, the lysyl-tRNA synthetase of E. coli is coded for by two genes, the constitutively expressed lysS and the heat-inducible lysU, which code for proteins sharing 88% sequence identity . Explanations for the evolutionary pressures resulting in two lysyl-tRNA synthetases have included a possible physiological role for LysU in the synthesis of the unusual 5'-5"' linked dinucleotide diadenosine tetraphosphate (Ap4A) under stress conditions, or for LysU being less easily inhibited than LysS by cadaverine, a breakdown product of lysine that accumulates under stress . Most class II aminoacyl-tRNA synthetases are able to synthesise Ap4A, although LysU is most effective, and Ap4A concentration increases dramatically under various stressing conditions .
A series of high resolution crystal structures of LysRS have been published, revealing the enzyme structure at various steps along the aminoacylation reaction pathway [13-16]. LysU is a homodimer, with an extended dimer interface. Each monomer has a molecular weight of 58 kDa and contains 504 amino acids. Monomers are organized into three domains: a smaller N-terminal domain which binds the tRNA anticodon, a larger C-terminal catalytic domain which contains the active site, and an α-helical insertion domain, as illustrated in Figure 1(a). The insertion domain recognises the cognate tRNA acceptor and is unique to the different class II enzymes . The active site is a large cavity on one side of the central β-sheet, and is closed by surrounding loops. The amino acid and ATP ligands, along with three magnesium ions are situated within the cleft. Three sequence motifs characteristic of class II are also found in this region. Motif 1, with a strictly conserved proline forms part of the dimer interface. Motif 2 comprises two antiparallel β-sheets connected by a long flexible loop (the motif 2 loop, shown in green in Figure 1) that only acquires an ordered conformation upon ATP binding. Many residues of this motif interact with ATP. Residues of motif 3 are located in the β-sheet where they interact with the nucleotide phosphates and Mg2+ ions, hence stabilizing the bent conformation of ATP . These structural studies show that LysU optimally positions the substrates for the nucleophilic attack of the lysine carboxylate on the ATP α-phosphate. Reaction occurs through an in-line displacement mechanism to form the lysyl-adenylate with release of pyrophosphate and elimination of a single Mg2+ ion . Onesti's most recent work on the structurally homologous LysS  revealed a major conformational change on binding of lysine, involving the ordering of the two loops and a reorientation of the insertion domain to close up the active site around the substrate. In contrast, ATP binding and the adenylation reaction caused minor structural changes, confined to the motif 2 loop.
Figure 1. Ribbon representations of LysU (a) The 2.4 Å resolution crystal structure of the LysU dimer viewed down the molecular 2-fold axis, with the substrates lysine (pink) and AMP-PCP (cyan) shown in a space filling representation. Mg2+ ions are shown in gold. One complete monomer is shown in blue, whereas the other is coloured according to the domain structure. The N-terminal domain is shown in orange. Within the larger catalytic domain, the conserved core which is common to all class II aminoacyl-tRNA synthetases is shown in red, while the insertion sequence specific to LysRS is shown in purple. In the red monomer, the flipping loop (yellow) is located above the lysine substrate, while the motif 2 loop (green) is situated above the ATP. Strands C2 and C3, along with their symmetry related equivalents form part of a 6-stranded antiparallel β-sheet centred about the two-fold axis. The loop between strands C2-C3 is positioned above and interacts with the motif 2 loop of the opposite monomer. (b) The average structure of the 1 ns FO trajectory. Monomer 2 has been superimposed upon monomer 1 by fitting the core catalytic domain. The same colouring scheme is used. The movement of several loops and domains is apparent, resulting in asymmetry. The motif 2 loop of monomer 1 (green) shifts further over ATP substrate along with the C2-C3 loop of monomer 2. The helix-loop-helix insertion motif also shifts by a few angstroms. Less apparent from this viewpoint is the variation in position of the 4-helix bundle of the insertion, and a rigid body rotation of the N-terminal domain. This figure was created using Molscript  and RASTER3D .
In this paper, we examine these binding events from a thermodynamic perspective, focusing on the ATP binding process. As part of our ongoing investigation into the binding properties of LysU mutants, benchmark molecular dynamics simulations of the wild type enzyme recently revealed an unexpected asymmetry (Hughes et al, manuscript in preparation). In a series of simulations, local nonbonded interactions between ATP and the protein were consistently stronger in one active site. A global asymmetry that was not eliminated by successive improvements to the model was also observed, as shown in Figure 1(b). The intriguing possibility that this asymmetry might be genuine and of functional importance, led us to pursue a parallel biochemical and computational study of the binding events of LysU. Such a study was intended to complement the crystallographic data by adding further mechanistic and thermodynamic detail to the binding events. In this paper we describe how isothermal titration calorimetry, circular dichroism and intrinsic tryptophan fluorescence studies data give compatible and complementary results to molecular dynamics and free energy calculations – whilst lysine binding takes place at both sites in the dimer, nucleotide binding is only thermodynamically favourable at one site, revealing a clear asymmetry in the enzyme. The mechanism by which binding information is communicated from one active site to another is explored.
Molecular dynamics simulations were recently carried out on the LysU dimer with ATP and lysine present at both active sites (Hughes et al, manuscript in preparation). These calculations represent the first examples of molecular dynamics simulations of a dimeric aminoacyl-tRNA synthetase in which the entire dimer is fully solvated in explicit solvent and simulated without restraints. Due to the large size of aaRSs and their dimeric nature, previous simulations were focused on a single monomer, with the application of harmonic restraints to atoms outside a 20 Å radius of one active site [18-20]. In LysU, the two active sites are approximately 30 Å apart. Our unrestrained simulations showed that the enzyme develops a significant asymmetry, both in terms of local active site interactions with ATP and also in a global sense as illustrated in Figure 1(b). The motif 2 loop shown in green in Figure 1(b) (residues 261–273) and the loop between strands C2 and C3 which forms part of the dimer interface were shifted by between 1.5–2 Å, so that the motif 2 loop was pulled further over the active site in monomer 1. The C2'-C3' loop (where the prime indicates an element of monomer 2) appeared to lie directly above and move in concert with the motif 2 loop of monomer 1. Further differences were seen in the N-terminal domain which seemed to undergo rigid body rotation, as well as in the insertion domain, which showed variations in orientation of the 4-helix bundle and in the helix-loop-helix motif. Differences were also observed in the active site nonbonded interactions of the ATP ligand, with generally better van der Waals interactions at site 2. Interestingly the flipping loop (yellow in Figure 1(b), residues 215 to 217) that lies above lysine was unaffected, and the active site interactions of lysine appeared to be equivalent in the two monomers. The presence of asymmetry, and its persistence in a number of models, led us to consider the possibility that it might be genuine, and effect different binding affinities.
Here, the reported 1 nanosecond simulation above is termed the fully occupied (FO) model, since ATP is present in both active sites. The initial structure was obtained from the published coordinates of LysU with lysine and ATP analogue AMPPCP present in both active site pockets . Three magnesium ions are associated with the ligands in each cleft. For the investigation of cooperativity described here, two half-of-the-sites (HS) models were simulated in addition for 750 ps. The ternary complex of LysU with lysine and AMPPCP was again used as starting structure however the nucleotide and its associated magnesium ions were removed from the second site to form the HS1 model, and from the first site, to construct the HS2 model. The validity of each new model was verified by calculating the average backbone RMS deviation (1.1 Å in all cases) and analysing the ligand-protein interactions to ensure active site integrity (data not shown).
Isothermal titration calorimetry (ITC), intrinsic tryptophan fluorescence and circular dichroism (CD) techniques were employed to study the ATP and lysine ligand binding events.
Fluorescence Studies of Ligand Binding
LysU has two tryptophan residues per monomer. W115 is in a loop between helix 4 and β-strand A4 in the N-terminal OB domain, W364 is in a loop between helix 12 and helix 13 of the insertion domain 4 helix bundle. These residues were used as intrinsic fluorescent probes to study structural changes on ligand binding. Results are shown in Figure 2, and thermodynamic parameters are reported in Table 1. All experiments were performed with 10 mM MgCl2 in 25 mM TRIS-HCl pH 8.0 at 20°C, excitation was at 295 nm in order to prevent interference from tyrosine residues, and emission was observed at 340 nm. It was not possible to easily attribute particular fluorescence behaviour to particular tryptophan residues, so gross change in fluorescence emission of the enzyme at 340 nm was used as a probe of structural changes. The fluorescence intensity of the tryptophans increased by 14% when LysU was saturated with lysine, however no fluorescence change was observed upon the addition of ATP or AMPPCP to either LysU or the LysU-lysine complex. This suggested that whereas a global structural change may occur on lysine binding, only local structural changes resulted from nucleotide binding. Therefore tryptophan fluorescence was only considered suitable for observing the lysine binding event. Importantly, no fluorescence change was observed on lysine titration in the presence of 20 mM EDTA, suggesting that magnesium ions were required for lysine binding. Magnesium ion titration alone did not affect intrinsic tryptophan fluorescence. Two types of experiments were performed, firstly, the dissociation constant was elucidated under Michaelis-Menten conditions of low enzyme concentration and results were fit using a simple single site model (Figure 2, panel a). Second, the Klotz tangent method was used to calculate the stoichiometry by using an enzyme concentration significantly above that of the value of the binding constant (Figure 2, panel b). The dissociation constant of lysine binding to LysU was measured to be 7.96 ± 1.40 μM with a stoichiometry of 0.95 ± 0.1 lysine per monomer (Table 1). No allosteric effects of lysine binding were observed. Similar results were observed on addition of 100 μM ATP prior to the titration of lysine, suggesting that lysine binding preceded nucleotide binding.
Figure 2. Fluorimetric titration of LysU (a) Fluorimetric titration of LysU with lysine observing change in fluorescence under Michaelis-Menten conditions. The enzyme is at a concentration (0.5 μM) significantly below the dissociation constant. The continuous line is a simple Michaelis-Menten single site fit to the dataset. (b) Fluorimetric titration of LysU with lysine under stoichiometric conditions. The enzyme is at a concentration (17 μM) significantly above the value of the dissociation constant. The two straight lines represent fits on the linear regions of their data, the intercept results in the position of the Klotz transition represented by the dotted line.
Table 1. Experimentally determined thermodynamic parameters of binding events.
Isothermal Titration Calorimetry Studies of Ligand Binding
In contrast to fluorescence studies, ITC proved applicable to the study of both lysine and nucleotide binding. Experiments were carried out in 25 mM TRIS-HCl pH 8.0 at 20°C. Both binding events were shown to be exothermic, the downward peaks in Figure 3 indicating heat release on injection of ligand into the enzyme solution. The binding of lysine to LysU is shown in Figure 3, panels a and b. The dissociation constant of the LysU-lysine complex was found to be within experimental error of that measured by fluorescence (Table 1), and the calculated stoichiometry confirmed that binding of lysine occurs at both active sites of the dimer. No lysine binding took place if MgCl2 was replaced by 20 mM EDTA, thus confirming the requirement for the metal ion. ATP was found to be inappropriate as a ligand to study the nucleotide binding event, owing to ATP hydrolysis and Ap4A synthesis taking place and hence complicating the calorimetric data. Of the two analogues AMPCPP and AMPPCP tested in place of ATP, only AMPPCP was found to bind tightly to LysU. The latter has also been used to solve a crystal structure of the LysU•lysine•ATP ternary complex, and has been shown to bind in much the same manner as ATP. The results of AMPPCP binding are shown in Figure 3, panels c and d and Table 1; crucially the nucleotide was found to display half-of-the-sites behaviour by only binding at one of the active sites of the dimer.
Figure 3. ITC data of ligand binding to LysU (a) The lower trace shows a calorimetric titration profile of aliquots of lysine being injected into 6 μM LysU (concentration of monomer) in the presence of 10 mM MgCl2 in 25 mM TRIS-HCl pH 8.0 at 20°C. The upper trace is the control of 225 μM lysine being injected into the same buffer under identical conditions in the absence of LysU. (b) A least squares fit of the data, the heat absorbed per mole of titrant versus the ratio of concentration of lysine to the monomeric concentration of LysU. (c) The lower trace is a calorimetric titration profile of aliquots of 250 μM AMPPCP being injected into 20 μM LysU (concentration of monomer) in the presence of 10 mM MgCl2 and 1 mM lysine in 25 mM TRIS-HCl pH 8.0 at 20°C. The upper trace is the same titration in the absence of lysine. (d) A least squares fit of the data, the heat absorbed per mole of titrant versus the ratio of concentration of AMPPCP to the monomeric concentration of LysU.
Nucleotide binding was found to only take place in the presence of lysine (Figure 3, panel c) suggesting an ordered mechanism of substrate addition i.e., lysine in both active sites followed by nucleotide in a single active site.
CD Studies of Ligand Binding
The differential absorption of circularly polarised light was used as a probe to study the conformational changes of LysU on ligand binding (Figure 4). It can be clearly seen that amino acid binding caused a significant change in the tertiary structure of the enzyme, whilst magnesium ion and nucleotide binding events caused far smaller structural changes to the global conformation of the protein. These studies also confirmed the order of binding, lysine binding caused a large structural change even without nucleotide, whilst a smaller nucleotide-induced structural change was only observed in the presence of lysine.
Figure 4. CD studies of LysU in a variety of ligand bound states. The continuous black line (A) shows the ellipticity of 2 μM LysU without ligands; the red dashed line (B) shows 2 μM LysU + 10 mM MgCl2; the green dashed-dotted line (C) shows 2 μM LysU + 10 mM MgCl2 + 1 mM lysine; the blue dotted line (D) shows 2 μM LysU + 10 mM MgCl2 + 1 mM lysine + 250 μM AMPPCP; the brown dash dot dot line shows 2 μM LysU + 10 mM MgCl2 + 250 μM AMPPCP(E).
Overall, our combined experimental binding studies were considered consistent with this ordered mechanism of substrate binding to the dimer, with observable asymmetry developing during the nucleotide binding event:
LysU•LysU → Lysine•LysU•LysU•Lysine → AMPPCP•Lysine•LysU•LysU•Lysine
Calculation of Binding Free Energies
The computational approach adopted here was to employ molecular dynamics simulations in conjunction with free energy calculations. The Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) method [21-23] was applied to determine the free energy of binding ATP at both active sites of the FO model, and the relevant active site of the half-of-the-sites models. In this approach, energy terms of the bound and free species were calculated at intervals along the trajectory, and averaged to obtain the calculated binding free energy.
Definition of Complex, Protein and Ligand Species
Binding free energies were calculated as the difference between the free energy of the enzyme-ligand complex and the free ligand and protein. The "complex" species in each case was derived from snapshots of the solute along the simulated trajectory. The definition of the "protein" and "ligand" species was more involved. As lysine and magnesium have been shown to bind prior to the nucleotide, it was assumed that lysine and a single magnesium ion are always associated with the protein (see Materials and Methods section). The remaining two magnesium ions are treated as part of the ligand ATPMg2. The protein species was therefore defined as the entire dimeric protein with associated lysines, and in the FO model, included the ATP and three magnesium ions present in the alternative site. Average energy terms and the free energies of ATPMg2 binding to both sites of the FO model are detailed in Table 2. The definition of each of the energy terms is described in the Materials and Methods section, except for two. The polarisation and electronic energy term Epolar+elect is simply the sum of the electrostatic and solvent polarisation contributions to the overall energy. G*binding is the free energy of a molecule without the entropy contribution:
Table 2. Energy terms for binding of ATPMg2 to both sites of the FO model
The calculated binding free energies ΔGbinding of -3 and -13 kcal/mol are quoted as estimates, due to the approximate nature of the solute entropy contribution. Experimentally, the absolute binding free energy was found to be -8.068 ± 0.321 kcal/mol, and binding only occurred at one site. The estimated entropy term for the ligand, TΔS ≈ 28 ± 5 kcal/mol, was thought likely to be the main cause of the deviation from the experimental value. It is interesting to note that if the entropy contribution were 3 kcal/mol larger at 31 kcal/mol (which is within the error range), then the absolute binding free energy estimates would change to 0 kcal/mol in site 1, and -10 kcal/mol in site 2, improving correspondence with the calorimetric data. In this case a calculated ΔG*binding value of -31 kcal/mol would correspond with a zero experimental binding energy.
For the purposes of this work, the difference in binding free energies of the two active sites were of more interest, as a discrepancy was anticipated from the previous modelling study (Hughes et al, manuscript in preparation). Given that the entropy change was assumed to be the same for binding at the two sites, comparing ΔG*binding values is equivalent to comparing ΔGbinding, the absolute binding free energy difference. As the ΔG*binding values shown are averages over the trajectory, the associated errors were calculated more precisely. Importantly, significant differences in the relative binding free energies, ΔG*binding, of the two sites were observed, in agreement with the experimental findings. The statistical significance of the 10 kcal/mol binding free energy difference was established using a standard t-test.
The cause of the large binding free energy difference was further examined (Table 2). Binding at site 2 occured more strongly, in agreement with the structural analysis of nonbonded interactions in the modelling study. The discrepancy in binding affinity was found to be largely due to differences in van der Waals interactions (ΔEvdw) which appear to be more favourable in the second site. Contributions from the dominant electrostatic and solvent polarisation terms largely cancelled in the overall ΔG*binding. This behaviour is not unexpected, as a large effect of solvation is to screen electrostatic interactions. Kollman has noted the importance of the balance between the continuum model Poisson-Boltzmann energies and gas phase interaction energies in the accurate prediction of binding free energies with the MM-PBSA method . Likewise, the errors associated with these terms are correlated, and therefore largely cancel to reduce the overall standard error in ΔG*binding.
Time Evolution of Binding Free Energies
The variation of instantaneous ΔG*binding values along the equilibrium trajectory for ATPMg2 binding to the FO model is shown in Figure 5. While the ΔG*binding values for binding at site 1 fluctuated about roughly the same average i.e. -31.34 ± 0.69 kcal/mol, the ΔG*binding values for binding at site 2 were more varied. During the interval delimited by arrows in Figure 5 (middle panel), the binding energies deviated substantially from the average of the rest of the data. Throughout this period, the ΔG*binding was seen to fluctuate quite significantly. In this instance, the gross averaging over all trajectory data served to increase the uncertainty in the final binding free energy of site 2, and diminished agreement with experiment.
Figure 5. Time evolution of calculated properties ΔG*binding as a function of time for each point in the MM-PBSA analysis for ATPMg2 binding to site 1 (top panel) and site 2 (middle panel) of the FO model. The solid line indicates the average binding energy at each site, in site 2 this excludes the interval delimited by arrows. Bottom panel: time evolution of the distance between Arg269' carbocation and ATP γ-phosphate. The correlation of ΔG*binding at site 2 to this interaction distance is apparent.
A potential advantage of free energy calculations over experiment is that one can find a structural cause of particular binding free energy changes or differences. The binding free energy trace shown in the middle panel of Figure 5 was considered to bear a similarity to an interaction between Arg 269' (where the prime indicates the residue is in the second monomer) and ATP that was observed in the simulations (Hughes et al, manuscript in preparation). The interaction distance is shown in Figure 5 (bottom panel). The side chain of Arg 269' is located at the entrance to the active site, next to the class conserved positively charged residue His270 and during the time period 557 ps to 696 ps, it was observed to interact variously with the γ-phosphate of ATP, thus enhancing the binding of ATP to this site. It is likely that this residue, part of the motif 2 loop, may play a role in the dissociation of the pyrophosphate and Lys-tRNALys products. Once reaction has occurred, the mobility of the loop and of the arginine side chain could expedite dissociation from the active site. The simulation appears to represent an attractive illustration of this process, although as no reaction can occur, the arginine has to leave "empty-handed".
Binding energies calculated during this interval would not correspond to a ground state enzyme-ligand complex, but rather to a quasi-transition state or quasi-product state, and therefore should not be included in the average binding free energy. When the data in the interval 557 – 660 ps were removed, the remainder of the trajectory averages at ΔG*binding = -31.73 ± 0.74 kcal/mol for binding at site 1, and -38.20 ± 0.56 kcal/mol for binding at site 2. The latter corresponds to an absolute binding free energy of ca -10 kcal/mol, which is in better agreement with experiment. Further improvement was obtained if the entropy estimate of 31 kcal/mol was used, resulting in a binding free energy of approximately -8 kcal/mol. Statistical analysis showed the two data sets to have significantly different means (at the 99% confidence level), even when the middle interval data were removed.
For consistency with experimental findings that the ATP analogue was found to bind at only one site, and to explore the possibility that the asymmetry in the model was due to negative cooperativity, two half-of-the-sites (half-sites) models were simulated. If a negatively cooperative mechanism operates, then we surmised that ATP bound strongly at the second site of the FO model would prevent strong binding at the first site. For instance if ATP were placed only in the first site, then it should be bound strongly, and make good interactions with the protein. Likewise with ATP bound at the second site. If only one ATP were bound per dimer, then whichever site were occupied should be able to adjust to the presence of ligand, and therefore make strong interactions. The calculated binding energies of the HS1 and HS2 models along the 750 ps trajectory are shown in Figure 6.
Figure 6. Time evolution of the half-sites models ΔG*binding as a function of time for each point in the MM-PBSA analysis for ATPMg2 binding to the two half-sites models. The solid lines indicate the average binding energies at site 1 (HS1 model) and site 2 (HS2 model) over the time period 400 – 750 ps.
ATP binding in the HS2 model was always strong, but ΔG*binding stabilised around -40 kcal/mol after an initial period. The binding affinity for ATP was expected to be strong in HS2, since it was constructed from the FO model, which showed stronger binding at site 2. Of more interest is the binding affinity of the HS1 model. Initially the calculated ATP binding energy increased, reaching a maximum at around 320 ps. Thereafter the binding energy decreased, and converged with the binding energy for the HS2 model around 400 ps. The average binding energies for the remainder of the trajectories were both favourable: -8 kcal/mol in site 1 and -11 kcal/mol in site 2, and were of the same magnitude as the strong binding ligand in the FO model. This result seems to corroborate an anti-cooperative model, as in the absence of ATPMg2 in the alternative site, the binding free energy can reach its full potential. The slightly lower average value of binding to the first monomer could be due to the initial unfavourable binding state, longer simulations might remove this statistically insignificant difference in average ΔG*binding values.
Further evidence for anticooperative behaviour was found by examining the FO simulation. The gas phase molecular mechanical interaction energy between the ligand and Arg269 was calculated for each snapshot along the trajectory (Figure 7). It is apparent that a decrease in the interaction energy between ATP and Arg 269 at site 1 occurred during the same time interval as the increased binding at site 2 due to the interaction between Arg269' and ATP. Seemingly a stronger interaction at one site resulted in weaker interactions at the other. This is consistent with an "alternate sites" mechanism of negative cooperativity .
Figure 7. Time evolution of the molecular mechanical interaction energies between ATP and Arg 269. A decrease in the interaction energy between ATP and Arg 269 at site 1 (black trace) is correlated to an increased interaction energy between ATP and Arg 269 at site 2 (grey trace).
Analysis of the conformational changes that occurred in the half-of-the-sites models was expected to reveal the mechanism of information transfer. In Figure 8 monomer 1 is shown in blue, while monomer two is in pink, the darker shades are used to represent the starting structure, while the final snapshot is shown in a lighter shade. The HS1 model was examined first in Figure 8(a). Site 1, which contains ATP (black), exhibited a significant movement (up to 2.3 Å) of the motif 2 loop further over ATP, with some buckling. Packed directly above the motif 2 loop of monomer 1 was the C2'-C3' loop of monomer 2. Although strands C2' and C3' maintained their position, this loop also shifted up and left by up to 2 Å, thus maintaining the nonbonded interactions between the loops. Directly above lysine, the flipping loop also shifted left and down by 1.7 Å. The nearby insertion helix-loop-helix motif appeared to move by a similar amount and direction. The 4-helix bundle of the insertion domain rotated slightly towards the active site. The flipping loop and the two insertion motifs have been shown to act in a cooperative manner to close up the active site upon lysine binding . Thus is reasonable to assume that they are acting in a cooperative manner here.
Figure 8. Conformational changes in the half-of-the-sites models. Divergent stereograms showing a snapshot of the HS1 model (a) and the HS2 model (b) superimposed upon the catalytic core of the starting structure. The view is down the twofold axis as in Figure 1(a). Monomer 1 is coloured in blue and light blue, while monomer 2 is shown in shades of pink. The darker shade is used to represent the starting structure, while the final snapshot is shown in light pink or light blue. In both cases the dimers are oriented such that the active site containing ATP (shown in ball and stick representation, in black) is the left hand subunit. Very similar conformational changes are observed for the two half-of-the-sites models. Conformational changes can be observed in the motif 2 loop (labelled M2), the C2-C3 loop (labelled C2), flipping loop, an insertion domain loop and the 4-helix bundles of both monomers (locations as given in Figure 1). In the left hand monomer, the helical bundle rotates slightly toward ATP while in the right hand monomer it rotates away from the cleft. The N-terminal domain rotates as a rigid body in both monomers. This figure was created using Molscript  and RASTER3D .
Crossing the dimer interface, the motif 2 loop of monomer 2 moved further over the active site but to a lesser extent (1.2 Å) and without buckling. The flipping loop also shifted by a similar amount, however it did not move downwards, as observed in site 1. The helix-loop-helix motif twisted substantially causing the loop to flip downward and the helical bundle moved in the opposite direction in this monomer, rotating as a unit away from the active site to a greater extent.
Both N-terminal domains moved to some extent, apparently as rigid bodies. The anticodon binding domain of one monomer was in direct contact with the catalytic domain of the other, specifically the helix-loop-helix motif. Hence, rotation of the β-barrel domain could be correlated to the structural changes around the active site of the alternate subunit; that the anticodon binding domain of monomer 1 moved to a greater extent than that of monomer 2 seemed to support this argument. Alternatively the movement could simply reflect the flexibility of the smaller terminal domain.
Very similar changes were observed in the HS2 model, illustrated in Figure 8(b), except that the monomer in which ATP is present was reversed, and so the changes described above applied to the other subunit. One difference is that the helix-loop-helix motif in the "unoccupied" monomer 1 did not twist to the same extent as in the HS1 model, however the outward rotation of the helix bundle and the movement of the flipping loop were identical. This is extremely encouraging evidence for a cooperative mechanism, as the two half-of-the-sites models undergo the same conformational changes, and resulted in equivalent binding free energies.
In this work, we have performed both modeling and biochemical experiments in order to gain a better understanding of the mechanism of the lysyl-tRNA synthetase, LysU. An important aspect is that an initial modelling study provided a hypothesis of asymmetry that was subsequently confirmed by experiment. Further modelling studies have extended this hypothesis and provided important insights into the functioning of LysU.
The utility of the MM-PBSA method in predicting association energies and understanding binding processes in LysU has been demonstrated. Post-processing free energy methods have gained in popularity due to their speed advantage over the more rigorous free energy perturbation or thermodynamic integration approaches. In fact, the traditional techniques are not applicable to the present problem – the calculation of ATP binding free energies in a dimeric enzyme. Conventional methods rely on the calculation of relative free energies of two (or more) ligands binding at the same site, and cannot be applied to the same ligand binding at two different sites. Furthermore, it is not feasible to include the entire dimer in the more rigorous approaches, due to the large size of LysU. Previous thermodynamic integration calculations on asparagine and aspartic acid binding to the dimeric class II aspartyl-tRNA synthetase (AspRS)  have been limited to a single active site, with stochastic boundary conditions preventing movement in the rest of the enzyme, and hence any cooperative behaviour. In contrast, the simulations reported here are performed on the fully solvated dimer, with no restraints. This affords greater flexibility to the modelled enzyme and permits the observation of asymmetry and cooperative behaviour.
Asparagine binding by mutant AspRS has also been studied by approximate free energy calculations, using the Poisson-Boltzmann free energy method (PBFE) . Again the study focused upon binding at a single active site. The PBFE method is similar to MM-PBSA but neglects non-electrostatic interactions to the binding free energy. As has been noted for the FO model, the electrostatic terms are very similar, and binding free energy differences arise largely from the differences in the van der Waals contributions. The omission of this relatively inexpensive term may result in incorrect ranking of binding affinities in PBFE calculations.
At the outset of this study, no applications of the MM-PBSA method to a metalloprotein had been reported, however Donini and Kollman have reported studies on the zinc-dependent matrix metalloproteinases  in which the difficulties associated with the inclusion of metal ions are highlighted. In this work, binding energies were calculated in reasonable agreement with experiment, in the presence of three magnesium ions.
An ordered binding mechanism
A combination of isothermal titration calorimetry, fluorescence and circular dichroism has revealed the order and stoichiometries of the binding events and allowed the determination of associated thermodynamic parameters. In union with previous crystallographic data on the apo and lysine bound LysS structures  as well as the ATP, AMPPCP and adenylate bound LysU structures [13,14,16], this further work provides a comprehensive description of the substrate binding mechanism. Lysine appears to bind before the nucleotide in an ordered magnesium dependent process. The requirement for an ordered mechanism is suggested by the large conformational change that accompanies lysine binding, as measured by tryptophan fluorescence and CD. LysU required the presence of lysine to crystallize , however detailed changes brought about by lysine binding have been obtained by a comparison of the apo and lysine bound LysS crystal structures . As the overall structure of LysS is very similar to that of the highly homologous LysU, a common mechanism may be presumed. Onesti et al have observed that lysine binds in an induced fit mechanism, triggering a reorientation of the four helix bundle (H10 to H13) of the insertion domain. As W364 is a residue within the rotating insertion domain, it is likely to be responsible for the fluorescence changes observed during the binding of lysine to LysU. The structural studies also show that binding of lysine imposes significant order on the flipping loop (yellow in Figure 1) and the loop between two helices H14 and H15 of the insertion domain helix-turn-helix motif (shown in purple in Figure 1). The overall effect is to close up the active site around lysine . Similar observations have been reported for the Thermus thermophilus LysRS . These changes are likely responsible for the loss of entropy as measured by isothermal calorimetry, which shows the process to be enthalpically driven. A stoichiometry of two lysines per dimer is observed by ITC in agreement with the fluorescence data. Whilst the LysU and LysS structures do not reveal any magnesium ions in their lysine-bound states, our experiments have shown that magnesium is required for lysine binding.
Binding order appears to be enzyme specific, as it is not conserved even within the LysRS , however a recent series of papers has elucidated the mechanism of the Bacillus stearothermophilus lysyl-tRNA synthetase [29-31]. As in our experiments, it was observed that lysine binding causes a significant change in fluorescence whilst ATP only causes a small change that was not quantifiable. An ordered binding mechanism was found, with two moles of amino acid binding prior to ATP.
Asymmetry in nucleotide binding
Nucleotide binding elicited no change in intrinsic tryptophan fluorescence, and by inference no large structural changes. This is in agreement with the published crystal structures of the LysU ATP or AMPPCP bound states which show the only significant effect of nucleotide binding is to impart order on the motif 2 loop . These structures also show that there is little difference between the LysU-lysine-AMPPCP and the LysU-lysine-ATP states. This therefore supports our argument that AMPPCP is a good substitute for ATP in studying the binding reaction. The same analogue has been used in the ITC experiments described here, where it was not possible to use ATP due to hydrolysis and Ap4A synthesis under the same conditions. Calorimetry established the nucleotide binding process to be both enthalpically and entropically favourable, but perhaps of more interest is the stoichiometry of the nucleotide binding which was found to bind in the ratio of 1 per dimer. This important result was not anticipated from any previous biochemical data of LysU, but is compatible with the asymmetry observed in molecular dynamics simulations (Hughes et al, manuscript in preparation).
Other class II synthetases display asymmetric binding behaviour or half-of-the-sites activity. In the B. stearothermophilus LysRS only one mole of aminoacyl adenylate is formed per dimer, despite the presence of two moles of lysine. The functional asymmetry is therefore acquired after amino acid binding. Yeast LysRS, another homodimer, also exhibits half-of-the-sites binding with respect to aminoacyl adenylate as well as tRNALys . AspRS, from bakers yeast, shows cooperative behaviour, and has two different ATP binding constants that differ by several orders of magnitude (μM and mM) . No evidence for asymmetry or differential binding behaviour can be found from crystallographic studies of the lysyl-tRNA synthetase. This may be due to the effects of crystal packing, or simply that under conditions of high ATP concentrations, binding at the weaker ATP site forces symmetry on the molecule that is not relevant to its in vivo mechanism.
Molecular dynamics based free energy calculations also highlight the different binding characters of the two monomers. Structures from the previous simulations (the FO model) have been the subject of free energy calculations described here. The MM-PBSA method has been used to calculate ATP binding free energies along the trajectory.
A clear difference in calculated binding affinities has been observed with only one site binding strongly. The absolute free energies reported here are approximate due to the uncertainty in the solute entropy term. A further compounding feature is that the experimental value (-8.068 ± 0.321 kcal/mol) was measured for AMPPCP and not ATP. It has already been established that AMPPCP is a good analogue of the bound ATP hence binding affinities are expected to be in the same range. Supporting this assumption, a dissociation constant corresponding to a ΔGbindingof ca -7 kcal/mol for ATP association to B. stearothermophilus LysRS has been reported .
The large difference between the calculated binding affinities (ca 10 kcal/mol) has been examined by comparing ΔG*binding values which are independent of the solute entropy and hence more accurate. The difference in binding energies between the two sites of the FO model is statistically and biologically significant. Binding at the second site was found to be more favourable. In conjunction with the experimental study, these findings led us to propose that LysU may operate in a negatively cooperative manner.
Several mechanisms have been proposed to explain half-of-the-sites activity or strong negative cooperativity . The choice of functional active site can be random, or guided. In the flip-flop model, proposed by Lasdunski, binding at the first subunit reduces the affinity of the second subunit for the same ligand . Catalysis at the first site revives the binding affinity of the second site. The binding energy of the second ligand can be used to facilitate further catalysis at the first site, or help in product release. The "alternate site" model described by Boyers is a more general model with the same principle of ligand-induced subunit cooperation being an integral part of the catalytic process . Alternatively, a ligand-induced conformational change could lead to repeated utilisation of only one site .
The molecular dynamics studies have been carried out with ATP present in both active sites, and indicate only a preference for one site, furthermore the preferred binding region appears to be random. To further explore the possibility of cooperative behaviour, two half-of-the-sites models have been simulated, by removing ATP from a single monomer of the FO model. The binding energies of both models converge on the same mean despite the HS1 model being the initially weak site, as it was in the FO simulation. Similarity between this mean and the affinity of the strong binding site of the FO model indicates that the full binding potential can only be realised at one of the two active sites. This was seen as evidence for ligand induced anti-cooperative mechanism. Further evidence for negative cooperativity has been found in the time evolution of binding energies or gas phase interaction energies of the FO model. The mobile motif 2 loop residue Arg 269' was seen to enter the active site of the second monomer and interact with ATP for a period. During this interval, the interaction energy between ATP and Arg 269 is significantly enhanced in site 2 and weakened in site 1. Somehow, improved binding at one active site results in weaker binding at another active site 30 Å away. This behaviour is strongly reminiscent of an alternate sites model, especially in light of the role that Arg 269 is postulated to play in the removal of pyrophosphate from the active site once the first reaction has occurred.
An alternating sites mechanism, and half-of-the-sites activity, requires a dimeric (or tetrameric) configuration. In general, the class I synthetases are monomeric while most class II enzymes are dimers, heterotetramers or tetramers . Two of the three conserved motifs in class II enzymes (motifs 1 and 2) are located at the dimer interface, and dimerisation was shown to be required for activity in ProRS and AspRS . Several dimeric class I synthetases also show evidence of half-of-the-sites activity in the adenylation reaction. Perhaps the earliest and best known example is the tyrosyl-tRNA synthetase (TyrRS) from B. stearothermophilus, which binds one mol of tyrosine and forms only one mol of adenylate per dimer . A second mol of tyrosine and ATP only binds after the tyrosyl adenylate has formed, hence an interacting sites mechanism was proposed for TyrRS activity . The crystal structure of TyrRS is symmetric  and tyrosine is present in both sites of the E:Tyr complex . To resolve this apparent discrepancy, Ward and Fersht created truncated mutant heterodimers, which showed that TyrRS is asymmetric in solution, and that this asymmetry is frozen-in and randomly distributed between active sites with no detectable interconversion over several minutes . It was concluded that the enzyme asymmetry, rather than negative cooperativity, was responsible for the observed half-of-the-sites activity . Given the nature of the mutations, which involved removing an entire anticodon binding domain of the dimer, it is plausible that the lack of interchange between sites is due to the removal of residues involved in a cooperativity pathway.
The closely related tryptophanyl-tRNA synthetase from beef pancreas is negatively cooperative in the formation of the enzyme-adenylate complex . The second mol of adenylate was shown to form two orders of magnitude more slowly. In the presence of pyrophosphatase (or tRNATrp), two mol of adenylate are formed without any apparent cooperativity. It was concluded that the anticooperativity was due to the presence of pyrophosphate in the first active site. This was an interesting finding, as it implied that asymmetry could be "trapped" in an alternating sites mechanism. Many kinetics experiments on aaRS are carried out in the presence of inorganic pyrophosphatase or tRNA, and would therefore not show negative cooperativity. In our ITC experiments on LysU a non-hydrolysable ATP analogue was used, allowing the hidden asymmetry present in the mechanism to be uncovered. Similarly in the molecular dynamics simulations, no reaction can occur and the pyrophosphate cannot diffuse from the active site. In a recent paper, Retailleau et al conclude that there is strong crystallographic evidence that the B.s. TrpRS is an asymmetric dimer under most circumstances and suggest that a second ATP binds with a much lower affinity .
Possible functional roles of negative cooperativity and a postulated alternating sites mechanism
Functional asymmetry is a fairly common feature of the dimeric aaRSs, therefore it is not unreasonable to assume that it might play an important role. Several reasons have been proposed for active sites to interact in a negatively cooperative manner. One advantage would be to increase reactivity at one subunit  or to regulate substrate binding in multimeric proteins . Another function proposed is more specific to aminoacyl-tRNA synthetases: interacting active sites might act as a safeguard against errors in genetic translation [24,41]. Qiu et al have suggested that cooperativity may be necessary for efficient product dissociation and release in aaRS .
We propose that LysU might operate in the following alternating sites mechanism: initial ATP binding to the first site causes a closing up of that cleft and prevents strong binding at the second site. Once adenylate formation and pyrophosphate dissociation has occurred at site one, ATP binding can occur strongly at site 2. In this mechanism, a conformational change or release of strain caused by pyrophosphate dissociation from one monomer is transmitted across the dimer interface to allow strong binding at the other monomer. Reaction followed by release of pyrophosphate at this site in turn revives the ATP binding affinity of the first. Such a mechanism would account for the half-of-sites binding observed by experiment, and is also consistent with the negative cooperativity mechanism proposed for TrpRS . It is also relevant that in the molecular dynamics simulations Arg 269 only enters the strongly bound ATP site of the FO model. The simulated models were examined to see how this information could be transferred. The most direct pathway between the two active sites is through the motif 2 residues Arg 262 and Phe 261 of one monomer across the dimer interface to their symmetry equivalent Phe 261' and Arg 262' however given the nature of the conformational changes observed in the half-sites models a more distributed pathway can be proposed. The movement of the motif 2 loop, flipping loop and insertion domain elements in monomer 1 (illustrated in Figure 8) all serve to close up the active site slightly. This suggests some release of strain upon the removal of ATP from the second site. In closing up the reaction cleft of the first monomer, a shift of the motif 2 loop could be transmitted via the C2'-C3' loop, to strand C2'. From here, the signal could be transmitted across to strand C2, and hence to the C2-C3 loop above the motif 2 loop of the second monomer, as well as to the flipping loop of the second monomer, which is just downstream of strand C2'. The flipping loop is in close proximity to the twisted helix-loop-helix motif. This causes the helix bundle to rotate away from the active site in a concerted manner. In this way, the closing of one active site around ATP is coupled to the opening up of the site in the other subunit, and hence a reduced binding affinity.
Such a mechanism could improve catalytic efficiency in two ways. The same conformational changes that open and close the active sites could prime a particular monomer for tRNA binding and release. Concerted conformational changes also ensure the correct order of substrate binding and reaction. It is apparent that cooperative structural dynamics are an important feature of aminoacyl-tRNA synthetase activity.
Previous molecular dynamics simulations suggested an asymmetry in nucleotide binding to LysU. Further biochemical and computer experiments described here have elucidated a functional asymmetry in this homodimer. A functional asymmetry appears to be a common feature of many oligomeric aminoacyl-tRNA synthetases and may play an important role in increasing catalytic efficiency. We have suggested a manner in which this could be achieved by LysU operating in an alternating sites mechanism.
Materials and Methods
Adenosine 5'-triphosphate disodium salt (ATP), α,β-methyleneadenosine 5'-triphosphate lithium salt (AMPCPP), β,γ-methyleneadenosine 5'-triphosphate sodium salt (AMPPCP), and lysine were all purchased from Sigma-Aldrich. Purity of nucleotides were confirmed as >95% on a Mono-Q HR 5/5 FPLC column using a 0 to 1 M sodium chloride gradient (Pharmacia). Magnesium chloride was purchased from BDH Laboratory supplies. TRIS buffers were purchased from Anachem Ltd.. All solutions were prepared with metal-free deionised water which had been distilled and passed through a MilliQ system (Millipore).
LysU Expression and Purification
A new, controllable expression vector for lysU named pADH2 was designed. The lysU gene was amplified using pXlys5 as a template, incorporating SacI and HindIII sites onto the 5' and 3' ends of the gene respectively. A manipulation vector, pADH1, was constructed from the lysU-PCR product incorporation into pGEM3, then finally the lysU gene was cloned into pTrc99A, an IPTG-inducible vector. The resultant vector, named pADH2, was transformed into E. coli lysU deletion strain lysU2-17A . Purification of the overexpressed protein was as previously described . Yields of LysU were 160 mg/l bacterial culture. Purity of LysU was estimated as >95% by SDS-PAGE.
Fluorescence Measurement of the LysRS-ligand Interaction
The interaction between LysU and ligands was studied using a Shimadzu RF-5301 PC spectrofluorimeter. Measurements were made at 20°C in a stirred 700 μl volume cuvette. Tryptophan fluorescence was excited at 295 nm (5 nm bandwidth) and measured at 340 nm (5 nm bandwidth). Titrations were performed adding 1 μl aliquots of ligand solution sequentially. Binding constants were calculated using enzyme concentrations of 0.5 μM, stoichiometry values were calculated using enzyme concentrations of 20 μM (concentrations are of monomer). Linear least squares regression fitting were fitted using the Origin 5.0 data analysis program, using a simple single site Michaelis-Menten model. Circular Dichroism (CD) Spectroscopy. CD spectra were observed using a Jasco J-715 spectropolarimeter and a 1 mm pathlength cell. Enzyme concentrations were 1 μM, volume change on ligand addition were less than 0.5%. Measurements were made at 20°C in 10 mM TRIS-HCl pH 8.0 using a 1 nm bandwidth and 4 second response time. Spectra were obtained as an accumulation of 4 runs.
Isothermal Titration Calorimetry (ITC)
The heat released on LysU interaction with various ligands was measured using a Microcal VP-ITC calorimeter, and data was fitted using the Microcal Origin software. The reaction cell was maintained at 20 ± 0.2 °C during all titrations. Typically enzyme concentrations were 15–20 μM (of monomer), and ligand concentrations 200–250 μM, and were in 25 mM TRIS-HCl pH 8.0 and 10 mM MgCl2. Titrations consisted of a 1 μl injection followed by either 26 10 μl injections or 54 5 μl injections. 1 μl injections were added over 2 seconds, 5 μl injections over 10 seconds, and 10 μl injections over 20 seconds. In all cases subsequent additions took place at 4 minute time intervals. Titration results were corrected for heat of ligand dilution.
Molecular Dynamics Simulations
Development of the FO model will be described elsewhere (Hughes et al, manuscript in preparation) In brief, the X-ray crystal structure of LysU from E. coli in a complex with lysine and the inhibitor AMPPCP (PDB accession code 1E22) was used as the starting structure for the molecular dynamics simulations . A dimer was created from the monomeric asymmetric unit by application of the 2-fold symmetry transformation, and AMPPCP was replaced by ATP. Thus the complex consists of two monomers, containing the ligands lysine, ATP and three associated magnesium ions in each active site. The original monomer is referred to as monomer 1, and the generated monomer as monomer 2. Similarly, the active sites of monomers 1 and 2 are referred to as site 1 and site 2. Residues 1–10 of the N-terminal domain, and residues 154–160 corresponding to a flexible region between the N-terminal and C-terminal domains were not resolved in the crystal structure due to the absence of electron density, these regions have been omitted from the model. The complex was charge neutralised by the addition of 33 sodium counterions, and then solvated by a cubic box of TIP3 waters extending at least 7Å in each direction from the solute and counterions, to create a total of 99 232 atoms. The sander module of the AMBER6 software suite  employing the Cornell et al force field  was used for minimisation and dynamics. The system was minimised with 1000 cycles of steepest descent and conjugate gradients minimisation, followed by heating to 300 K over 15 ps and equilibration at 300 K for 105 ps. Production dynamics followed for 885 ps, giving a total simulation time of 1005 ps. Periodic boundary conditions and a 12 Å cutoff for nonbonded van der Waals interactions were applied. The particle mesh Ewald method was used to treat long-range electrostatic interactions during dynamics. A 1.5 fs timestep was used with constant pressure and temperature conditions. The temperature and pressure coupling parameters were 1.0 and 2.0 ps respectively. Snapshots were saved every 0.15 ps during dynamics. This model was labeled the fully occupied FO model, as ATP was present in both active sites. In addition to the full dimer, two half-of-the-sites (HS) models of LysU were constructed and simulated using an identical protocol, for 750 ps. The lysine ligand and Mg 1 were present at both active sites, with ATP, Mg 2 and Mg 3 present only in the first monomer of the HS1 model and only in the second monomer of the HS2 model. To demonstrate that asymmetry is not an artefact of the simulation a fourth model, FOrev, which displays a reversed binding site preference was also simulated, by randomly changing the initial solvation conditions.
Binding Free Energy Calculations
Binding free energies ΔGbinding were estimated using the MM-PBSA method. The procedure is well described in the literature [21-23], and summarised here. Binding free energies are calculated as the difference in the free energies of the enzyme-ligand complex and the unliganded protein and free ligand:
ΔGbinding = Gcomplex - (Gprotein + Gligand)
The free energy of each species is obtained from a combination of "gas phase" molecular mechanics energy with solvation free energy and solute entropy terms:
Gmol = EMM + Gsolvation - TSsolute
EMM = Eelectrostatic + Evdw + Einternal
Gsolvation = Gpolar + Gnonpolar
EMM is the "gas phase" molecular mechanical energy, which has contributions from electrostatic interactions (Eelectrostatic) van der Waals interactions (Evdw) and internal degrees of freedom (Einternal). The solvation free energy Gsolvation has a solvent polarization component Gpolar and a cavity term Gnonpolar.
In this "post-processing" free energy method, energy terms are calculated as averages over solute conformations obtained as snapshots from molecular dynamics simulations performed in explicit solvent, after removal of solvent and counterions. These degrees of freedom are implicitly incorporated through the use of continuum models, thereby reducing computational effort . Averaging over many solute conformations improves the accuracy of the relative binding energies . The most difficult term to calculate is the rotational, vibrational and translational contribution to the binding free energy TΔS, typically estimated by normal mode analysis in the harmonic approximation. Due to the computational cost this is typically only performed on a few snapshots hence a larger uncertainty is associated with this term. Because of this, and the approximations inherent in the solute entropy term, absolute binding free energies calculated by the MM-PBSA method are not as accurate as those from the traditional thermodynamic perturbation or integration approaches. Despite this, good agreement with experiment has been reported for the association of biotin with avidin (ΔGbinding within 3 kcal/mol of the experimental value) , and the oncoprotein Mdm2 with the N-terminal stretch of the p53 tumor suppressor protein (within 5 kcal/mol of experiment) . Surprisingly good agreement with the experimental binding free energy (within 1.6 kcal/mol) was obtained for the binding of an inhibitor to HIV-1 RT .
The MM-PBSA method has also been successfully to predict relative binding free energies quickly, allowing a series of ligands to be ranked in terms of their binding affinity to a common receptor [50-52]. In such cases, if a similar magnitude of entropy change is expected for similar sized ligands, the free energies can be compared in the absence of the TΔS entropy term.
The present application involves the calculation of both absolute and relative ATP binding free energies in a dimeric enzyme. Our fluorescence studies showed that magnesium is required for lysine binding, but we were unable to ascertain the number of magnesium ions present. It was anticipated that calculated absolute binding free energies might enable us to establish the number of magnesium ions associated with ATP during the binding step. The more accurate relative binding free energies of the two sites can subsequently be compared without the entropy contribution to determine if a significant difference is observed.
For the calculation of absolute binding free energies, two possibilities for the ligand were considered. Test MM-PBSA calculations were carried out for the binding of ATP:Mg2 (denoted ATPMg) and ATP:Mg2:Mg3 (denoted ATPMg2). Mg 1 is situated deep within the cleft, close to the lysine ligand, and was thus considered the most likely to be required for the lysine binding step observed by fluorescence experiments. It is also the metal ion with the highest occupancy in the reported crystal structure. Mg 3 displays the next highest occupancy and penetration within the cleft. Mg 2 was expected to be the least strongly bound, as it is eliminated on formation of the pyrophosphate, and was observed to move the most of all magnesium ions in the molecular dynamics simulation. These tests showed that using ATPMg2 as a ligand resulted in binding free energies of the correct magnitude, whereas binding free energies of ATP with a single magnesium ion were 20–30 kcal/mol larger than experiment. Thus ATPMg2 was used as the definition of ligand for all the reported MM-PBSA calculations.
For the fully occupied model FO the corresponding "protein" species were the ternary complexes with lysine and the remaining magnesium ion Mg 1 present in the site under consideration while the alternative site contained lysine, ATP and three magnesium ions. Two approaches to generating structures of the free protein and ligand species were possible. Conformations could be obtained from individual trajectories of the three species or from a single trajectory of the complex. The former approach is expensive, requiring at least twice the simulation time, and results in greater standard errors. The latter approach is cheaper and results in smaller standard errors due to correlation of the errors from individual species. In selecting an approach, careful thought was given to the nature of the binding process. ATP binding has been shown not to cause any large conformational changes, thus the use of a single trajectory to describe the enzyme-ligand complex and the unliganded "protein" should be valid. The same would not hold true of lysine binding, hence this process was not submitted to MM-PBSA analysis.
Snapshots were selected at 4.2 ps intervals from the molecular dynamics simulations of the complexed enzymes, resulting in 200 structures for the FO model, and 150 for the half-sites models. EMM the molecular mechanical energy was calculated with the anal module of the AMBER molecular modelling software package. The Cornell et al forcefield  was applied as for the dynamics simulations, however infinite cutoffs were used for all interactions. The polarisation contribution to the solvation free energy, Gpolar was evaluated by continuum electrostatics, using the DelPhiII program , which solves the standard linearised form of the Poisson-Boltzmann equation by finite difference techniques (FDPB). A grid spacing of 2 grids/Å that extended by 10% in all directions beyond the dimensions of the solute was used.
The choice of solute dielectric constant (εint) is largely empirical . For proteins, low values between 1 and 4 are typically chosen to reflect the lack of polarisability of the protein environment, although higher values might be more appropriate for the study of certain processes . In this work, the solute dielectric was chosen as 1 for consistency with the AMBER forcefield and the exterior solvent dielectric constant was 80. Test calculations with higher values of solute dielectric (εint = 2–4) showed poor agreement with experiment for both ATPMg and ATPMg2, thus supporting our use of εint = 1. Tests also showed that 2000 iterations were sufficient for convergence to within 0.05 kcal/mol. Cornell et al charges, and scaled AMBER van der Waals radii were used for all atoms. The nonpolar contribution to the solvation free energy was calculated as a linear dependent term in surface area:
Gnp = γSA + b
The empirically determined parameters γ and b are 0.00542 kcal/Å2 and 0.92 kcal/mol respectively. SA is the solvent accessible surface area calculated with the program MSMS .
The entropy term was the most difficult to estimate. Normal mode analysis in the harmonic approximation was used to estimate the entropy loss during complex formation. The structure of the first monomer was obtained from a single snapshot taken at 435 ps. After removal of solvent and counterions this conformation was used to create complex, protein and ligand species for separate normal mode analyses. The structures were first energy minimised to a gradient of 10-3 kcal/mol.Å before normal mode analyses with the nmode module from AMBER6. A dielectric constant of 4rij (where r is the distance between the ith and jth atom) was used to mimic solvent screening. This procedure was performed only once for each ligand (ATPMg2 and ATPMg), due to the prohibitive memory and time requirements of the calculation. The error associated with a single calculation cannot be determined but is estimated at less than 5 kcal/mol. It is also assumed that the entropy change upon ATP binding is the same for both monomers, within this error range.
Mg2+ Radius Parameterisation
A potential difficulty with the solvation calculations was presented by the metal ion. Solvation energies are extremely sensitive to the metal ion radius. Test calculations carried out using the van der Waals radius for the Mg ion in the FDPB calculations resulted in incorrect positive binding free energies. Kollman has noted that for divalent metal ions "the use of van der Waals radii from molecular mechanics can lead to a significantly too positive calculated ΔGbinding presumably because continuum models overestimate the desolvation penalty of moving a divalent ion from a high to a low dielectric medium". It was clear that a suitable solvation radius needed to be derived through some sort of parameterisation procedure.
A logical approach was to find the radius that allows the FDPB method to reproduce experimental data for aqueous solvation of the magnesium ion. A series of solvation free energy calculations were carried out on a single magnesium ion. The solvent polarisation energy was combined with a cavity term (ΔGnonpolar) to compare with the experimental hydration free energy of -455.9 kcal/mol . The best agreement with experiment was obtained with a radius of 1.436 Å (Table 3). Using a simple Born model of solvation and the same hydration free energy yields a calculated radius of 1.44 Å in good accord with the PB calculations. This value was used for all Poisson-Boltzmann calculations.
Table 3. Calculated hydration free energies for various Mg2+ radii
Analysis of Statistical Errors
The largest statistical error is that associated with the entropy term, and is estimated to be 5 kcal/mol. Therefore the absolute binding free energies calculated in this work can only be considered to be estimates. Fluctuations in the electrostatic (Eelectrostatic) and solvent polarisation (Gpolar) terms are the next largest source of statistical error, however the effect of solvation is largely to screen the gas phase electrostatic interactions, so that the standard error in the sum of these terms Epolar+elect is much smaller. This cancellation of errors is important for the accurate prediction of binding free energies. Standard t-tests were used to establish statistical significance of the binding free energy differences. To account for possible statistical inefficiency, the data have also been averaged in blocks of length 29.4 ps to estimate the true standard error of the mean. These standard errors were used for the t-tests. Standard errors quoted elsewhere are not block averaged.
SJH carried out the molecular dynamics simulations and free energy calculations, while JAT designed and carried out the biochemical experiments. ADH designed the controllable expression vector for lysU. ADM and IRG provided intellectual input, mentorship and financial support. SJH and JAT drafted the manuscript; all authors read and approved the final manuscript.
aaRS, aminoacyl-tRNA synthetase; LysRS, AspRS, ProRS, TyrRS, TrpRS: specific aminoacyl-tRNA synthetases are referred to using their three letter amino acid code, followed by RS; ATP, adenosine triphosphate; AMPPCP, β,γ-methyladenosine 5' triphosphate; AMPCPP, α,β-methyladenosine 5' triphosphate; Ap4A, 5',5"-diadenosine polyphosphate; ITC, isothermal titration calorimetry; CD, circular dichroism; MM-PBSA, molecular mechanics Poisson-Boltzmann ; FDPB, finite difference Poisson-Boltzmann.
We thank Jonathan Swinton of Proteom for help with the statistical analysis, and the Mitsubishi Chemical Corporation for their continued support of the Genetic Therapies Centre. The authors acknowledge the EPSRC for funding support for JAT, the Department of Chemistry, Imperial College London and the Mitsubishi Chemical Corporation for funding support for SJH.
Yaremchuk A, Tukalo M, Grotli M, Cusack S: A succession of substrate induced conformational changes ensures the amino acid specificity of Thermus thermophilus prolyl-tRNA synthetase: comparison with histidyl-tRNA synthetase.
Nucleic Acids Research 1991, 19:3489-3498. PubMed Abstract
Proc Nat Acad Sci USA 1983, 80:7496-7500. PubMed Abstract
Journal of Physical Chemistry B 1997, 101:8349-8362. Publisher Full Text
Journal of the American Chemical Society 2000, 122:3909-3916. Publisher Full Text
Journal of the American Chemical Society 1999, 121:8133-8143. Publisher Full Text
Perspectives in Drug Discovery and Design 2000, 18:113-135. Publisher Full Text
Archontis G, Simonson T, Karplus M: Binding free energies and free energy components from molecular dynamics and Poisson-Boltzmann calculations. Application to amino acid recognition by aspartyl-tRNA synthetase.
Cusack S, Yaremchuk A, Tukalo M: The crystal structures of T. thermophilus lysyl-tRNA synthetase complexed with E. coli tRNA(Lys) and a T. thermophilus tRNA(Lys) transcript: anticodon recognition and conformational changes upon binding of a lysyl-adenylate analogue.
EMBO Journal 1996, 15:6321-6334. PubMed Abstract
Takita T, Ohkubo Y, Shima H, Muto T, Shimizu N, Sukata T, Ito H, Saito Y, Inouye K, Hiromi K, Tonomura B: Lysyl-tRNA synthetase from Bacillus stearothermophilus. Purification, and fluorometric and kinetic analysis of the binding of substrates, L-lysine and ATP [published erratum appears in J. Biochem. (Tokyo) 1998 Jun;123(6):1218].
J Biochem (Tokyo) 1996, 119:680-689. PubMed Abstract
Takita T, Hashimoto S, Ohkubo Y, Muto T, Shimizu N, Sukata T, Inouye K, Hiromi K, Tonomura B: Lysyl-tRNA synthetase from Bacillus stearothermophilus. Formation and isolation of an enzyme-lysyladenylate complex and its analogue.
J Biochem (Tokyo) 1997, 121:244-250. PubMed Abstract
J Biochem (Tokyo) 1998, 124:45-50. PubMed Abstract
Kern D, Lorber B, Boulanger Y, Giege R: A peculiar property of aspartyl-transfer RNA synthetase from bakers yeast – chemical modification of the protein by the enzymatically synthesized aminoacyl adenylate.
Journal of Molecular Evolution 1995, 40:499-508. PubMed Abstract
Biochemistry 1975, 14:5-12. PubMed Abstract
Biochemistry 1975, 14:13-18. PubMed Abstract
Journal of Molecular Biology 1987, 194:287-297. PubMed Abstract
Biochemistry 1987, 26:8031-8037. PubMed Abstract
Biochemistry 1988, 27:1041-1049. PubMed Abstract
Retailleau P, Huang X, Yin YH, Hu M, Weinreb V, Vachette P, Vonrhein C, Bricogne G, Roversi P, Ilyin V, Carter CW: Interconversion of ATP binding and conformational free energies by tryptophanyl-tRNA synthetase: Structures of ATP bound to open and closed, pre-transition-state conformations.
Biochemistry 1993, 32:13658-13663. PubMed Abstract
Pearlman DA, Case DA, Caldwell JW, Ross WS, Cheatham TE, Debolt S, Ferguson D, Seibel G, Kollman P: Amber, a Package of Computer-Programs for Applying Molecular Mechanics, Normal-Mode Analysis, Molecular-Dynamics and Free-Energy Calculations to Simulate the Structural and Energetic Properties of Molecules.
Computer Physics Communications 1995, 91:1-41. Publisher Full Text
Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA: A 2nd Generation Force-Field for the Simulation of Proteins, Nucleic-Acids, and Organic-Molecules.
Kollman PA, Massova I, Reyes C, Kuhn B, Huo S, Chong L, Lee M, Lee T, Duan Y, Wang W, Donini O, Cieplak P, Srinivasan J, Case DA, Cheatham TE 3rd: Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models.
Wang J, Morin P, Wang W, Kollman PA: Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA.
Kuhn B, Kollman PA: Binding of a diverse set of ligands to avidin and streptavidin: an accurate quantitative prediction of their relative affinities by a combination of molecular mechanics and continuum solvent models.
Berthet-Colominas C, Seignovert L, Hartlein M, Grotli M, Cusack S, Leberman R: The crystal structure of asparaginyl-tRNA synthetase from Thermus thermophilus and its complexes with ATP and asparaginyl-adenylate: the mechanism of discrimination between asparagine and aspartic acid.
Journal of Applied Crystallography 1991, 24:946-950. Publisher Full Text
Acta Crystallographica Section D-Biological Crystallography 1994, 50:869-873. Publisher Full Text