De Novo modeling in cryo-EM density maps with Pathwalking.

While the original implementation of the Pathwalking protocol was far better than purely manual model building, it did require a considerable amount of human interaction. This means the ultimate accuracy of the model is highly dependent on the skill of the person doing the modeling. The main interactive steps in the original version of Pathwalking includes identifying secondary structures using the SSEHunter [ 38 ], setting a threshold for the density map before pseudoatom placement, and finally, evaluating the identified paths and setting constraints. Even with the semi-automated approach of Pathwalking, these manual steps could be relatively time-consuming for large proteins. This is further complicated by the fact that multiple models may be required for large structures, making accurate map segmentation a critical issue. In the new version of Pathwalking, we have addressed these problems by introducing additional quantitative assessments and automation. A fully automated modeling process can be performed using the command: e2pathwalker_auto .py density_map . mrc -- natoms =N -- denthr =T Here N is the number of residues in the protein and T is the threshold of the iso-surface value in the density map. The program will calculate a default threshold value if --denthr=T is not provided, but a carefully selected threshold value that separates features is still recommended especially at lower resolution. This script will run the programs in our new Pathwalking workflow using the default parameters, and generate the output for each step. It should be noted that the single e2pathwalker executable does not assign the sequence to the trace, rather the it is simply an enumerated poly-alanine C-alplha backbone trace. The EMAN2 companion program, e2seq2pdb.py, will do both a forward and reverse sequence assignment to the trace. It is also important to note that, like the original version of Pathwalking, no explicit sequence information is utilized. As such, no registration of strucural details such as SSEs or bulky sidechains are used in modeling. However, this type of information can be used to independently validate models generated from Pathwalking. While the fully automated protocol can generate decent models when the input density map is well resolved, users may also opt to run each program in the protocol step by step to adjust parameters and access advanced options. This is particularly useful when attempting to apply additional structural constraints. All output models are in standard pdb format so users can view the results in molecular visualization software, such as UCSF’s Chimera [ 40 ]. A detailed tutorial can be found online at http://blake.bcm.edu/emanwiki/Pathwalker . In the following sections we detail the algorithmic and automation improvements in each step of our new protocol ( Fig. 1 ). 2.1. Pseudoatom placement In the first step of the Pathwalking protocol ( Fig. 1B ), the density map is first mapped onto a Cartesian grid with a grid spacing of 3.8Å, the average distance between C α atoms in peptides. While there is no reason to expect that we are going to identify precise C α locations rather than random locations along the protein fold, this scale is representative of the level at which we need to characterize the structure. N pseudoatoms are first placed in the center of cubes with the N highest mean densities, where N is the number of amino acids in the protein. K-means clustering is then performed using these points as initial cluster centers and all the points above a certain threshold in the density as samples. The threshold is set to be the lowest mean density of the N grid points by default. The k-means method has been modified so that the new center positions in each iteration are a density weighted average of the member points rather than an unweighted average. This pseudoatom seeding process terminates when the relative motion of cluster centers falls below a small threshold value. 2.2.

Path detection

