Effect of Charge Distribution in a Modified tRNA Substrate on Pre-Reaction Protein-tRNA Complex Geometry.

An important aspect of molecular mechanics simulations of a protein structure and ligand binding often involves the generation of reliable force fields for nonstandard residues and ligands. We consider the aminoacyl-tRNA synthetase (AaRS) system that involves nucleic acid and amino acid derivatives, obtaining force field atomic charges using the restrained electrostatic potential (RESP) approach. These charges are shown to predict observed properties of the post-transfer editing reaction in this system, in contrast to simulations performed using approximate charges conceived based upon standard charges for related systems present in force field databases. In particular, the simulations predicted key properties induced by mutation. The approach taken for generating the RESP charges retains established charges for known fragments, defining new charges only for the novel chemical features present in the modified residues. This approach is of general relevance for the design of force fields for pharmacological applications, and indeed the AaRS target system is itself relevant to antibiotics development.

Introduction

Advancement of computational and natural sciences could provide answers to the most significant and complicated questions, usually related to the study of mechanisms or estimation of the probability of the process. Complex systems containing nucleic acids require individual approaches in order to explain and interpret some experimental data. Evidently, the development of a force field should find the middle ground between its simplicity and accuracy to speed up the calculation without loss of the fidelity. The common challenge for such balancing is a treatment of nonbonded interactions, which are one of the most sensitive elements of molecular mechanics calculations due to the growing data array to process on each step of molecular dynamics (MD). The Amber standard force field is one of the most used for MD, especially for systems, which include small molecules and nucleic acid components. Amber was initially developed for a range of parameters between nucleic acid bases interactions. For standalone nonstandard residues like ligands, there are plenty of software and services to calculate missing parameters. To simulate a novel molecule, it is just a matter of assigning the desired atom types and reliable charges. However, due to the optimization process, which depends on which method is used for the generation of charges, it is necessary to apply those that maintain the force field integrity.

The study aims to determine a convenient and advanced method of appropriate force field parameterization of an aminoacyl-tRNA fragment (“charged tRNA”) to automate and facilitate MD simulations. In a fundamental aspect, such problems are associated with the study of aminoacyl-tRNA at several levels of decoding the genetic information, from the synthesis of an aminoacyl-tRNA molecule to the biosynthesis of the polypeptide chain on the ribosome. It is of interest to study erroneously synthesized aminoacyl-tRNAs, as well as those aminoacyl-tRNA, which are formed with non-proteinogenic natural amino acids (norvaline, norleucine, homocysteine, D-amino acids, etc.). On the other side, the model of an aminoacyl-tRNA synthetase (AaRS) complex is a peculiar playground for comparing different approaches and finding a compromise between performance and accuracy. The problem of derivation of new force field parameters also is needed for RNA derivatives containing a wide range of rare nucleic acids bases. In tRNA, a number of nonstandard bases (for instance, modified by methylation or the inclusion in the sixth position of the adenosine derivative of threonine, etc.) take part in the decoding accuracy of genetic information and several other regulatory functions in the cell. In the case of DNA, it is epigenetics, modification of DNA in various pathological processes, etc.

Previously, several laboratories published their studies, with incompatible partial charge distribution for the substrate and therefore incorrect substrate motion and binding modes.[1,2] Here, we propose a new model of aminoacyl-tRNA fragment parameterization, based on electronic structure calculations and RESP (restrained electrostatic potential) charge fitting,[3] followed by the MD simulations in an explicit solvent compatible with the Amber99sb force field.[4] The validation of the force field parameters and the atomic charges was performed against experimental data. Finally, the results from our MD simulations met and explained an existing biochemical data. Our model comprises both protein and nucleic components suitable for simulation in a composite all-atom force field.

In this paper, we combine ab initio quantum mechanical (QM) methods and MD simulations to derive and validate force field parameters for several substrates of aminoacyl-tRNA synthetases (AaRS) compatible with the Amber99 force field family. The molecule of interest is an aminoacyl molecule formed with the amino acid and 3′-terminal adenosine of tRNA. Our goal is to understand the putative mechanism of the editing of the aminoacyl molecule formed with the incorrect amino acid and 3′-terminal adenosine of tRNA. We also present some details of the modeling process.

