Large-scale simulations of the ribosome: a new landmark in computational biology

Within computational biology, all-atom simulation is the most computationally demanding field, in terms of compute load, communication speed, and memory load. Here, we report molecular dynamics simulation results for the ribosome, using 2.64 × 106 atoms, the largest all-atom biomolecular simulation published to date. The ribosome is the largest asymmetric biological structure solved to date to atomic resolution (2.8 Å). While simulations requiring long-range electrostatic forces have been previously restricted to much smaller systems, breakthroughs in electrostatic force calculation and dynamic load balancing have enabled molecular dynamics simulations of large biomolecular complexes. The LANL Q Machine played a key role in enabling such large simulations to be performed. The LANL Q machine displays approximately 85% parallel scaling efficiency for the ribosome system on 1024 CPUs. Using the targeted molecular dynamics algorithm, we have simulated the ratelimiting step in genetic decoding by the ribosome. The simulations use experimentally determined ribosome structures in different functional states as the initial and final conditions, making our simulations entirely and rigorously consistent with these experimental data. The simulations have identified candidate 23S rRNA nucleotides important for the accommodation of tRNA into the ribosome during protein synthesis.


Introduction
Computational biology is a relatively new field in comparison to computational physics. While systems biology, computational neuroscience and bioinformatics are each potential high performance computing applications of the future, all-atom biomolecular simulations (both quantum calculations and simulations using classical molecular dynamics potentials) are high performance computing biology applications of the present and near-future. Here we describe the first million-atom simulation in biology. By comparison, the first million particle simulations in materials science and cosmology were performed over a decade ago. 1; 2 Today, multibillion-particle simulations are performed in these fields. 3; 4 The calculation presented here uses ~ one million CPU hours on 768 CPUs, or about 10% of the 13.88 TFLOP/s LANL Q Machine, enabling us to simulate approximately 22 ns of physiological time of a large biomolecular complex. 5 Targeted molecular dynamics drives the system from one experimentally determined state to another experimentally determined state, while simultaneously accounting for thermal fluctuations. 6 The simulations allow us to interpolate structures that occur between the two states during conformational changes important for the function of the biomolecular complex. These interpolated structures are then used to predict which regions of the complex are key to the biological function. The predictions can then be tested experimentally by site-directed mutagenesis in vitro or by a variety of single-molecule fluorescence resonance energy transfer (FRET) techniques.
Although this simulation is the largest of its kind, it represents a small fraction of the simulation time required to observe these conformational changes spontaneously, as opposed to driving the system using targeted molecular dynamics. Rapid kinetic experiments show rates of spontaneous conformational changes for this system on the order of ~ 100 ms, or seven orders of magnitude greater than the physiological time simulated in our current simulations. Based on simulations of this system using the Q machine, we estimate that a 100 ms simulation would require ~ 1000 Peta-FLOP/s using the NAMD code, or ~ 100 PetaFLOP/s, using a more scalable version of GRO-MACS 7 or perhaps LAMMPS. 8 Even including a factor of 10 speed-up for an unforeseeable improvement in molecular dynamics algorithms, a 10 ms simulation of this system would still require 10 PetaFLOP/s, about 2 orders of magnitude greater than current capabilities. Nonetheless, the 22 ns simulation of the large biomolecular complex presented here is useful to experimentalists, who are already testing our predictions. The simulation serves as a proof-of-principle for simulations of large-scale biomolecular complexes.
The biomolecular system simulated here is called the ribosome. The ribosome is a giant molecular complex central to all living things and ancient in evolutionary terms. The job of the ribosome is to read genes and synthesize the proteins that the genes code for. The ribosome performs a complex 'look-up table' operation, reading sequences based on the 4-letter alphabet used by DNA and RNA, and converting these sequences into new sequences written in a 20-letter alphabet (i.e., the amino acid sequences which constitute proteins). The ribosome accomplishes this feat with a suite of adapter molecules called transfer RNAs (tRNAs), which were first predicted to exist by Francis Crick 9 , co-discoverer of the DNA double helix structure. The ribosome is a nano-scale computer that acts as the central information processing unit of the cell and is very much analogous to the central processing unit (CPU) of the cell. The ribosome is so fundamental to life that many parts of this molecular machine are identical in every organism ever sequenced, suggesting that these parts existed in their current form in the last common ancestor of life on earth. 10 These 100% evolutionarily conserved portions are called universal bases.
We have simulated the key step by which the ribosome decodes genes using the ASCI Q Machine and identified 18 universal bases important for decoding. 5 Until now, static snapshot structures of the ribosome have been available. 11 Limitations in time resolution and spatial resolution prevent experimental imaging of the ribosome in motion in atomic detail. Our simulations allow us to image the motion of transfer RNAs inside the ribosome occurring during decoding.

