Identifying and engineering ion pairs in adenylate kinases. Insights from molecular dynamics simulations of thermophilic and mesophilic homologues.

Molecular dynamics simulations were performed to study thermal stabilization of proteins via electrostatic interactions of ion pairs. Dynamic motions of four ion pairs previously proposed to be important in thermal stability of adenylate kinase from the thermophile Bacillus stearothermophilus were monitored during the simulation. One of the four ion pairs identified in the crystal structure, Lys180-Asp114, was not maintained in close contact suggesting that the ion pair does not contribute to thermal stability. Among the other three ion pairs, the ion pair Arg116-Glu198 was proposed to be the most important for stability. To predict behaviors of the ion pairs when engineered into a mesophilic homologue to increase stability, in silico mutants of adenylate kinase from the mesophile Bacillus subtilis were generated, and their molecular dynamics simulations were carried out. The ion pairs in the mutant simulations displayed similar behaviors to those in the simulation of the thermophilic protein. To validate the results of the simulations experimentally, the same mutants were produced in vitro and their thermal stabilities were measured using differential scanning calorimetry. In agreement with the simulations, the Lys180-Asp114 did not result in any increase in stability by itself or additive effect with other ion pairs, whereas a mutant with the Arg116-Glu198 exhibited the highest stability among the mutants having one of the four ion pairs. These results provide specific knowledge about stability in adenylate kinases and more generally suggest that molecular dynamics simulations can provide valuable information for identifying and engineering ion pairs.

Because proteins from organisms living at high temperatures generally have increased stability, they have drawn attention for both practical and fundamental reasons. From small scale experiments to large scale industrial processes, it is often required or advantageous to use proteins that are stable and function at high temperatures, and thermophilic proteins have been used in various fields directly or as models for protein engineering (1,2). They also provide a unique opportunity to understand molecular basis of protein stability (3,4) and its relation to dynamics and function (5,6).
To illustrate the molecular mechanism of increased stability of thermophilic proteins, numerous studies have been performed by comparing them with their mesophilic homologues. There is no unifying mechanism found for the increased stability. Rather, it is apparent that different thermophilic proteins have different strategies for temperature adaptation and use combinations of subtle intramolecular interactions to gain the enhanced thermal stability (7,8). The principal interactions may include electrostatic interactions between charged residues (9), hydrogen bonding (10), and hydrophobic interactions (11).
Among these, electrostatic interactions have been of particular interest for various reasons. Although hydrogen bonding and hydrophobic interactions are relatively difficult to define, it seems easy to identify favorable electrostatic interactions in thermophilic proteins by detecting ion pairs (also called salt links or salt bridges) with a simple distance cutoff assuming that most stabilization results from attraction between two oppositely charged residues. With relatively simple models, it is also possible to find important residues in electrostatic stabilization by considering interactions with all the other charged residues, not only with oppositely charged residues within a certain distance cutoff (12,13). When engineering mesophilic proteins for increased stability, modifying the electrostatic interactions between charged residues has advantages over altering other interactions. Because charged residues are usually found on surface of proteins, but their interactions are able to contribute to overall stability, engineering electrostatic interactions by mutating surface residues can minimize unwanted effects such as disruption in structure or function, which is often found in protein mutagenesis, especially when mutating core residues. Furthermore, the decrease in the dielectric constant of water with increasing temperature enhances the importance of electrostatic interactions in thermostability.
The same simple electrostatic model used for identifying important residues in thermal stabilization of thermophilic proteins can also be used in the process of engineering mesophilic proteins (14). Using the model, residues involved in destabilizing interactions with the same charged residues can be selected and mutated to neutral or oppositely charged residues for removing destabilization or further stabilization. Although this strategy was successful in several studies (15)(16)(17), it is only applicable to the cases where unfavorable electrostatic interactions are already existing. It is more desirable to create additional favorable interactions by introducing new charged residues.
Adenylate kinases from the thermophile Bacillus stearothermophilus (AKste) 1 and the mesophile Bacillus subtilis (AKsub) are excellent targets to study protein stabilization. Although they come from organisms in the same genus and have a high level of sequence identity (74%), they exhibit significant difference in stability (18,19). Thus, the differences between them most probably result from temperature adaptation unlike other examples where pairs of proteins are evolutionarily distant and differences related to temperature may be substantially masked by random or other genetic drift. Their crystal structures have been solved and compared with each other to reveal the molecular mechanism for the increased stability of AKste (20,21). Thermophilic AKste seems to use additional favorable electrostatic interactions of charged residues, and four ion pairs in AKste were proposed as critical for stability enhancement by electrostatic attraction (21).
In this work we used thermophilic AKste and mesophilic AKsub as models to analyze and engineer protein stabilization via electrostatic interactions of ion pairs. The molecular dynamics (MD) simulations of the two adenylate kinases were performed to observe their dynamic behavior including the motions of the ion pairs. Based on the results of these MD simulations, AKsub mutants with added ion pairs were generated in silico, and their MD simulations were carried out to study ion pair engineering. Finally, the AKsub mutants were produced in vitro and tested for their thermal stability to compare the results with those from the MD simulations.