Results and Discussion

Common Techniques of Parameterization

The electrostatic potential (ESP) of a molecule is a description of charge distribution and is regularly used in medicinal chemistry, modeling, and computational chemistry.[5,6] One of the ways to circumvent some force field limitations is a RESP approach used to assign partial charges. This is based on an accurate quantum mechanically calculated MEP (molecular electrostatic potential) minima using an atom-centered point charge model of Amber force fields. To compute missing charges of unparameterized compounds, the following two procedures were applied to make topologies suitable for the Amber99 force field family:

(1) Bond-charge correction (BCC) or the AM1-BCC model,[7] which uses an improved semi-empirical method for the generation of molecule-specific point charges.

(2) Computing RESP charges[8] to obtain point charges for atoms in molecules, which are compatible with Amber force fields. The procedure requires calculation of the molecular ESP at the quantum-chemical level for the target molecule.

The first method is suitable for standalone molecules with an integer charge, like aminoacyl-adenylates,[9,10] or standalone aminoacyl molecules,[11] and it is implemented in the Antechamber module of the AmberTools package.[4] However, this algorithm is not suitable for chain monomers and their modification (acetylation, methylation, non-canonic nucleotides, covalently bonded cofactors). The ESP for such molecules could be derived from quantum mechanical calculations through GAMESS and Gaussian packages. However, previously it was shown that electrostatic potential on the molecular surface is conformationally dependent and is formed with unphysically high charges on atoms inside the molecule.[12] To avoid this issue and to demonstrate the importance of correct charge mapping, we processed our nonstandard residues using an RESP charge model. GAMESS and Gaussian software can produce several orientations for the lowest-energy native-like conformations. As a result, different sets of RESP/ESP charges can be calculated for the same geometry depending on the orientation in space. To solve the problem, the rigid body re-orientation algorithm, which is implemented in the R.E.D. III server, was applied to each minimized structure right before the MEP calculation to get reproducible RESP charges. This step is equally important to the search of the ligand-binding conformation or optimization and equilibration of the system.

Parameterization of the Aminoacyl-tRNA Fragment

Even though oligonucleotides with terminal 5′-phosphates can be simulated with either force fields, one of the most notable features of Amber force fields is the designation of the terminal nucleotides and polarizability of atoms, which is strongly affecting the stability of nucleic acids.

Aminoacyl-tRNAs are specific substrates of AaRSs, formed with an ester linkage between the carboxy group of an amino acid and the 3′OH group of a nucleotide. According to the General Amber Force Field parameters, the internal nucleotide in the chain consists of a phosphate group, certain base, and a sugar. One of two hydroxyl oxygen atoms (O3′) is engaged in forming a phosphodiester bond with the next monomer; thus, it is deprotonated. Such state is characterized by a negative integer charge of −1′. To provide structural stability along with the MD simulation, the phosphate group of a 5′-terminal nucleotide is replaced with hydrogen (5HT).[13] Due to the well-refined charge distribution model it affects only charges on the nearest joined atoms, however, increasing the net charge up to −0.3081′. In turn, the O3′ atom of the 3′-terminal nucleotide becomes protonated (Figure ). This increases the net charge of the 3′-terminal nucleotide up to 0.6919 and compensates the −0.3081, as the sum of these partial charges is equal to −1. Thus, the problem of terminal residues’ modification arose from its non-integer charges.

Figure 1

All chained and nonterminal RNA nucleotides in the Amber force field could be represented as a set of fragments: a phosphate group, nucleoside, and its O3′ oxygen (A). The sum of these charges equals −1. However, elimination of a phosphate group from the 5′-terminal nucleotide leads to the formation of a fractional net charge on this nucleotide (B), which is compensated with automatic 3′-terminal nucleotide modification (H3T hydrogen capping of O3′). The color scheme of the atom types corresponds to the Supplementary Information (Table S1).

