Fitting Low-Resolution Cryo-EM Maps of Proteins Using Constrained Geometric Simulations^☆

Constrained geometric simulations

Our multi-scale flexible fitting tool is based on the recently-developed constrained geometric simulation algorithm FRODA (Framework Rigidity Optimized Dynamical Algorithm) (20), which is implemented as a module of the FIRST (Floppy Inclusion and Rigid Substructure Topography) software package. FRODA is a Monte Carlo-type algorithm in which diffusive motion in biomolecules is simulated as motion of rigid clusters within the biomolecule by rotation about bonds, and steric clashes are avoided. In short, FRODA models a biomolecular structure as a collection of atoms, each of which belongs to one or more rigid “ghosts” that serve to guide the motion of the atoms. At each step, the atoms are thrown from their ghosts by ∼0.1 Å, and the atoms and ghosts are iteratively re-fit to one another until a new valid conformation is found. A more detailed description of the FRODA method can be found in the Supplementary Material accompanying this article or in Wells et al (20). The rigid clusters are identified using FIRST, which identifies rigid and flexible regions in the structure using a graph-theoretical algorithm known as the pebble game. A detailed description of FIRST and the pebble game approach can be found in (23–28).

In this study, the direction of the atomic throws (and therefore of all subsequent motion) was random, being guided only by the rigid ghosts and by changes in the scoring function. The throw directions could also be chosen in a nonrandom way, such as movement along an elastic normal mode or bias toward a particular target structure.

Scoring

To use FRODA to fit an atomic structure to a density map, a constrained Monte Carlo simulation is used to maximize the real-space correlation coefficient

which quantifies the degree of overlap between the experimental density ρ^exp and the simulated density ρ^sim. This avoids having to normalize the intensity of the experimental density. The correlation coefficient ranges from 0 (no overlap) to 1 (perfect overlap). At each step in the simulation, ρ^sim is calculated based on the current FRODA conformer, and the correlation coefficient is evaluated. The decision to keep or reject the new conformer is made based on the Metropolis criterion (29), with the correlation coefficient playing a role analogous to the energy in a standard Metropolis simulation. In other words, a random move that improves the correlation between ρ^exp and ρ^sim (i.e., ΔC > 0) will always be accepted, and a random move that decreases the correlation (ΔC < 0) will be accepted with a probability of e^ΔC/s, where the annealing scale s plays a role similar to that played by the thermal energy kT in the usual Metropolis algorithm, and for simplicity we will refer to s as a temperature. To avoid having to carefully select an appropriate value of s, we used a dynamic temperature adjustment scheme that adjusts the value of s during the course of a simulation to maintain a desired acceptance rate, usually 50%.

Generating density maps

To calculate a simulated density map based on a high-resolution structure, we start with scattering factors that were calculated for naturally-occurring elements using a multiconfiguration Dirac-Fock code (30), and conveniently expressed as a sum of a few Gaussians. Electron densities for individual atoms were obtained from the scattering factors using three-dimensional Fourier transforms (Fig. 1), and the densities of each individual atom in a high-resolution structure were summed to create molecular density maps.

Simulating experimental data

To make these density maps more similar to those obtained experimentally, we convolved this ideal map with a Gaussian function, yielding noise-free maps at finite resolution. To obtain physically-relevant values for the resolution, we fitted a three-dimensional Gaussian to the function obtained by starting with a single real-space δ-function (corresponding to a constant in k-space) that is Fourier-transformed back into real space, neglecting any contribution from k > k_max. The width of the resulting real-space peak was used to set up a direct correspondence between the resolution obtained in a diffraction experiment (2π/k_max) and the width σ of the Gaussian used for convolution. This procedure yields results that look qualitatively similar to experimental electron density maps at a given resolution. The choice of a Gaussian kernel for the broadening of the density map was not inevitable; other tools for the generation of density maps from crystal structure coordinates allow for the usage of a variety of kernels (31). The Gaussian kernel was chosen primarily because of the fact that the approximate scattering factors (30) were expressed in terms of a sum of Gaussians, which made Gaussian convolution particularly convenient. Other more complex schemes could also be incorporated into this method.