EXPERIMENTAL PROCEDURES
System Setup-The crystal structures of AKste and AKsub (Protein Data Bank entries 1ZIO and 1P3J, respectively) were used as starting models for MD simulations. The five C-terminal residues not resolved in the AKsub structure were constructed based on those from the crystal structure of the homologous Bacillus globisporus adenylate kinase (AKglo, Protein Data Bank entry 1S3G). The four crystallographic water molecules coordinating a magnesium ion in both structures were maintained in the models. The psfgen package of the program VMD (22) was used to make mutations and to add hydrogen atoms. The wild type and mutant adenylate kinases were then solvated in a 51.4 ϫ 58.2 ϫ 63.3 Å 3 box consisting of TIP3 water molecules (23) and neutralizing ions (Na ϩ and Cl Ϫ ) at physiological concentration using the solvate and autoionize packages of the VMD program (22).
Simulation Configuration-The MD simulations were carried out with the program NAMD (24) and the CHARMM27 force field (25). The force field for Ap 5 A was generated based on that for ATP. For the zinc binding site, partial charges and intramolecular force constants were taken from Eriksson et al. (26), and the average bond lengths and angles in the three crystal structures (AKste, AKsub, and AKglo) were used as the equilibrium values. The interaction between the magnesium ion and its coordinating atoms was considered as a tight electrostatic attraction, and no covalent bond or angle was defined. A cutoff distance of 12 Å for van der Waals interactions was assumed with switching function starting at 10 Å. Periodic boundary conditions were used with a flexible cell. The particle mesh Ewald summation (27) was used to compute long range electrostatic interactions without a cutoff. An integration time step was 1 fs. A constant temperature (300 or 330 K) was maintained by performing Langevin dynamics (28). The Nose-Hoover Langevin piston method (29,30) was used to keep a constant pressure (1 atm).
System Equilibration and Simulation-With the fixed crystallographic atoms, the systems were minimized for 1000 steps and equilibrated at 300 K for 20 ps. Another minimization (2000 steps) and equilibration at 300 K (100 ps) was then performed with all atoms free to move. For simulations at 330 K, the systems equilibrated at 300 K were further equilibrated at 330 K for 100 ps. A Langevin damping coefficient of 5 ps Ϫ1 was use during equilibration. The final equilibrated configurations at each temperature were used as starting points for simulations. All the simulations were then carried out for 1.2 ns with a damping coefficient of 1 ps Ϫ1 . The computation of each 1.2 ns simulation took ϳ8 days with a computer equipped with dual 2 GHz AMD Athlon processors. The coordinates were saved at every 500 time steps, and those of the last 1 ns were used for analysis.
Analysis of Simulations-All analyses were done for residues 1-212 excluding the five highly disordered C-terminal residues. The root mean square deviation (RMSD) of atomic positions was used as a numerical measure of the differences between two structures. On the other hand, the root mean square fluctuation (RMSF) is calculated as an indicator of flexibility. The RMSF is defined for each instantaneous structure in a trajectory. The reference structure used in computing RMSF is calculated by aligning and averaging all the instantaneous structures for the last 1 ns of the simulations. In each ion pair the shortest distance between oppositely charged atoms (amino nitrogen and carboxylic oxygen) was considered to be the distance of the ion pair.
In Vitro Measurement of Thermal Denaturation Midpoint (T m ) of AKsub Mutants-To produce AKsub mutant genes, site-directed mutagenesis was performed using mismatched primers in polymerase chain reactions. The mutations were verified by sequencing the entire mutant genes. As described previously (21), the mutant proteins overexpressed in Escherichia coli were purified by a two-step procedure involving affinity chromatography and gel filtration, and their T m s were measured by performing differential scanning calorimetry.