To generate a set of satisfactory models we computed charges for aminoacyl-tRNA fragments using the RESP method. The quality of the prediction was estimated with a relative root mean square (RRMS) fit, a way to compare the ESP models derived from single molecule structures and under the charge constraints. Finally, R.E.D. algorithm integrates charges into a separate force field library, ready for use in molecular dynamics simulations.[14]

We prepared two residues, namely, a nucleotide and amino acid, with recalclulated partial charges. As a result, we obtained the unified charge model for the amino acid-bound nucleotide with significant changes related only to the atoms involved in aminoacyl ester formation. Two atoms of the nucleotide, C2′ and O2′, have undergone significant changes, as they are involved in aminoacyl ester bond formation. At the same time, the general model of charge distribution and the total charge of the chained base retained their original form. With regards to the amino acid residue, which becomes the terminal residue in the chain, its initial charge of +1 turned into +1.3081 (Figure ). We also attempted to create a model where the last nucleotide from the 3′-end tRNA covalently bound to the amino acid was represented as a single residue with the output charge of −.6919. However, this approach resulted in a high relative root mean square (RRMS) value (greater than 0.56) and unstable behavior of the substrate during MD simulation.

Figure 2

Mechanism of charge fitting using R.E.D.III algorithms. Modification of an initial structure (3′-terminal adenosine or RA3) with default charge distribution to obtain a novel terminal nucleic acid structure (RA-M) (A). It contains a capping H3T hydrogen atom, like terminal nucleic residues, and lacks another capping hydrogen atom on O2′. In doing so, it has a net charge of −1, inherent to nonterminal nucleic residues. Net charge values are linked to each monomer and are colored in red. Combination of charge assignment to the group of atoms (M and Z) in the input structures and subsequent elimination of the group from the output structures applied to nucleic acid and amino acid residues (B). The final charge distribution model of a schematic aminoacyl-tRNA molecule (C).

In fact, there was no significant difference between single-point energy calculation of multiconformation (or orientations offered with RBRA) for the amino acid part of the aminoacyl. Compact and nonpolar substituents, like 1-methylethyl ether and OMe, gave the smallest error for amino acid-based charge fitting. In turn, acetate and glycine-substituted nucleotides showed the lowest RRMS values and best results in MD. Some of the best (the most dynamically stable) charge and atom typing parameters of aminoacyl components are represented in the Supporting Information (Table S1).

Example of Application in a Fundamental Study

Biologically, any AaRS is responsible for the two-step aminoacylation reaction in the aminoacylation site and several editing mechanisms (Figure ). The first step initiates activation of amino acid by ATP to form aminoacyl adenylate. Chemically, this process is represented with inorganic pyrophosphate cleavage and further phosphoester bonding between AMP’s phosphate and carbonyl of the cognate amino acid. The fidelity of the substrate recognition process is provided with a pre-transfer editing mechanism, which prohibits the reaction based on the size and shape of an amino acid. The second step includes a subsequent transfer of an amino acid moiety to the 3′-CCA end of the tRNA molecule, so called aminoacylation, forming aminoacyl-tRNA. This transesterification process results in bond rearrangement and connection of a OH group on the ribose of an 3′-terminal adenine base to a carbonyl carbon of the amino acid molecule. In some individual cases, the more precise post-transfer editing mechanism (for example CP1 domain of Leucyl-tRNA synthetase) can bind and hydrolyze wrong aminoacyl-tRNA to avoid integration of similar, but noncognate amino acids, into the newly synthesized protein.

Figure 3

General representation of LeuRS functional sites of the catalytic core and the CP1 domain. The tRNA molecule is shown in pre- (orange) and post-transfer orientation of the 3′-termini.