After pseudoatom generation, path detection is performed using a TSP solver ( Fig. 1C ). Pathwalking is compatible with both the Concorde [ 36 ] and LKH [ 37 ] TSP solvers. Here, the paths are evaluated by biochemical constraints and the fitness to the density map. A distance matrix is created in which the distance between two pseudoatoms is expressed in terms of the cost function to be an actual C α -C α “bond”. This cost function is calculated using the distance between the two atoms expressed as a deviation from the optimal 3.8Å C α -C α distance and the density value between the atoms. (1) { ( d - 3.78 ) 2 , d < max _ dist , ( d - 3.78 ) 4 , d ≥ max _ dist + w ∗ ( 1 / avg _ den ) ∗ { 2 , avg _ den > den _ thr , 1 , avg _ den < = den _ thr Here d is the distance between the two pseudoatoms, and avg_den is the average density of all voxels on a straight line between them. max _ dist and den_thr are the maximum distance and density threshold that can be set by the user. This distance matrix is then used as the input for the TSP solver, which attempts to find the lowest cost path of the protein through the density map. 2.3.

Path refinement

A particularly complicated problem for Pathwalking is the fact that different features in the map may have significantly different relative density values. If a single density threshold is used, the path can be either disconnected or highly ambiguous. This is particularly evident with SSEs at near atomic resolution; at higher thresholds helices are typically well resolved but coils and strands are usually not discernible. When lowering the density threshold, SSEs often merge together, again making it hard for accurately discriminate density features. To overcome these issues, we implemented an iterative path refinement step based on automatic SSE identification. Secondary structure elements are firstly identified based on the initial path from Pathwalking as well as geometric criteria which are further detailed in the following sections ( Fig. 1D ). Pseudoatoms in these secondary structures are then fixed to their ideal positions. After this, the density corresponding to these SSEs are removed from the density map and the remaining pseudoatoms are seeded again on the modified density map ( Fig. 1E ). A new path is then computed based on these pseudoatoms, including those “locked” from the previous round of secondary structure elements identification. By repeating this process, more and more secondary structure elements are identified, in addition to refining both pseudoatom placement and path identification. Model refinement stops when there are no more SSEs found in an iteration ( Fig. 1F ). 2.3.1. Secondary Structure Identification: α-helices After an initial path is generated, SSE detection first identifies long helical/cylindrical structures along the path. Here, helical shape is defined as a subpart of the path that has all its pseudoatoms within a range of distances to the central axis of these points. For each pseudoatom, we sum the number of pseudoatoms extending in both directions that maintains an α -helical profile, and set this value as the helical score of this pseudoatom ( Fig. 5C ). Potential helices are identified using this score. To confirm these putative helices, the algorithm generates ideal pseudoatoms for a helical path at these sites and fits them to the density map; potential helices can be identified by the fitness of ideal helices. For the confirmed helices, pseudoatoms in the helices are replaced by the ideal α helices. The length and position of these helices are then optimized based on the density map. After the helix fitting step, a density map with helical density removed is generated. The program then fix the pseudoatoms on helices and re-seeds the rest pseudoatoms. The path tracing and helix detection is repeated until no new helices are found. 2.3.2. Secondary Structure Identification: β-sheets While helix identification can generally be performed at resolutions as low as 6–8Å in cryo-EM density maps, β -sheet tracing is typically best suited for maps at better than 4.5Å resolution, where the strands in the sheet begins to resolve. For sheet identification, we calculate a score for each pseudoatom indicating the cost function of it belonging to a β sheet ( Eq. 2 ). The score is calculated based on the angle between the bonds connecting the pseudoatom and the other bonds near the pseudoatom in the 3D space, so that the score for the pseudoatom is higher when it is in a region where the bonds tend to be parallel to each other, suggesting a potential β sheet ( Fig. 5A ). (2) ∑ ( k nearest atoms ) PDF Gaussian ( distance , sigma ) ∗ abs ( cos ( angle between bonds ) ) Then, the score is smoothed along the backbone chain using a Gaussian moving average, and the algorithm searches for the longest continuous chain in which the pseudoatoms have a score higher than a threshold (mean+std by default). Once a chain has been identified, it is elongated from both ends based on the position of nearby pseudoatoms. After that, the score for all the rest of the pseudoatoms is calculated again while the bonds near and parallel to the existing chains are weighted higher. The process is performed iteratively until the score or length of the next new chain is below a certain threshold. 2.4. Path modification After the refinement completes, the program searches the final path for “crossover bonds”, points along the path that should be continuous but exhibit a marked departure in their current direction and forming non-protein-like structure. In most cases, these crossover bonds can be fixed by simply switching the order of residues in the backbone path. If the reordered path still does not meet the geometry criteria, the related pseudoatoms can be randomly shaken in the local region to obtain a better path. In addition to these non-protein like features, pseudoatoms can be placed in large sidechain densities in the pseudoatom generation step. These pseudoatoms can be identified by the density around the atom, length and angle of bonds and the relationship with the main chain. These sidechain atoms are then eliminated from the backbone model. Along with the aforementioned Pathwalking utilities, we have provided a supplemental script for common structural modification in C α model. In this step, users should evaluate their model and refer to the program documentation regarding these types of modification steps. 2.5. Termini determination In addition to iterative path improvement, a new step has been incorporated into Pathwalking to better localize the backbone termini. In most instances, the position of terminal pseudoatoms in the density map is unknown before modeling, which can also be complicated by the fact that in many cases the C and N termini are near each other in 3-D. The new algorithm identifies the most likely termini based on the geometry of the pseudoatoms. Two phantom points are initially placed in the model as start and end points, and the distance to each pseudoatoms is set to 0. Since the TSP solver is designed to find the path with the minimum total length, it will connect a pseudoatom to this phantom termini when it is the most unlikely to be connected to any other atoms. Finally, the program removes the phantom points and assigns the pseudoatoms directly connecting them as the termini. Although the program can detect the termini automatically, the direction of the path is still ambiguous, which can be resolved by the user using sequence and secondary structure information. As an added feature, this approach may also be used to model multiple subunits in the same density map that cannot be properly separated/segmented. Again, an extra phantom point is added for each subunit, where adding the phantom point is equivalent to allowing a gap in the calculated path. The TSP solver will find the most separable bond to put the gap in order to minimize the global distance. As such, adding phantom points in the seeding step allows for the tracing additional subunits. In addition, this technique can also be used to build models on density maps with poorly resolved regions. Users can treat a poorly resolved density map as multiple subunit so the program can automatically avoid density breaks/gaps, and build the backbone for only the well resolved parts of the map.

Supplementary Material supplement