The presence of noise in experimental data is inevitable, and can provide serious limitations to the accuracy of flexible fitting methods. To assess the ability of this method to cope with noisy data, we also generated theoretical density maps with a controllable amount of noise. If the density measured in a particular voxel is N, then we assumed that the noise level in that voxel will be proportional to the square root of N. More precisely, the noise was modeled using a Gaussian distribution with σ proportional to the square root of N; the constant of proportionality can be varied to control the amount of noise present in the simulated data, and an overall signal/noise ratio can be calculated after the noisy map has been created.

Structure preparation

Simulations with FIRST/FRODA require that hydrogen atoms be present in the structure. Hydrogen atoms were not present in any of the structures obtained from the PDB (Protein Data Bank, http://www.pdb.org), and were added using Reduce (available at http://kinemage.duke.edu) (32). FIRST/FRODA imposes rather stringent requirements for stereochemical quality, and the crystal structure for GroEL (PDB accession code 1KP8) contained some steric clashes that prevented FRODA from starting successfully. To resolve these clashes, we subjected the GroEL structure to a conjugate gradient minimization using the CHARMM force field in the molecular dynamics package NAMD2 (33). To avoid any large-scale changes during the minimization, each atom was tethered to its initial position by a harmonic restraint with spring constant k = 5.0 kcal/mol/Ų.

Simulation details

The apo-GroEL simulation was run with a resolution of 10.3 Å in three separate phases. In the first phase, every individual rigid body was thrown in FRODA, allowing for optimal sampling of large-scale, collective motions. In each step of the second and third phases, each rigid body was thrown with a probability of 0.1 and 0.01, respectively, whereas those rigid bodies not being moved actively made only minor adjustments in response to the movement of their neighbors. This allowed us to sample localized motions more effectively in the later stages of the fitting, as required for an optimal fit. The ATP-bound GroEL simulations were run at a resolution of 14.9 Å, using the same three-step scheme to decrease throw probability. Because these simulations involved larger motions and therefore a potentially greater computational cost, we increased the simulation speed by increasing the spacing of the map grid during the first two phases. The original cryo-EM data was on a grid with a spacing of ∼1.9 Å; by ignoring every other point we could increase the spacing to ∼3.7 Å and decrease the information content of the map by 1/8, resulting in a significant speedup. This speedup, of course, comes at the cost of a lowered signal/noise ratio; for this reason the full grid was used during the final phase of the simulation. Each simulation phase was stopped when the correlation score stopped improving. A fourth simulation phase (in which the throw probability would have been further decreased) was not used because the improvement between the second and third phases was small. A similar sort of procedure could be applied for other systems, but will probably be most effective in large structures with many separate rigid clusters.

Resolution and grid spacing

As a first test case, we selected adenylate kinase (ADK). X-ray structures of ADK have been solved in both its closed form (PDB accession code 1AKE) and its open form (PDB code 4AKE). The root mean-square difference between these two structures is relatively large (∼7.2 Å) and yet the system is small enough (214 residues) to allow for fast computations. Most of the conformational change between the open and closed forms involves motion about a hinge domain (Fig. 2 A) and the improved fit to the density map can clearly be seen.

In any constrained geometric simulation using FIRST/FRODA, the choice of the noncovalent constraint network is crucial (34). If too many noncovalent constraints are included, the protein will be overly rigid and flexible fitting will become impossible. If too few are included, the protein will be overly flexible and the geometric integrity of secondary structure elements will be lost, eliminating some of the advantages of the constrained geometric approach. In general, however, including fewer constraints leads to a higher final value of the correlation coefficient, so that the choice of the noncovalent network involves a compromise between detailed fitting and chemical integrity. In this case, we have chosen a noncovalent constraint network that maintains most of the α-helices in ADK as rigid units, whereas most of the random coil regions remain flexible. Based on a direct all-atom targeting (20) using this constraint network, the closed conformation should be able to come within about RMSD ∼1 Å of the open structure without violating any of the noncovalent constraints included.

This method requires comparison of the target density map with a simulated map calculated based on the moving structure. This map is calculated by placing calculated electron densities on each nonhydrogen atom in the structure, adding them together, and using Gaussian convolution to broaden the map to an arbitrary resolution. To test the effects of resolution on flexible fitting, we have fitted the closed form of ADK to density maps generated using the open form, with the resolution ranging from 0 to 13.3 Å in 0.67-Å steps. The fitting in the opposite direction, starting from the open crystal structure and fitting to the closed form, was also tested and works equally well. The correlation function is evaluated on a regular cubic grid, with the choice of grid spacing reflecting a balance between speed (favored by a large grid spacing) and accuracy (favored by a small grid spacing). For each resolution, we repeated the calculation with 1, 2, 3, and 4 Å grid spacings. For the fitting at each resolution and grid spacing, we measured the improvement of the correlation during the course of the simulation, as well as the root mean-square deviation (RMSD) from the open form of ADK. The RMSD-to-target is a particularly important measure, because it indicates how well the targeting routines are able to reproduce the “real” structure behind a low-resolution map. This sort of validation is only possible for theoretically-generated data.

The results show that, for nonpathological cases, a roughly linear relationship can be established between the accuracy of fitting and the resolution of the density map. For most simulations, the RMSD and the correlation score roughly track one another, and each reaches an asymptotic value eventually (Fig. 3). The asymptotic value of the RMSD will tend to be higher for lower-resolution maps, and lower for high-resolution maps. This result is what we would expect intuitively. Higher-resolution maps have higher information content and are able to constrain more effectively the final fitted structure (Fig. 4).

Another important result concerns those cases in which successful fitting cannot be carried out. The linear behavior described above breaks down when we attempt to fit a map with a resolution that is comparable to the grid spacing. In these cases, the simulations do not converge, and the simulated structure does not move toward its target. Below this threshold, however, the choice of grid spacing seems to be of little importance. We can therefore obtain optimal fitting performance by choosing a grid spacing that is somewhat smaller than the length scale imposed by the cryo-EM resolution.

The same set of trial simulations was conducted using lactoferrin (Fig. 2 B). Lactoferrin is somewhat larger than ADK, with 691 residues, and undergoes a conformational change of ∼6.5 Å. Plots of the all-atom RMSD and correlation are contained in the Supplementary Material accompanying this article. The results obtained with lactoferrin were qualitatively similar to those with ADK, with the fitting process failing as the grid spacing approaches the resolution length scale. This fitting method should therefore be used with caution on maps in this resolution range. The difficulty is associated with the highly structured density at resolutions <5 Å, where the density is significantly higher at the center of each atom. This leads to a “corrugation” effect in which it is difficult for the trial density to glide over the target density. At resolutions greater than ∼5 Å (the relevant region for cryo-EM), the density is much smoother and the problem disappears.

Noise tolerance

To effectively fit atomic models to real cryo-EM data, it is essential that the method be able to tolerate noise in the map. To evaluate the robustness of the fitting algorithm in the presence of noise, we have used theoretical density maps to which noise had been added, as described in Materials and Methods. ADK was fit at 3.3 Å resolution with signal/noise (S/N) ratios ranging from 10 to 1 (Fig. 5). High levels of noise were detrimental to both the speed and accuracy of fitting, but modest noise levels (S/N > 1) seemed to have very little impact on the fitting results. This gives some confidence that this method is robust enough to withstand reasonable amounts of experimental noise.

Effects of missing density

Another potential issue that could arise when working with experimental cryo-EM maps is that regions of density may be missing in the experimental data. To test the effects of missing density on our fitting method, we constructed a map at 6.6 Å using the open form of ADK, except that the C-terminal helix had been removed (Fig. 6). We then attempted to fit the complete closed form of ADK into this partial map. The result was that, rather than dangling outside, the C-terminal helix was pressed into the density, distorting the surrounding region and forcing some other secondary structure elements slightly out of the way. This result could have been anticipated and is what we would expect intuitively. This also underscores the importance of knowing what is present in an experimental ED map before attempting flexible fitting. If some portions of the structure are missing, the resulting fit could be seriously distorted. One possible way to avoid this problem would be to choose parameters that give larger rigid clusters, effectively decreasing the number of degrees of freedom and avoiding internal distortions. ADK has a high degree of intrinsic flexibility, however, and drastic limitation of the degrees of freedom would be physically unmotivated.

If concerns exist that the structure being used for fitting may contain more (or fewer) atoms than the molecule used to generate the density map, one appropriate solution might be to subtract a map generated based on the final fitted structure from the target map. Such a difference map would show significant density only in those areas where local errors existed in the fit. The structures shown in Fig. 6 suggest that errors introduced by missing density will be localized around the faulty region in the density map, so a large amplitude in the difference map could point to problem areas in the density map or the atomic model. It should also be emphasized that these difficulties do not nullify one of the major advantages of real-space fitting (3), which is that single members of a multi-component complex can be fitted into a map of the complete complex even if their neighbors are missing.

Fitting experimental density

To test the ability of our method to work with experimental cryo-EM data, we studied a conformational transition in the bacterial chaperonin GroEL. GroEL is a multimeric complex formed from two heptameric rings, with a total of 7350 residues. Each GroEL monomer can be divided into three domains (35), the equatorial domain (residues 1–133 and 409–526), which forms strong contacts with monomers in the opposite ring, an intermediate domain (residues 134–190 and 377–408), and an apical domain (residues 191–376) that interacts with the co-chaperonin GroES. Misfolded or unfolded proteins bind nonspecifically to the hydrophobic interior of unliganded GroEL. On binding ATP, GroEL undergoes a conformational change that enables binding to GroES. This traps the misfolded substrate protein in an environment more conducive to proper folding; for a recent review of GroEL function see Lin and Rye (36).

GroEL has been the subject of numerous structural studies using both x-ray crystallography and cryo-EM. We have used a 2.0 Å crystal structure (37) of GroEL-(KMgATP)₁₄, which was prepared for the simulations as described in Materials and Methods; this was the highest-resolution structure of the complete GroEL complex that was available in the Protein Data Bank. This high-resolution structure was then fit into density maps of GroEL both before (apo-GroEL) and after (GroEL-ATP) the ATP-induced conformational transition (38). Some of the results of this fitting are depicted in Figs. 7 and 8 and are described in detail below. The “intact domain” and “secondary” curves on these figures refer to two different fitting schemes using the GroEL-ATP map, and the “apo-GroEL” curve refers to the fitting results from the uncomplexed GroEL map.

To assess the usefulness of our flexible fitting method, it is necessary to compare our results with what can be obtained through conventional rigid-body fitting. The apo-GroEL map was very similar to the GroEL crystal structure, and we were able to rigidly align the crystal structure with the map without modification. The crystal structure and the map were aligned by hand, and then the fit was optimized by maximizing the correlation function with only translational and rotational degrees of freedom. The best fit between the crystal structure and the apo-GroEL map gave a correlation coefficient of C = 0.8943.

We then carried out flexible fitting on the apo-GroEL map. The conformational change seen was very slight and involved an improvement of the correlation coefficient to C = 0.9013. Because the map was so similar to the initial conformation given by the crystal structure, we maintained the network of noncovalent interactions that were present in the initial structure. This noncovalent bond network maintained each individual α-helix as a separate rigid unit, and these rigid units were connected extensively to one another by flexible coil regions and sidechains, which served to maintain the overall shape of the complex. Fig. 7 shows a plot of the RMSD between residues in the crystal structure and the fitted GroEL structures, averaged over the 14 subunits in the complex. The apo-GroEL fitting showed most of its motion in the apical domain, somewhat less in the intermediate domain, and rather little motion in the equatorial domain.

To understand the intersubunit variations within the fitted structure of the multimeric GroEL complex, we measured the structural variation within the ensemble of 14 monomer structures at each residue (Fig. 8). For each residue, the RMSD was calculated and averaged over all unique intersubunit pairs. Carrying out this analysis for the crystal structure gives us a baseline of variation among nominally identical subunits, around 0.5 Å. For the fitted structures, the plot is extremely rough, and we can gain more insight by smoothing it out so that each point represents the average value for 50 adjacent residues. The most striking feature of the smoothed curve for apo-GroEL is its similarity to Fig. 7. The residues that experience the most motion during the fitting process are also those that show the greatest variation when fitting is complete.

The density map of GroEL complexed with ATP showed a significant conformational change from the apo-GroEL form. In their work describing this structural change, Ranson et al. (38) obtained an approximate fit to the density map by breaking the individual GroEL monomers into distinct domains, then fitting each of these domains separately into the map. This allowed us to estimate that GroEL makes a highly concerted conformational change with an RMSD of ∼6.3 Å. This change involves an axial extension of ∼14% (from ∼150 to 161 Å), a 4% increase in the radius of gyration (from 63.9 to 66.4 Å), and a counterclockwise outward rotation of the apical domains of each subunit. The best fit we could obtain between their domain-fitted structure and the cryo-EM map gave a correlation coefficient of C = 0.9507.

The same flexible fitting procedure used with the apo-GroEL map was then carried out on the GroEL-ATP map. The 2.0 Å crystal structure was again used as a starting conformation, but the larger conformational change involved made the choice of noncovalent bond network less clear than in the apo-GroEL case. We carried out the fitting with two different sets of noncovalent interactions. In the first, which we will refer to as the intact domain case, we calculated the noncovalent interactions for both the initial crystal structure and the structure obtained by rigid-domain fitting, and kept the interactions that were present in both cases. Any interactions that would be disrupted during the structural transition were therefore omitted, and those that are present at both endpoints were assumed to be present throughout the transition. In the second case, we maintained only those hydrogen bonds internal to secondary structure elements (i.e., α-helices and β-sheets) as identified in the crystal structure using STRIDE (39). The noncovalent interactions included in the secondary structure case were a subset of those included in the intact domain case, and generally allowed for more movement of the complex due to the smaller number of constraints. In the intact domain case, we were able to improve the correlation coefficient from an initial value of C = 0.9253 to a final value of C = 0.9571. In the secondary structure case, the correlation reached a final value of C = 0.9687, indicating that the increased mobility allowed us to find a better fit to the density map. Both flexibly-fit results obtained a better correlation than the rigid-body result of C = 0.9507.

The motion experienced by each residue during the fitting process is shown in Fig. 7. The motion in the secondary structure case shows several peaks that are not present in the intact domain case, especially in the equatorial domain. Many of these arise from the fact that the hydrogen bonds included were often not sufficient to maintain strict rigidity of secondary structure elements. For example, the last turn of an α-helix often was not part of the same rigid cluster as the rest of the helix and unwound in several cases, leading to a large RMSD for a few adjacent residues. The smoothed version of these curves show that both GroEL-ATP fitting simulations experienced significant motion in the apical domains, a second prominent peak in the N-terminal portion of the intermediate domain, and significantly less motion in the equatorial domains.

Analysis of the intersubunit variation for the results of the GroEL-ATP fitting simulations (Fig. 8) shows that, for both fitted structures, the largest variability is present in the apical and intermediate domains. The situation is similar to that seen with apo-GroEL; the residues that experience the most motion during fitting show the largest variability in the results. This variability is not a result of asymmetry in the cryo-EM map, which was refined with explicit sevenfold symmetry. Instead, these deviations from symmetry give a sense of the ambiguity inherent in fitting an atomic model to low-resolution data and allow us to estimate the precision of the resulting fit. In principle, it would be possible to force FRODA to make moves that are consistent with this sevenfold symmetry, thus producing a perfectly symmetric result. This was beyond the scope of this study.

Careful examination of the fitted structures produced by FRODA shows that a large number of intra- and intersubunit contacts were broken and re-formed during the course of the fitting simulations. In general, more contacts were broken than were re-formed; this stems in part from the more extended structure of the ATP-complexed GroEL and in part from the fact the FIRST does not optimize the formation of new interresidue contacts. New contacts were generally formed between the same regions of the protein subunits in which contacts were lost, indicating that the structural transition in GroEL involves deformations of the individual subunits more than the formation of entirely new contacts.

One can also compare the disruption and formation of protein–protein contacts in the flexibly fit structures with the contact changes in the structure obtained by rigid fitting. Relatively few intraresidue contacts are broken or re-formed in the rigid-fitting case relative to the flexible fitting. This is to be expected because the individual domains were treated as rigid units in the rigid fitting, but internal motions were allowed in the flexible fitting. The more intriguing result was that the changes in intersubunit contacts are almost identical for the two fitting methods. This suggests that, due to length scale limitations imposed by the map resolution, the constrained geometric fitting did a better job of capturing the relations between subunits than the internal dynamics of individual subunits.