The crystal structures of T. thermophilus,P.horikoshii LeuRSs were taken from the Protein Data Bank (PDB IDs: 1OBH and 1WKB).[15] To determine the mechanism of amino acid selectivity at the CP1 domain, a docking procedure and MD simulation were performed. In the case of aminoacyl-tRNA fragment modeling with the original Gromacs GMX force field, the main challenge has been met in a correct ether bond parameterization between nuclei and amino acid moieties, which resulted in excessive flexibility of the fragment in the binding site of the CP1 domain (Connecting Peptide 1/editing domain). Backbone conformation of nucleic acids has six torsion angles; thus, a simulation of nucleic acids with increased flexibility is far more complex than those of proteins. Along with a backbone angle, these specific features are completely important during the protein–nucleic recognition. The same is true for aminoacyl recognition and subsequent water molecule approximation to the bond between an amino acid and nucleotide preceding the hydrolysis process. Initial configuration of the entire system demands a very precise treatment for the long-range electrostatic interactions, while the water motions around the ligand, which precedes a nucleophilic attack, does not take long. Therefore, the system should be carefully designed and well equilibrated, as the internal strain of the tRNA contour could disrupt any interactions of amionacyl-tRNA motion with the CP1 domain of LeuRS.

LeuRS is a characteristic enzyme representative of the class I of AaRSs, aminoacylating the 2′OH atom of ribose and possessing the editing activity. Normally, leucyl-tRNA should not be hydrolyzed with the CP1 domain. However, the mutation of the T252A residue increases the rate of such event, allowing to bind leucyl- and isoleucyl-tRNA with its subsequent division in primary components.[16,17]

The identification of this known mechanism can take place by a two-stage reaction. On the first stage, the hydrolyzed molecule should take the right geometry, two pairs of H-bonds (Thr247 with a carbonyl oxygen/Asp347 with an amino group of the amino acid) stabilize and activate the plane of a carbonyl group. Then, a water molecule (W1) should attack the activated carbonyl carbon, being activated with another water molecule (W2), an assistant molecule, preferably forming H-bonds with surrounding amino acids of the binding site. To construct and simulate leucyl-tRNA in the editing state, bound to the wild-type and T252A mutant CP1 domains, we used an R.E.D. III charge-fitting procedure and an alternative method of simple topology combination described in Hagiwara et al.’s article.[1] The last one is supposedly simple combination of atom types, bond orders, and charge values into the same section of the residue topology database.

When the R.E.D. III algorithm was applied, the conformation and water accessibility of the ester bond in the leucyl-tRNA fragment significantly depended on the residue in the 252 position of the CP1 domain. In the WT protein, the molecule of leucyl-tRNA is weakly interacting with Asp347 and rarely with Thr247 (Figure A,C), and a side chain of leucine is exposed to the binding site. At the same time, T252A mutation creates an additional space, which is sufficient for the location of leucine and an appropriate orientation of the ester’s bond plane (Figure B). H-bond interactions with Thr247 and Thr248 increase the probability of a nucleophilic attack on the carbonyl carbon of the ligand due to stabilization of the geometry and pulling of the electron cloud density. Simultaneously, Asp347 forms a strong interaction with the amino group of leucyl-tRNA with the decrease in the number of degrees of freedom. The results of MD analysis, namely, graphs of interaction energies and RMSD of the ligand and the CP1 domain are shown in Figure S1.

Figure 4

Most stable examples from MDs of CP1 with Leu-tRNALeu, prepared using charge fitting (A, B). The MD demonstrated the dependence between a mutation in the 252 position, an interaction map and conformational changes (C). A charge fitting model of Leu-tRNALeu preserved the initial angle of amino acid turn toward the ribose group, forming H-bonds with T247 and D347. These cause the increase of water accessibility and stable distances from ligand atoms to the binding site residues (D). The total number of water molecules in appropriate location toward the ester plane during MD, W1; the frequency of water pairs (in attacking and assisting modes) detected during MD, W1 + W2; frequency of the occurrence of initial coordinates of the reaction (a pair of water molecules and the correct geometry and interaction map of the protein–ligand complex), W1 + W2 + T247.