Materials and Methods
To measure the performance of the NAMD code on the LANL Q machine, we performed scaling studies on several systems with different numbers of atoms. The Q Machine currently has 2048 HP Alphaserver ES45 nodes, each with 4 EV6 1.25 GHz CPUs and an 8 MB cache. The nodes are connected by a Quadrics interconnect with ~2 µs latency and 300 MB/s bandwidth.
The systems simulated include: (1) the small subunit of the ribosome (N atoms ~ 1.07x10 6 ) 12 ; (2) the 70S ribosome (N atoms ~ 2.03x10 6 ) 13 (Fig. 1); (3) larger 70S ribosome system with larger messenger RNA (N atoms ~ 2.64x10 6 ) 5 . Ions were placed randomly in a box around the solute at concentrations of 0.1 M KCl and 7 mM MgCl 2 . Ions were equilibrated with the NAMD 14 molecular dynamics simulation code and AMBER 15 force field using a continuum water model with 5 Å radius for the ions to ensure electrostatic energy convergence. For the whole ribosome simulations, the ion-solute sys-tems were equilibrated for 10 ns. Electrostatic energy equilibrium was reached in approximately 1.5 ns. The ion-solute systems were embedded in a TIP3P water solvent box using the solvate routine. 16 The ion-solvent-solute systems were then minimized. The solute and ions were gradually heated via constant volume molecular dynamics and temperature coupling from temperature T =10 K to T=300 K over 200 ps, while keeping the solute fixed in place. The molecular dynamics simulations use particle mesh Ewald electrostatics with a grid spacing of ~ 1 Å, cutoff = 9 Å, multiple time steps, and NAMD. The system was then equilibrated with respect to volume running at constant pressure, P= 1 atm, for 200 ps using Langevin-Nose-Hoover pressure coupling. The solute was restrained with harmonic positional restraints at 200 kcal/mol Å 2 which were gradually lowered to 1 kcal/mol Å 2 over ~ 1.2 ns. The restraints were set to 0 with the exception of bases near the large subunit proteins, which could not be modeled in the case of 70S ribosome simulations. In total, the equilibration time was approximately 1.6 ns. Figures were made with VMD 17 and Molmol.
18 Figure 1. The ribosome complex after the initial equilibration stage of the simulation for the smaller of the two whole-ribosome systems. Blue, 16S rRNA; magenta/pink, 23S rRNA; purple, 5S rRNA; cyan, ribosomal proteins; grey, peptidyl tRNA; green, magnesium ions.
Performance measurement results are reported in FLOP/s, where the number of floating point operations per cycle for each system was determined by the perfex monitoring utility. Once equilibrated, restart simulations were run for 200 cycles, load balanced by complete reassignment based on measurement over cycles 100-200, and load balanced again by refinement based on measurement over cycles 300-400. The execution time was measured over steps 2100-2440 after restart. For single processor simulations, restart runs were performed starting at step 2020. The execution time was measured over steps 2100-2440.

Performance Results
Near-linear scaling is obtained as a result of the large ratio between the compute time and the communication time. Figure 2 displays the performance as a function of the number of processors.
Our peak performance was achieved for the whole ribosome system (N atoms ~ 2.64x10 6 ) at 437 GFLOP/s. The relative performance as a function of processors has a speed-up of ~867 for the larger whole ribosome system (~84.67% efficient), where efficiency = speed-up/N procs , speed-up = t 1 /t Nprocs , t 1 is the execution time on a single processor without MPI and t Nprocs is the execution time on N procs processors (Fig. 2). Studying the scaling as a function of the number of atoms, we observed that the "turnoverpoint" in speed-up occurs near N procs =256 and 5.73x10 4 < N atoms < 9.22x10 4 . We also observed that NAMD requires at least 2 GB/process for simulations with N procs > 768 and N atoms > 2x10 6 .