Stability of Trajectories in the Simulations-
The RMSD values between the final minimized structure and the starting crystal structure for AKste and AKsub were small, 0.37 and 0.46 Å respectively. This indicates that the force field is sufficiently accurate. The total and backbone RMSD values of the instantaneous structures during the 1.3 ns (0.1-ns equilibration and 1.2-ns simulation) from the final minimized structures are plotted in Fig. 1. The initial changes in the 300 K equilibrations were caused by heating the systems from 0 to 300 K. The changes in the 330 K were smaller and the RMSDs did not start from zero because the final equilibrated structures at 300 K were used as starting coordinates for equilibration at 330 K. All of the RMSD values seemed to stop increasing during the 100-ps equilibration, but they began to increase as the actual 1.2-ns simulations were started. This is probably because the Langevin damping coefficient was switched from 5 to 1 ps Ϫ1 . However, it is clear that the systems reached equilibrium at 300 ps after which the trajectories were saved for analysis.
The radii of gyration of AKste and AKsub were also stabilized after 300 ps. The average values for the last 1 ns of the simulations were 17.07 and 16.99 Å for AKste at 300 and 330 K, and 16.72 and 16.78 Å for AKsub at 300 and 330 K, respectively. Their standard deviations were all less than 0.1 Å. This indicates the stability of the overall dimensions of the proteins during the simulations.
During the simulations the integrity of the zinc and magnesium binding sites was well maintained. The tetrahedral geometry in the zinc site was preserved. In AKsub the bidentate coordination of Asp 153 to the zinc ion was also preserved. The octahedral coordination around the magnesium ion was stable indicating that treating the interaction between the magnesium ion and its coordinating atoms as a tight electrostatic interaction was suitable.
To compare dynamics of the simulations with those in the crystal structures, the RMSF value of each residue during the last 1 ns was calculated and compared with crystallographic B factors. The calculated residual RMSF values correlate well with the experimental B factors. The linear correlation coefficients are 0.42, 0.56, 0.53, and 0.29 for AKste at 300 and 330 K and AKsub at 300 and 330 K, respectively. This shows that the simulations were fairly accurate in describing the dynamic feature of the systems.
Overall Flexibility of the Systems during the Simulations-The average RMSF value during MD simulation is considered as a measure of the overall flexibility of the system. The total and backbone RMSF values were calculated during the last 1 ns of the simulations, and the average values are listed in Table I. The total RMSF values were always higher than those for the backbone atoms in each simulation. This result confirms the less flexible nature of the backbone atoms. For both AKste and AKsub, the simulations at 330 K displayed higher average RMSF values than those at 300 K indicating increased flexibility of both systems at the higher simulation temperature. AKste and AKsub had similar average RMSF values at 300 K.
AKsub showed higher average RMSF than AKste at 330 K, but the difference was too small to conclude a higher flexibility of AKsub than AKste at 330 K.
Ion Pairs in the Simulations-From the structural comparison with mesophilic AKsub, it was suggested that electrostatic interactions of ion pairs in thermophilic AKste may contribute to its increased stability at high temperatures (21). Among many ion pairs found in AKste with 4-Å distance cutoff, only four ion pairs were identified as critical ion pairs for thermostability based on the distance between two residues in the sequence and the comparison with AKsub. They are Lys 19 -Glu 202 , Arg 116 -Glu 198 , Arg 131 -Glu 156 , and Lys 180 -Asp 114 (Fig.  2). They are all bridging the distant regions of a polypeptide (more than 10 residues), and their corresponding residues in AKsub do not form any ion pairs because at least one of the two residues is substituted. It is reasonable to choose ion pairs formed by residues from distant regions of a polypeptide because they are most likely not maintained when the protein is unfolded and thus contribute to stability by increasing free energy difference between the folded and unfolded states.
The four critical ion pairs were monitored during the last 1 ns of the simulation at 330K. The 330K simulation was selected because ion pairs are known to contribute to protein stability mostly at elevated temperatures due to reduced desolvation penalty for formation of ion pairs at high temperatures. It is also expected that the difference between AKste and AKsub would become more obvious at 330 K than at 300 K because AKsub starts to be denatured and inactivated in between 300 and 330 K but AKste remains stable at 330 K (18,19,21).
The distances of the four critical ion pairs over time are plotted in Fig. 3. The most striking finding is that a short distance for the ion pair Lys 180 -Asp 114 was not maintained during the simulation. There were moments when the two residues came within a relatively short distance (less than 6 Å) during the simulation, but they fell apart eventually. This suggests that the ion pair Lys 180 -Asp 114 may not be kept within a short distance in solution and thus not contribute to thermal stability of thermophilic AKste at high temperatures. Among the other three ion pairs maintained within short distances, the ion pair Arg 116 -Glu 198 was maintained within 4 Å without any transient separation, whereas the two ion pairs Lys 19 -Glu 202 and Arg 131 -Glu 156 spent considerable times away from the short distances seen in the crystal structure, although the separations between the pairing residues were recovered during the simulation. This results in the higher average distances for the Lys 19 -Glu 202 and Arg 131 -Glu 156 ion pairs than that of the ion pair Arg 116 -Glu 198 (Table II) and suggests that the Arg 116 -Glu 198 may be more important for stability than the other two ion pairs.
Simulations of AKsub Mutants-It is most likely that if the four ion pairs were engineered into AKsub, they would show similar behaviors to those in AKste considering the high structural similarity between AKste and AKsub, but this is not guaranteed. To test this proposition, an AKsub mutant (AK-subm1) having the four critical ion pairs was generated in silico, and its MD simulation was performed. The simulation was carried out as the same manner as the AKste and AKsub simulations and showed similar degrees of stability (data not shown). As expected, the four engineered ion pairs in AKsubm1  at 330K simulation displayed similarities to those in AKste (Table II). The two residues in the ion pair Lys 180 -Asp 114 were clearly separated from each other during the simulation, whereas the other three ion pairs were maintained within relatively short distances. This suggests that the ion pair Lys 180 -Asp 114 may not contribute to the stability of the AKsub mutant, and thus engineering only the other three ion pairs would be the more efficient way to increase thermostability of AKsub by introducing ion pairs. The three ion pairs would probably exhibit the same pattern when engineered in AKsub whether the ion pair Lys 180 -Asp 114 exists with them or not. It is, however, possible that the Lys 180 -Asp 114 pair and the other three ion pairs behave differently in the absence of each other. To exclude this possibility two more AKsub mutants were designed. The second AKsub mutant (AKsubm2) has only the three ion pairs that were kept in short distances. The third construct (AKsubm3) is a mutant with only Lys 180 -Asp 114 . Their simulations were performed and validated as for the others (data not shown). As shown in Table II, the maintenance of the short distance in the three ion pairs (Lys 19 -Glu 202 , Arg 116 -Glu 198 , and Arg 131 -Glu 156 ) and the dissociation between Lys 180 and Asp 114 were independent. The three ion pairs were kept within a short distance when the Lys 180 -Asp 114 was not with them, and the Lys 180 and Asp 114 were also separated without the other three ion pairs.
In an attempt to further dissect which is more important in stability among the three ion pairs maintained within relatively short distances, the simulations of three more AKsub mutants (AKsubm4, AKsubm5, and AKsubm6) that have only one of the three ion pairs (Lys 19 -Glu 202 , Arg 116 -Glu 198 , and Arg 131 -Glu 156 , respectively) were carried out and validated (data not shown). As listed in Table II, each ion pair was maintained within short distances in each mutant. It also seemed that the ion pair Arg 116 -Glu 19 8 was more closely maintained than the Lys 19 -Glu 202 and Arg 131 -Glu 156 ion pairs not only in the mutants having the multiple ion pairs but also in those having only one ion pair, although it was not as dramatic as in the case of the ion pair Lys 180 -Asp 114 .
With the analysis of the ion pairs in all the mutant simulations, it is tempting to draw simple conclusions about direct influences between the four ion pairs. For example, it seems that the absence of the Lys 180 -Asp 114 may tighten the other three ion pairs when comparing AKsubm1 and AKsubm2. However, in another simulation of AKsubm2 for verifying this, the average distances were 3.35, 4.17, and 5.98 Å for Lys 19 -Glu 202 , Arg 116 -Glu 198 and Arg 131 -Glu 156 , respectively. It is not clear whether the Lys 180 -Asp 114 is affecting the other three ion pairs or the differences are within the error range. Thus, it would be more appropriate to interpret the results of the MD simulations collectively rather than individually. Our hypothesis from the MD simulations is that Lys 180 and Asp 114 are not kept close to each other, and thus they do not increase thermostability by favorable electrostatic interaction, whereas the ion pair Arg 116 -Glu 198 is the most tightly maintained among the other three ion pairs and is the most important in protein stability.
In Vitro Measurement of Thermal Stability of AKsub Mutants-To verify the hypothesis generated from the simulations, in vitro experiments were performed. The six AKsub mutants studied in the simulations were produced and their T m s were measured by differential scanning calorimetry. In Table III, T m s of the six AKsub mutants are listed with differences from that of the wild type AKsub. Both AKsubm1 and AKsubm2 showed ϳ5°C increases of their T m s from AKsub, and the AKsubm3 had essentially the same T m as AKsub. This clearly indicates that having the three ion pairs Lys 19 -Glu 202 , Arg 116 -Glu 198 , and Arg 131 -Glu 156 is contributing to thermal stability, but the ion pair Lys 180 -Asp 114 is not useful for engineering stability in agreement with the result from the MD simulations.
The AKsubm2 appeared to be more stable than AKsubm1,  whereas the AKsubm3 showed almost the same T m as the wild type AKsub, suggesting a negative effect of the Lys 180 -Asp 114 pair on thermal stability when engineered together with the other three ion pairs. However, it is not clear that the reduced stability of AKsubm1 results from direct influence of the Lys 180 -Asp 114 pair on electrostatic interactions of the other three ion pairs because the difference in T m between AKsubm1 and AKsubm2 is not large enough to draw a definitive conclusion, and it is also possible that the destabilization may be caused by a more complicated indirect mechanism. Among the three AKsub mutants with one of the three ion pairs closely maintained in the MD simulations AKsubm5, having the ion pair Arg 116 -Glu 198 , showed the highest T m , and this is consistent with the hypothesis from the MD simulations that the ion pair Arg 116 -Glu 198 is the most tightly maintained in the three ion pairs and is the most important in stability. Surprisingly, AKsubm6 containing the ion pair Arg 131 -Glu 156 did not display an increased T m compared with the wild type AKsub. Although the inaccuracy of the MD simulations cannot be excluded as a source of the disagreement, it is possible that the ion pair is actually maintained but does not contribute to stability. In fact, it was found in the study of chimeras from AKste and AKsub that in certain areas of the adenylate kinases, residue substitutions from AKsub to AKste did not increase stability. 2 Rather, some of the substitutions decreased stability probably by disrupting complementary changes in AKsub because AKsub may have a different mechanism from AKste for local stabilization in the areas, and exchanging only one partner of cooperative interaction may result in destabilization.