However, the method of simple concatenation of topologies turned out to be somehow unrepresentative and inappropriate for the mechanism study, particularly for the absolutely similar stability of either substrates and inability to interpret our in-house and already published biochemical data.[18,19] Incorrect charge distribution in the model, derived from the simple topology combination, causes an increased rotation of the amino acid radical independently of the residue in the 252 position and loss of interaction with T247 (Figure S2). Application of the combined topology of N-Leucine with a non-terminal adenosine did not demonstrate any significant differences in the binding mode of Leu-tRNALeu in either wild-type or mutant proteins. In general, Leu-tRNALeu conformation residues are the same during the MD simulation regardless of mutation, but the number of water molecules near the carbonyl carbon of leucine increased in the case of T252A simulation. The interaction of Leu-tRNALeu with amino acids of the binding site also wasn’t affected with mutation; in both cases, Thr248 and Asp347 formed strong H-bonds with the 3′OH group of ribose and the amino group of leucine, respectively. At the same time, the Thr247 residue, which was proven to be critical for the hydrolysis process, did not form contacts with the carbonyl oxygen of leucine. Thus, even in the presence of water molecules the geometry of the aminoacyl poorly suites the reaction requirements (Figure S2 and Video S1).

Each result obtained with RESP charge assignment was tested with several AaRS systems. We applied the method to compare both norvalyl-tRNA and isoleucyl-tRNA in the editing site of the native and mutated CP1 domains of LeuRS (Figures S3 and S4). To generate transferable RESP point charges, some averaging of the model was necessary, because this approximation strongly affects conformational variation of studied amino acid moieties. Based on the dual norvalyl- and isoleucyl-tRNA model, we managed to predict a common binding mode in LeuRS from bacteria and archaea, however, with different directions of the water attack, as it was reported previously.[20,21]

Conclusions

One of the most important findings of these in silico experiments on aminoacyl-tRNAs using the RESP charge deviation method is the development of accurate and adjustable topologies, which are not possible to create with different third-party software. The specificity of the Amber force field is rendered challenging to reproduce a correct charge distribution for such molecules like aminoacyls, which comprise both amino acid and nucleic moieties. Since it is a chained terminal residue and its net charge is not an integer, it causes difficulties to parameterize the molecule according to force field rules with any software intended for calculation of single molecules. The great advantage of this study lies in the fact that all results were consistent with in-house and already published in vitro tendencies.

A practical implementation of the protocol and its advantages over other described methods were proven in the study devoted to the decrease of the post-transfer editing reaction efficiency in mutated LeuRS proteins. In the earlier studies, the importance of the canonical T252 position for the editing substrate recognition was shown.[16,17] The study of this mutation and computational reproduction of the evident difference between the wild type protein and T252A mutant should become a test system for the parameterization protocol. In an attempt to reproduce the result of in vitro findings, two charge distribution approaches were compared and assessed via an MD simulation method. Based on the described protocols applied to the object (a combination of defined Gromacs parameters from *.rtp (residue toplogy associated file format) database and derivation of the charge model with a charge fitting algorithm of R.E.D.III), we prepared two topologies for each of leucyl-, norvalyl-, and isoleucyl-tRNA substrates. The starting conformation of each substrate after the docking procedure was not changed during the topology preparation state. At the same time, the mutations we studied (T252A) could affect the interaction map of amino acid residues but not the nucleic moiety of aminoacyl (Figure S5).

The impact of each mutation was assessed by the calculation of the nucleophilic attack probability rate. Apart from simple visual motions stability of the H-bond network, changes of dihedral values and fluctuation of interaction energy indicators and of course biochemical data are highly significant for the subsequent study and analysis (Supporting Information).Thus, in addition to the simulation of the pre-reaction conditions, we determined a direct dependence between the approach of topology generation and the correlation with experimental data.