Science Results
The rate-limiting step during ribosomal decoding for the acceptance of cognate tRNAs is accommodation, where the aminoacyl tRNA moves from its partially bound to fully bound state on the ribosome. 19 We have implemented the targeted MD algorithm 6; 20; 21 in explicit solvent to simulate the conformational change that occurs during accommodation. This method gradually reduces the rootmean-squared distance (RMSD) of the complex to the final state (determined experimentally) while allowing thermal fluctuations of the structure at any given RMSD. The targeted molecular dynamics simulations produce stereochemically feasible pathways with candidate interactions between the tRNA and the ribosome that can be tested via site-directed mutagenesis. The simulations were described in detail previously and will be summarized here. 5 To simulate single accommodation events, as opposed to spontaneous rates, several simulations (eight in total) were performed with durations of 2 ns each. The sensitivity of the results to time scale was tested by performing simulations of 1 ns and 4 ns. During the simulations, the body of the aminoacyl-tRNA relaxes from the kinked, partially bound state 22 to the native-like, fully bound state. An accommodation wall region is defined by the motion of the tRNA along the large subunit of the ribosome (Fig. 3), as the tRNA sweeps over the large ribosomal subunit during accommodation. Specifically, helix LH89 (Fig. 4) is positioned to act as a guide rail, ensuring that fluctuations in the tRNA elbow angle are sufficiently small to allow the 3'-CCA end to reach the peptidyl transferase center. The large subunit helices LH38 and LH69 are positioned to prevent the aminoacyl-tRNA from over-shooting its final equilibrium position.
Interestingly, a universally conserved loop (on LH92) blocks the entrance of the 3'-CCA end of the aminoacyl-tRNA into the ribosome's peptidyl transferase center. The CCA end flexes backwards relative to the motion of the tRNA body as the CCA end encounters the A-loop. Both the CCA end and the A-loop flex, allowing the CCA end to enter the peptidyl transferase center (Fig. 4). The flexing of the CCA end constitutes a second flex region on the tRNA, in addition to that discovered by Frank and co-workers. 22 The 4 ns validation simulation displayed similar behavior, demonstrating that 4 ns simulations offer little new information.

Discussion
With the explosive growth of computational power in recent years, biological systems are becoming more and more physiologically relevant. The ribosome simulations represent a significant increase in system size relative to previous simulations. Historically, biomolecular dynamics simulations originally simulated very short timescale dynamics (~ 10 ps) of small proteins (bovine pancreatic trypsin inhibitor 'BPTI-1' ~500 atoms) without solvent molecules due to limitations in compute power. 24; 25 Increases in computing power allowed inclusion of solvent for these systems (bovine pancreatic trypsin inhibitor 'BPTI-2' ~3100 atoms). 26 Fast multipole algorithms then achieved simulation sizes of 1.26 x10 4 (photosynthetic reaction center of Rhodopseudomonas viridis, 'RHOD'), 27 2.4x10 4 (POPC lipid bilayer patch), 28 3.6x10 4 atoms (estrogen receptor binding domain plus DNA complex) 29 . The fast multipole method was combined with a multiple-time-step method to simulate systems of 3.6x10 4 atoms. 30 A significant improvement in parallelization with dynamic load-balancing enabled a simulation of the solvated acetylcholinesterase dimer ('AChE') in absence of long-range forces. 30; 31 Finally, more sophisticated dynamic load balancing, replacing spatial domains by meta-domains based on compute-load, has produced simulations of > 3x10 5 atoms (f1ATPase macromolecular complex, 3.26x10 5 atoms). 32 The ribosome system (2.64x10 6 atoms) is approximately six times larger than the previous largest system simulated (DOPC) 33 . The ribosome simulations have not only helped to elucidate a crucial molecular mechanism for gene expression, but also open the door for simulations of other large molecular machines important for gene expression and drug design. Antibiotic drugs are less than 1/1000th the size of the ribosome, yet have the ability to diffuse into the most critical sites of this huge molecular factory and grind the inner workings to a halt. The simulation maps out regions where we expect the transfer RNA to rub against the ribosome. In particular, we have identified a corridor inside the ribosome that the transfer RNA passes through for the decoding to occur. This corridor is constructed almost entirely of universal bases, implying that it is evolutionarily ancient and critical for protein synthesis. The corridor represents a new region of the ribosome containing potential new antibiotic targets. Designing an antibiotic drug molecule that binds to one of these regions may degrade the function of the ribosome in harmful bacteria. The results represent a proof-of-principle for locating potential antibiotic targets on large molecular complexes.