DISCUSSION
Engineering ion pairs in an attempt to increase protein stability is not a simple task. Having positively charged and negatively charged residues in spatial proximity does not always guarantee tight ion pairs. Because proteins have finely tuned three-dimensional structures, naïve mutations for introducing ion pairs may rather disrupt already optimized local or global structures resulting in destabilization or inactivation. Studying ion pairs in thermophilic proteins can be useful either directly by introducing their ion pairs to mesophilic homologues or as general lessons for engineering thermal stability. If a structure of a thermophilic homologue of a target mesophilic protein is available and has ion pairs that the target does not have, one can replace corresponding residues in the target with the residues of the ion pairs in the thermophilic homologue. By doing this, it becomes more easy to decide which residues to mutate. It is also expected that the engineered ion pairs would act as in the thermophilic protein because of high similarity between a target and its thermophilic homologue.
It may seem easy to identify favorable ion pairs in thermo-philic proteins. However, choosing ion pairs in thermophilic proteins with a simple distance cutoff from their static structures may be misleading. Because most thermophilic protein structures have been solved by x-ray crystallography, probably due to their easy crystallization and because ion pairs are composed of hydrophilic charged residues, the majority of ion pairs are found on the surface of protein structures, which makes the residues susceptible to contacts with neighboring molecules in the crystal lattice. Therefore it becomes unclear whether the formation of ion pairs is natural or affected by crystal packing. In fact, three of the four critical ion pairs in AKste (Arg 116 -Glu 198 , Arg 131 -Glu 156 , and Lys 180 -Asp 114 ) have at least one residue whose side chain is contacting neighboring molecules in the crystal lattice within 4 Å. The flexible nature of the charged residues comprising ion pairs also makes the identification difficult. High solvent exposure and long side chains of charged residues often result in high flexibility and consequently ill-defined electron density with high B factors. However, it is not desirable to exclude residues having high B factors in search for true ion pairs because two residues can be flexible and form a tight ion pair at the same time by having concerted motions. MD simulation can provide insights to help solve these problems. Residues under the influence of crystal contacts can be relaxed during minimization and equilibration and likely mimic their behavior in solution during simulations. The flexibility issue is also overcome because it is possible to monitor motions of all the atoms all the time. It is therefore more straightforward to discern true ion pairs by selecting ones maintained during simulations rather than using a simple distance cutoff with average static structures. The results of MD simulations can also be beneficial to other approaches using the simple electrostatic models, because MD simulations can provide explicit information about geometry of charged residues (31).
Although it becomes easier to identify true ion pairs in thermophilic homologues based on MD simulations, it is not always guaranteed that the ion pairs identified in the thermophilic proteins remain tight when engineered in the target proteins. It is also possible that the structures of thermophilic homologues may not be available, and residues to mutate in target proteins should be decided by less reliable methods. In both cases, MD simulations can be helpful to predict the consequence of the mutations. If a structure of a target is known and residues are selected for mutations to build ion pairs, in silico mutations of the target can produce a series of engineered protein structures with different combinations of ion pairs. MD simulations of the mutant structures can be performed and analyzed for dynamic motions of the engineered ion pairs. Consequently the mutants in which the engineered ion pairs remained tightly in MD simulations can be selected for in vitro mutagenesis.
A major limitation of this strategy is that MD simulations only provide information about ion pair geometries, not about the overall stability of proteins, because it is not possible to accurately calculate protein stability in terms of thermodynamic quantities from results of MD simulations. It is in fact proposed that ion pairs can be tightly maintained without contributing to overall stability (32). It has been suggested that stability and flexibility of proteins are closely related (33). However, no definite correlations have been found between stability and overall flexibility in the MD simulations of thermophilic and mesophilic homologues (34 -36). It is possible that the time scales of the simulations are too short or the force fields are not sufficiently accurate to describe the flexibilities of the proteins. Moreover, recent experimental studies suggest 2 E. Bae and G. N. Phillips, unpublished data. that flexibility adjustment for temperature adaptation may vary locally (37,38). Thus, homologue proteins may not share similar global flexibility at their optimal temperatures and relations between stability and overall flexibility may not be expected. In our case, no clear correlation was observed between average RMSF and thermal stability for the pair of AKste and AKsub (Table I) or for the AKsub mutants (data not shown), and it was difficult to assess stability based on overall flexibility in MD simulations. Despite the limitation, it is clear that the strategy presented here can help identify critical ion pairs by explicitly monitoring motions of ion pairs or at least exclude false ion pairs, because a pair of separated residues can never contribute significantly to stability through electrostatic attraction. In our case, one of the four critical ion pairs selected from the static crystal structure of thermophilic AKste was suggested to be false based on its MD simulation. We were also able to propose that another ion pair is the most important in stability among the three ion pairs closely maintained in the simulation. Additional MD simulations of the AKsub mutants suggested similar behaviors of the ion pairs when engineered into the mesophile. From these results it was hypothesized that the Lys 180 -Asp 114 would not contribute to thermal stability whereas Arg 116 -Glu 198 would the most significantly in AKsub. This hypothesis was experimentally validated by producing the AKsub mutants in vitro and measuring their thermal stability by differential scanning calorimetry. This suggests that MD simulations can provide valuable information for identifying and engineering ion pairs and potentially replace substantial part of expensive and time consuming in vitro experiment.