The most promising and expensive quantum calculation can demonstrate unrealistic behaviour of the the incorrect pre-reaction geometry. Based on the obtained data, we are going to predict a more detailed mechanism of the hydrolysis reaction. Several of the most stable pre-reaction states from MD simulation of LeuRS from T.thermophilus could be treated as the initial coordinates for DFT QM calculation. It is reasonable to reduce the entire complex to the binding site forming amino acids (28 residues), a norvalyl-tRNA fragment, and four interacting water molecules to simulate the reaction.

Computational Methods

The Amber force field allows multiple nucleotide modifications, if it does not affect the net charge conditions, which strictly depends on the position of the element in a chain. We used the R.E.D.III server because it allows calculation of a non-integer charge for the target molecule applying a special constrain algorithm. This server was designed to generate non-polarizable RESP and ESP charges for new molecules/molecular fragments (with appropriate force field library format).

Initial Geometry of Aminoacyl

We used an optimized geometry of the entire complex to extract the coordinates of norvalyl-tRNA (PDB-ID: 2BTE) and docking poses for isoleucyl- and leucyl-tRNA (Figure ).[20,21] The same approach was used for the study of alanyl-tRNA binding in the prolyl-tRNA synthetase (PDB-ID: 2J3M).[22,23] All details for ligand–protein and nucleic acid–protein docking protocols with Gold CCDC[24,25] are already described in previous studies.[22] Selected conformations of an aminoacyl-tRNA fragment were prepared in a PDB format. Target files were processed with an Ante_R.E.D. 2.0 program, interfaced by an RESP server,[26] in the automatic mode. Thus, proper input structure geometry was converted to p2n type files to rigorously define the important elements.

Figure 5

Interaction map and important dihedral angle derived from 2BTE crystallographic data. An unnatural intermediate is formed with amide instead of an ester bond (labeled with N). Schematic representation of aminoacyl-tRNA fragment composition from an adenosine nucleotide and an amino acid N-terminal residue. Norvaline (NVA) and isoleucine (ILE) residues backbones with a protonated amino group, as it is required by the Amber force field.

Partial Atomic Charge Derivation

To deal with the task, we turned the last 3′-terminal nucleotide into a penultimate monomer of the chain, while an amino acid took its place. The only change in the charge on this nucleotide ought to apply to C2′–O2′ atoms, as they formed a bond with non-native nucleic acid connection atom types (the carboxy group of an amino acid). At the same time, the integer charge of the amino acid had to be constrained somehow to comply with the net charge of a 3′-terminal monomer or, maybe in some way, compensate the 5′-terminal net charge (Figure ). Charge distributions for each residue are shown in Tables S1 and S2of the Supporting Information. A color scheme of tRNA atoms was based on atoms appearing in the structure and its charge values. A group of atoms with similar charges and common for all models is colored in cyan. Some atoms, which possess slightly different charge values but do not significantly affect the entire charge distribution, are shown in yellow color. The atoms, which are absent in just one structure, are marked with green color. Atoms, which are critical for the structural integrity of charge changes, are shown with red color.

Accordingly, we generated two sets of molecules, containing either nucleic or amino acid as a core and substitutes of different natures. However, both sets shared through a common triad of atoms to mimic the ether/ester bond between the sugar and amino acid.

One set of molecules was built to calculate the charges on C2′ and O2′ atoms of the nucleotide, which are changing when the phosphodiester type bond is replaced with an ester bond. It was based on the adenine structure bound to methyl, and glycine and acetate groups bound to O2′ of ribose (RA-M, where RA is RNA adenine and M is a substituted group). Another set was generated from the amino acid (+NH3CHXCO, where X is the amino acid’s side chain) connected through the carboxyl oxygen to the group like OH, OMe, and ether bonded methylethyl and furan (Z). Both structures, +NH3CHXCO-Z and RA-M, contained the formal substituted group to be removed together with its part of the net charge. This part of the charge-fitting step was performed with the implementation of R.E.D.III-specific intramolecule (INTRA-MCC) and intermolecule (INTER-MCC) charge constraint options (Figure A).

INTRA-MCC (inside molecule) charge constraint of 0’was applied to the removal of the M-group of the RA-M residue (Figure A). The default charge distribution of the major part of the nucleotide remained intact, while significant changes are related only to the atoms involved in aminoacyl ester formation. The net charge of the nucleotide was still −1, but C2′–O2′ charges showed that it is a part of a larger molecule.

Similarly, the INTRA-MCC value of −0.3081 was applied to the Z group of +NH3CHXCO-Z, causing the appearance of a +1.3081 charge (Figure C). Simultaneously, we applied INTER-MCC to reproduce charge distribution, which is intrinsic to N-terminal residues, following the R.E.D.III documentation (chapter IV.2 of R.E.D. server: Examples & demonstration). Charge derivation for this fragment was carried out by setting an intermolecular charge constraint between the methyl group of methylammonium and the MeCO–NH group of atoms of the capped amino acid. Force field library building for this fragment involves removing all the atoms involved in these two constraints and adding new atom connectivity between the nitrogen atom of methylammonium and the alpha-carbon of the capped amino acid.

Geometry Optimization and Maintaining of the Force Field Library

We optimized all organic fragments (amino acid and nucleic acid) using the B3LYP/cc-pVTZ theory level, implicit solvent model (SCRF) theory level,[27] and the Gamess (generalized atomic and molecular electronic structure systems) software package (version 7.1.5) to obtain correct MEP. For the substituted residue fragment, several energy minima were selected, based on mimicking the docked conformation, suggesting a more realistic binding mode. A rigid body reorientation algorithm (RBRA) implemented in the R.E.D. server was applied for the geometry of each block containing the atoms C2′, O2′, and C as input depending on the model of interest. This step provides molecular orientation of the geometry and the reproducibility of the atomic charge values independently of the QM program or the initial structure. The Amber/Gromacs topologies for novel residues with all necessary parameters (bonds, angles, dihedrals, etc.) and charges were generated using the tLEaP module of AmberTools18.[28]

Molecular Dynamics Simulation

All MD simulations were run using the Gromacs ver 5.0.7 program with the all-atom type Amber99 force field.[29] Each system initially contained the solute protein–nucleic acid complex with 10 Ǻ of TIP3P water. Each system was minimized using 10,000 cycles of conjugate gradient minimization followed by 30 ps of molecular dynamics equilibration. To ensure that tRNA conformation is stable and 3′CCA (CCA 3′-terminal group of tRNA) will have a correct bond in the production run, MD simulations of 500 ns for LeuRS (from T. thermophiles) complexes with corresponding “discharged” tRNA molecules were carried out. The total number of atoms in simulated systems was in the range of 200 thousand for monomeric LeuRSs. A positional constraint was applied for the first 20 ns to a group of atoms from 3′-terminal adenosine (hydroxyls and planar ring atoms), and the editing site (SER227, ARG346, and LEU329) to fix a tRNA stem coordinates. After these relaxation steps, our complexes consisting of a protein and tRNALeu were applied with a harmonic function with a force constant of 500 kcal/mol per A2. To accurately switch from constrained MD to free MD, the force constant was reduced to 250, 125, 50, 25, 10, and 5 kcal/mol per A2 in six MD simulations each of 500 ps. After that, a free simulation of 20 ns MD was run at 333 and 310 K for thermophilic and nonthermophilic proteins, respectively, and the long-range interactions were evaluated using the PME (Particle–Mesh–Ewald) method. The temperature coupling mode using velocity rescaling was used together with the Parrinello–Rahman coupling algorithm. The Coulomb cutoff radius of 1.2 nm for the electrostatic and cutoff radius of 1.1 nm for Lennard–Jones interactions were applied. Analysis of the relative position of the substrate, all available water molecules and amino acids of the binding site, was performed with a purpose-written script in Python that integrates step-by-step analysis of all molecules and residues, surrounding an aminoacyl-tRNA fragment during MD simulation.[30] All results of the analyses were performed using Gromacs built-in tools and the last 10 ns were represented in graphs.