CS-ROSETTA

2.1. Energy

The main driving force of Rosetta is its energy function (Alford et al., 2017; Leaver-Fay et al., 2013, 2011; O’Meara et al., 2015). The energy function quantifies the static conformational state of a protein, which can be utilized to distinguish native-like structures from alternative folded and unfolded proteins (Alford et al., 2017). This function is a linear combination of terms that capture physically realistic structural features. Specifically, these terms describe a set of potentials that report on inter/intramolecular interactions, solvation, hydrogen bonding, electrostatics, side chain rotamer orientation, and dihedral preferences. During the structure modeling process, Rosetta can modify the contribution of each energy term using predefined (score3. wts, talaris2013.wts, ref2015.wts, and others) or customizable weights. This procedure may help steer the calculations toward a global minimum in the energy landscape for proteins of interest. We use the terms, energy function, energy, score, and score function synonymously in this chapter.

2.2. Pose

In Rosetta, a biomolecular system is represented using a C++ container class known as a Pose (Leaver-Fay et al., 2011). Pose holds information on the chemical composition (protein amino acid sequence), degrees of freedom (variable factors of a conformation), kinematics, and energy of the molecular system. For protocols such as ab initio (where structures are modeled solely from amino acid sequence information) (Bradley, Misura, & Baker, 2005), degrees of freedom of a Pose can be derived from an extended polypeptide chain. In alternative protocols, degrees of freedom are derived from a reference structure input manually by the user. Moves (see Section 2.5) are applied to a Pose during the process of sampling diverse conformations while searching for the native (or global optimum with respect to the score function). If the conformation of a Pose is modified using any mover, then its score function is updated accordingly.

2.3. Fold Tree

Each amino acid within a protein chain in Rosetta is characterized by two degrees of freedom known as internal coordinates (given by the torsion angles, φ and ψ) (Bradley & Baker, 2006; Leaver-Fay et al., 2011).

A small perturbation in the local conformation of a protein chain, induced by a mover (see Section 2.5), leads to the propagation of changes globally, called the lever arm effect. To eliminate this effect, a protein chain is denoted as a directed, acyclic graph called the fold tree, which is attached to the Pose. In a fold tree, the nodes indicate amino acids of a protein and the edges represent peptide bonds, which imply that any residue i is connected to its neighbors i − 1 and i + 1 by edges. Alongside peptide bonds, through-space bonds (or long-range connections) can be added as edges in this tree without the violation of acyclicity property of the graph. To preserve this property, an insignificant peptide edge is deleted, and this introduction of gap in the fold tree is called a chain-break or a jump. These edges govern the propagation of conformational changes during a move. The fold tree with discontinuous peptide edges is referred to as broken-chain fold tree.

2.4. Centroid and Full-Atom Representation

Protocols that require extensive sampling such as ab initio (Bradley et al., 2005) switch to coarse-grained representation, where each residue in a protein chain contains five backbone atoms (N, H, C^α, C, O) and one side chain centroid atom (CEN) characterized by its amino acid type. This simplified representation is also referred to as low-resolution mode. The use of centroid mode during early stages of structure determination boosts the performance of sampling methods by reducing the degrees of freedom introduced through side chain dihedral angles (Shen et al., 2008).

On the other hand, a full-atom representation (also referred to as high-resolution mode) includes backbone atoms in addition to all side chain atoms. The full-atom representation mode is typically used during the refinement stages of structure determination protocols. In this chapter, we interchangeably use terms for low-resolution with centroid and high-resolution with full-atom.

2.5. Mover

A mover is an entity that can modify or analyze any given Pose (Leaver-Fay et al., 2011). CS-Rosetta’s structure modeling protocols are built using several movers, which are referenced throughout this manuscript.

2.6. Topology Broker

Topology broker is a framework that enables construction of complicated simulation systems (Bender et al., 2016; Porter, Weitzner, & Lange, 2015).

Broker guides the resampling protocols used for structure determination by providing them information about (i) movers to apply, (ii) chain-breaks to add, (iii) representation modes to switch to, and (iv) construction of refined fold tree. Any knowledge about the protein target including structural restraints (from NMR data or any other source) or rigid regions in homologous structures from databases such as PDB (or protein data bank) (Berman et al., 2000) is retained by the broker in the form of claimers. Specifically, claimers define the accepted and unaccepted degrees of freedom for a given Pose. For instance, if rigid regions are defined for a Pose, then the broker will not accept any changes in those regions. Therefore, broker allows rapid protocol development, since the methods do not have to keep track of constraints that must to be satisfied by the system explicitly when perturbing the conformation.

3.1. Determination of Structurally Rigid Regions

Due to the observation that chemical shifts capture local structural features (Mielke & Krishnan, 2009), they are used to define rigid and well-conformed regions (typically composed of secondary structural elements such as α-helices and β-sheets) that aid the fragment selection process (Shen et al., 2008). In contrast, disordered or flexible regions of a protein (for example, at the N- and C-termini) exhibit higher levels of structural diversity and therefore decrease the efficiency of sampling protocols. To overcome this bottleneck, CS-Rosetta makes use of TALOS-N (torsion angle likelihood obtained from shift and sequence similarity) (Shen & Bax, 2013), which takes as input the amino acid sequence and NMR chemical shifts to characterize the regions of rigidity and flexibility using random coil index (RCI). Using the predicted RCI values, the protein sequence is trimmed along the termini to eliminate any flexible end regions.

3.2. Fragment Picking

CS-Rosetta uses a fragment picker module built within Rosetta to construct high-quality sequence-specific fragment library using chemical shifts and a protein structure database called “vall” (shipped with both Rosetta software suite and CS-Rosetta toolbox) (Gront et al., 2011). The vall database consists of multiple sets (referred as chunks) of protein chains from X-ray structures without any chain-breaks. Each chunk is split into sets of fragments of required size with overlapping amino acid sequences in a “sliding window” fashion to cover the entire protein chain. The selection of polypeptide segments of the target protein from vall into the fragment library is governed heavily by the secondary structural elements predicted using TALOS-N (Shen & Bax, 2013). The plausible candidate fragments selected from vall are scored based on target sequence profile, chemical shift compatibility, secondary structure similarity, the presence of residues in the allowed regions of Ramachandran map, and RMSD of φ and ψ dihedral angles. In addition, the fragment picker allows flexibility to choose structures from which the segments are derived. For instance, homologous structures can be excluded during fragment library construction, which is typically desired while running benchmark calculations for method development. For structure calculation using CS-Rosetta, 200 fragments of 3 and 9 amino acids in length are picked for each residue position in the target sequence. The sampled fragments in these libraries define the degrees of freedom (or backbone torsion angles) a polypeptide chain can adopt during the structure calculation process. Rosetta applies these degrees of freedom to a target chain being modeled using the fragment insertion mover. The larger fragments (nine AA) are used during early exploration phases of structure calculation, whereas small fragments (three AA) are utilized predominantly in the later refinement stages (see Section 4).

4.1. Abrelax

Ab initio relax (or Abrelax) in CS-Rosetta is a standard protocol that uses as input the amino acid sequence and NMR chemical-shift-guided segments and assembles them using a fragment sampler to model near-native protein structures (Fig. 1, violet) (Shen et al., 2008). This protocol is invoked through a mover called the “AbrelaxMover,” which consists of an initialization step and an apply step. The initialization step creates instances of protocols used for idealization, fragment sampling, loop closure, and structure refinement. The apply step is utilized to invoke each of the protocols using the instances created during the initialization step. The AbrelaxMover is initiated via BrokerMain, a file that contains the main function used to start the protocol.

Fragment sampling is carried out in the low-resolution mode in four distinct stages (I, II, IIIa, IIIb, and IV) (Fig. 1, violet). During each stage, fragment insertion moves are performed on the Pose. This process involves selecting a candidate segment from the fragment library at a given position and adjusting the torsion angles of the target protein chain to match the angles from the segment. The resulting Pose is evaluated using stage-specific score terms and weights (Fig. 1). The low-scoring candidates at each stage are retained for further sampling. During stage I, 2000 fragment insertion moves with segments of AA length 9 are performed using Monte Carlo method with a score term to penalize any steric clashes introduced during the moves. In stage II, the same number of moves are performed with score terms, including the hydrophobic core, pairwise intramolecular interactions, and secondary structure packing to capture the communication between residue-pairs. In stage III, 10 iterations of 2000 fragment insertions are used to promote sampling of both local (IIIa) and nonlocal (IIlb) β-strand pairings with the help of bias used for scoring. In the last stage of this protocol, 4000 cycles of small fragment insertion moves are carried out, and the resulting models are scored using terms contributed heavily by radius of gyration to form closely packed structures. These low-resolution structures further undergo postfragment sampling stages that include (i) loop closure, (ii) conversion to full-atom representation, (iii) energy refinement using relax (full-atom refinement of structures in the Rosetta force field) (Tyka et al., 2011). For the assessment of the resulting high-resolution models, experimental chemical shifts are compared with predicted chemical shifts from the models using SPARTA (shifts predicted from analogy in residue type and torsion angle) (Shen & Bax, 2007). The deviation between two chemical shifts is used to adjust the energy function (Shen et al., 2008).

The Abrelax protocol provides tuning parameters to increment the iteration cycles, modify the weights of score terms that are applied during each stage, and skip any phase if deemed necessary. Setting up structure calculations using Abrelax is made easy with the help of automated setup scripts (see Section 9). These scripts are used to initialize method parameters to standard values (based on a set of benchmark calculations carried out previously; Raman et al., 2012; Shen et al., 2008). These parameters are defined in files known as flag files. Here, we provide a brief description of flag files used during Abrelax calculations.

flags_abrelax typically consists of (i) option to include H^α atoms during centroid mode (-chemical ::patch_selectors CENTROID_HA). This option is incorporated if input NOE or RDC restraints contain H^α atoms, (ii) path to the topology broker setup file (-broker::setup), (iii) path to chemical shift file for the purpose of evaluating output structures (-evaluation:: chemical_shifts) if specified by the user, and (iv) a parameter to describe the number of structures to generate (-out::nstruct).

flags_denovo is used to define (i) weight correction terms for scores contributed by radius of gyration, and residues present in α-helices and loops, (ii) loop closure options such as -overwrite_filter_scorefxn to force specified weights rather than predefined weights that are applied to filter models when closing loops, (iii) fold tree chain-break options (for example, -jumps::ramp_chainbreaks, which is used to tune the bias applied to chain-break energy terms across multiple stages of simulation), (iv) logging information in the event of errors, and (v) fully qualified paths to input fragment and amino acid sequence files.

flags_nmr_patches file includes information about patch files, which describe weights for specific energy terms to be used during different stages of ab initio (-abinitio::stageXX_patch, where XX represents stages of the Abrelax protocol), full-atom relax (-score::patch), and loop closure (for instance, -patch_filter_scorefxn used to filter models when closing loops).

nmr_relax_patch defines atom_pair_constrai nt and rdc weight terms used to bias restraint-based scores (see Section 6) provided typically by NOE and RDC measurements that are applied during structure modeling and refinement.

nmr_patch is similar to and contains same set of parameters as the nmr_relax_patch file. The primary difference is that they are applied during different stages of structure calculation. For ab initio stages and loop closure, nmr_patch is used. For full-atom relax, nmr_relax_patch is used. Their usage can be found within flags_nmr_patches.

setup_init.tpb is a topology broker file, which contains definition of claimers, such as ConstraintClaimer, that is applied to the Pose during centroid and full-atom stages. The ConstraintClaimer defines degrees of freedom through structural restraints. The claimer definitions specify fully qualified paths to the constraint files and weights that are applied during scoring.

After setting up the protocol for Abrelax calculations, a directory consisting of input, run, and test folders is created. The input folder contains chemical shift, constraint, and topology broker files. The test and run folders carry all the flag files (listed in flag_list.txt) that will be used during calculation. While the test folder is similar to the run folder, it is used mainly to check if the calculations proceed without any interruption or errors. Abrelax can be run using the executable, minirosetta.mpi. <ext> (where <ext> is platform and compiler dependent) in parallel on a multiprocessor system using the MPI (message passing interface) framework. Here, each CPU performs independent exploration and refinement steps. The MPI framework (i) enables seamless communication between processors independent of platforms used, (ii) manages real-time distribution of workload across computer clusters, and (iii) ensures safe resource management. During the execution of the protocol, calculated structural models are appended to “decoys.out,” a compact representation used to optimize storage space and is also referred to as silent file in Rosetta. As the protocol progresses, it writes information messages, warnings, and errors in log files, which can be used to troubleshoot in the event of failure or premature termination.

4.2. RASREC

RASREC is an iterative Abrelax protocol, which exploits the features of intermediate structures to build self-consistent core of a protein that resembles the native fold (Lange & Baker, 2012). To accomplish this task, RASREC employs several optimization strategies and embeds them into distinct phases consisting of exploration, feature reuse (or resampling), and structure refinement (Blum, Jordan, & Baker, 2010; Bradley & Baker, 2006; Brunette & Brock, 2005, 2008; Qian et al., 2007). A subset of the optimization strategies is dedicated toward sampling nonlocal β-sheets due to the complexity involved in modeling them (Bradley & Baker, 2006).

RASREC is invoked through the Archive framework (see Section 5) via BrokerMain and runs across six stages (Fig. 1, tan). During early stages of the protocol, secondary structures are used to determine β-strand pairings, which allow sampling of diverse β-sheet topologies (Bradley & Baker, 2006). The strand pairings are represented using a fold tree, where the pair edges are added into the tree through fragment insertion moves. For this purpose, segments of β-strands are obtained from a database of β-sheets from high-resolution X-ray structures. Later, the fragments (both β-strand fragments and those from sequence-specific fragment library) are assembled using the four stages of ab initio. The low-scoring models are recycled to produce new libraries of structural fragments for use in the subsequent stages of assembly (Blum et al., 2010). The folds of the low-resolution models in the later stages are guided by the trajectories of successful candidates from early stages (Brunette & Brock, 2005). Following these steps, loops are closed in the low-resolution models, which are converted to full-atom representation mode and refined in the Rosetta force field to produce physically realistic structures.

RASREC is a highly parallelized protocol that makes use of the MPI framework. Specifically, it utilizes a workload manager to distribute batches of structure calculation jobs across computer clusters during multiple stages of resampling and refinement. The resampling stages are divided into four centroid (I–IV) and two full-atom (V and VI) modes (Fig. 1, tan). In the first two stages, secondary structural elements obtained from NMR chemical shift-guided fragments of length three (see Section 3) are used to enumerate all possible β-strand pairings. All the consecutive residues assumed to form β-sheets referred to as β-residues are considered in pairs to account for all possible peptide plane orientations and pleating. The use of nonlocal β-residues in these stages introduces chain-breaks in the fold tree (Bradley & Baker, 2006), which allows sampling of diverse and nonlocal β-sheet topologies. From these stages, a set of batches and two centroid pool directories are created. Each batch directory consists of intermediate structural models computed using ab initio (stages I–IV) together with enumerated strand-strand pairing data and custom parameter flag files used for calculation. The stage-specific centroid pool directories contain low-scoring decoys from respective phases and a STATUS file presenting structure acceptance history into the pool. Stages II and III resample β-topologies from a set of decoys found in the centroid pool of previous stages. The β-strand pairings are scored using features such as contact order (residue-residue interactions that are close in space), sequence separation (between residues located within strands), and frequency of occurrence of the strand pairs in decoys. Recurring β-strand pairings with high contact order are favored. The score contributions of selected strand pairs are summed to give the total score for the β-topology. In stage IV, the protein fold is guided by a restricted conformational search. Specifically, the annealing trajectories of successful candidates are used to restart the structure calculation from stage II (Brunette & Brock, 2005). Besides sampling complex nonlocal β-folds, other structural features are resampled during stages III to VI from the segments of low-resolution conformations (from early stages) determined using fragment assembly. The centroid pool directories (from stages III and IV) are populated with structural models for the later stages. In stages V and VI, incomplete low-resolution models from previous stages (due to the use of broken-chain fold tree to sample nonlocal β-strands) are closed. The chain-break (or loop) closure is carried out using two strategies, (i) idealization and (ii) loop rebuilding.

Loop closure by idealization attempts to bridge the broken chain by selecting physically realistic values of bond lengths and angles for each atom pair in that region. However, the introduction of such values may result in radical changes in the fold tree. To avoid these changes, high weight restraints are added to govern the lever arm effects and to ensure that the perturbations are < 0.1 Å. Nonetheless, it is not always guaranteed that this procedure will not disturb the identified fold, in which case RASREC opts for a second strategy to close loops. This strategy involves rebuilding the residues of the loop by performing fragment assembly followed by loop closure using cyclic coordinate descent (CCD) (Canutescu & Dunbrack, 2003; Wang, Bradley, & Baker, 2007). The score terms applied in the course of this stage rely on the number of residues to be rebuilt. Whereas chain-break energy terms are used while building smaller subunits (< 10 residues), additional terms contributed by residue-pair interaction are applied for larger loops.

The RASREC method setup is similar to that of the Abrelax protocol, where a directory consisting of flag files is created. In addition to the parameters used by Abrelax, RASREC makes use of flag files listed later.

flags_iterative is used to define control flags that are applied during different stages of resampling. For instance, -iterative:fullatom enables the method to run the full-atom structure modeling and refinement stages (V, VI), -iterative:accept_ratio describes the cutoff value to be used for switching to next stage. Here, if the number of structures chosen as plausible candidates falls below the cutoff, then the stage is assumed to be complete.

super_quick_relax_patch.txt is similar to nmr_patch or nmr_relax_patch files that contain weight definitions for scoring restraints derived typically from NOE or RDC measurements. The primary use of this file is to provide constraint score weights for quick side chain repacking within individual batch files.

flags_fullatom defines the parameters that can be used during the process of (i) closing loops such as -loops:idealize_before_loop_close (idealizes the structure before attempting to close loops) and (ii) structure refinement using relax.

flags_cs_rescore contains path to chemical shift file used for model evaluation together with weights for chemical shift score term applied during both centroid and full-atom modes. Bias terms (such as -iterative:centroid_before_quickrelax_weight and -iterative:fullatom_after_ quickrelax_weight) are also included for the purpose of filtering models before adding them to the structure pool for subsequent stages.

nmr_pool_patch has weight terms applied to restraint-based scores and chain-breaks. These scores are used to filter models before adding them into the structure pool.

flags_rasrec is the flag file used in RASREC instead offlags_abrelax (used by Abrelax). This file includes reference to flags_cs_rescore and topology broker (setup_init.tpb) files.

Unlike Abrelax, RASREC creates a set of target-dependent batch directories alongside stage-specific directories at the end of each of the six phases. The structures from the centroid pool stages are added into centroid_ pool_stageXX directories, where XX = {1, 2, 3, 4}. Similarly, the full-atom decoys from stages V and VI are added to fullatom_pool_stageYY directories, where YY = {7,8} (even though they are stages V and VI, they are called 7 and 8). The final structural models are available in the form of a silent file (decoys.out) in fullatom_pool folder within the run directory.

To highlight the difference between Abrelax and RASREC, we calculated the structure of 20kDa Northeast Structural Genomics (NESG) consortium target, SgR145 (Wunderlich et al., 2004). Both the structure calculations were setup using NMR chemical shift-guided fragments and previously measured NOE restraints (Lange et al., 2012) (BMRB ID 16806). We sampled a total of 1000 structures using Abrelax, whereas RASREC runs produced nine batches consisting of 100 decoys each before convergence. We found that the lowest energy models generated using Abrelax and RASREC exhibit C^α RMSDs of 2.4 and 1.6 Å from the native structure (PDB ID 2N47; Tang et al., 2015), respectively (Fig. 2). Moreover, the accuracy of side chain dihedral placement (χ¹) of buried residues calculated using 10 Ų solvent accessible surface area cutoff with respect to the reference (PDB ID 2N47) is approximately 65% and 69% in the lowest energy Abrelax and RASREC models, respectively (Fig. 2A and B). The convergence of 10 lowest energy structures generated using RASREC is much higher than the Abrelax models (Abrelax: 43% and RASREC: 76%) (Fig. 2C). Most importantly, the structures sampled using RASREC are far superior compared to Abrelax models as indicated by the energy funnel, suggesting that the optimization strategies implemented in RASREC exhibit great promise in finding near-native structures rapidly (Fig. 2D).

4.3. AutoNOE

Alongside de novo structure prediction, CS-Rosetta offers an algorithm, AutoNOE, for automated assignment of 2D, 3D, or 4D NOESY cross-peaks (Lange, 2014). The NOESY cross-peaks result from coupling of two proton nuclei with different chemical shift values during magnetization transfer (Keeler, 2010). The process of labeling these cross-peaks in the NOESY spectra is known as NOE assignment. Whereas information about atom labels can be obtained from a known structure or an assigned resonance list, AutoNOE makes use of assigned chemical shift list to carry out the NOE assignment task. These assigned NOESY cross-peaks give rise to interproton distance restraints that are employed by structure determination methods to guide the conformational search. In CS-Rosetta, the NOE assignment (or NOESY) module is inspired largely by alternative automated tools for NOE assignment including ARIA (Linge, Habeck, Reiping, & Nilges, 2003; Nilges, 1993, 1997), AutoStruct (Huang et al., 2003), and CANDID (Herrmann, Güntert, & Wüthrich, 2002). AutoNOE takes as input the assigned NMR resonances together with unassigned NOESY cross-peak lists and generates a set of NOESY assignments from which it derives NOE distance restraints for structure calculation within RASREC.

The objective of the NOESY module is to assign proton atom labels from the input chemical shift list to the NOE cross-peaks such that the values of frequencies in both these lists match within the user specified tolerance levels (Fig. 1, green). To achieve this, the method starts by mapping input chemical shift frequencies to the unassigned NOESY cross-peaks. The initial mapping results in a set of ambiguously assigned NOEs due to potential overlap in the spectra (Nilges, 1993, 1997). These ambiguous NOEs are filtered based on a series of criteria. At first, the diagonal peaks (these are the peaks, which have same frequency for the proton and carbon dimensions) are eliminated. Second, a set of scores are computed to determine the confidence in the assignments. Third, NOE cross-peak assignments are filtered using the thresholds determined based on the confidence scores. Finally, NOE cross-peak intensities are calibrated to generate distance restraints.

5.1. Structure Modeling Methods

Abrelax can be accessed by the broker_main() procedure present in the BrokerMain source file (Fig. 4). This procedure creates an instance of the AbrelaxMover class and presents it to the MPI-based job distributor (or JD) that handles distribution of jobs in a parallelized manner and resides in the jd2 namespace. The JD is responsible for assigning new jobs to the worker nodes (or processors) and run independent instances of AbrelaxMover on the Pose. JD works well when writing the output structures one at a time into a single file (decoys.out) and examining it before processing new jobs. However, it cannot handle unexpected file read/write (race condition) behaviors from processors that may occur during parallelization. To avoid this situation, MPI framework is used. In particular, MPIFileBufJobDistributor (a derived class from JobDistributor) decouples IO (input/output) and job distribution (Lange & Baker, 2012). A detailed working of the MPI framework in the context of CS-Rosetta is given in the supplementary information of Lange and Baker (2012). In this discussion, we have abstracted MPI-based job distributors with their parent class, JobDistributor.

The main sampling strategy of CS-Rosetta is in the FragmentSampler class, which provides implementations of methods that carry out stages I–IV of Abrelax. ConstraintFragmentSampler inherits from FragmentSampler and is used primarily to associate input constraints with the Pose and further overrides functions from its parent such as do_stage1_cycles() and prepare_stageYY() where YY = {1,2} to carry out stage-specific operations. After the completion of sampling, the Pose is passed to the topology broker, which ensures that it (Pose) respects all the claims defined in the setup_init.tpb file. The Pose is returned to AbrelaxMover, which (i) closes chain-breaks introduced during sampling, (ii) switches from centroid to full-atom representation, and (iii) performs energy refinement to produce realistic structure that is written into decoys.out file by the process managing IO.

The RASREC protocol like Abrelax can be triggered from BrokerMain (Fig. 4). RASREC interacts extensively with the Archive framework, which lives in the archive namespace under jd2. Archive was developed primarily to facilitate interfaces from the iterative protocols (through their parent, ArchiveBase class) to the controller class, ArchiveManager, that communicates with the job distributor. In particular, ArchiveManager receives messages such as QUEUE_EMPTY and JOB_COMPLETION from the job distributor, based on which, it invokes necessary interfaces from ArchiveBase. When the protocol begins execution, the job distributor will examine the queue for any pending jobs (which is empty). Therefore, it sends QUEUE_EMPTY message to the ArchiveManager, which subsequently calls the generate_batch() method from the IterativeBase class, a descendent of ArchiveBase (through polymorphism). During the process of batch generation, batch_XX directory containing all the necessary parameter and protocol files is created, where XX is the batch number. ArchiveManager writes the status of the batch in a BATCH_INFO file within each batch directory, which is used to retrieve the status of a calculation if the protocol stops prematurely. Whenever the current batch is saturated, the job distributor sends JOB_COMPLETION message to the ArchiverManager, which forwards the newly generated decoys to the read_structures() method in IterativeBase. The read_structures() method filters and adds decoys to the structure pool. The main functionality of RASREC is implemented in IterativeBase, which calls procedures such as gen_enumerate_pairings() (stages I and II), gen_resample_topologies() (stages II and III), gen_ resample_stage2() (stage IV), and gen_resample_fragments() (stages IV–VI) to perform stage-specific operations. Here, gen_resample_topologies() procedure invokes PairingStatistics to compute scores for β-strandpairs during β-sheet topology resampling. IterativeBase is accessed using a wrapper class, IterativeAbrelax, which references IterativeCentroid (stages I–IV) and IterativeFullatom (stages V and VI) classes (both derive from IterativeBase) based on the representation mode of the protocol. The number of batch directories to be created is decided in real time, therefore, the job distributor, specifically MPIArchiveJobDistributor (a subclass of MPIFileBufJobDistributor) enables addition of batches to the job queue spontaneously. The decoys from batches are randomly selected and stored by HedgeArchive in centroid_pool_hedge directory. These archived structures are used during stage V of RASREC to promote extensive sampling of the full-atom energy function.

Evaluation of the decoys is carried out based on the energy terms weighted using custom or predefined bias by EvaluatedArchive, a derived class from ArchiveBase and NormalizedEvaluatedArchive, which inherits from EvaluatedArchive. The aggregation of weighted energy terms is represented by score_final tag for each structure in the silent file (decoys.out). In addition, the weighted sum of scores in the intermediate stages undergo on-demand updates and are made rapidly accessible under _archive_ select_score tag in the silent file carrying structures to subsequent stages (or the structure pool).

5.2. NOE Assignment Module

The NOE assignment method within the noesy_module namespace can be invoked from a publicly available application, r_noe_assign.<ext> (where <ext> is platform and compiler dependent) or from RASREC through the generate_batch() method of the IterativeBase class (Lange, 2014). The public app is called by the CS-Rosetta toolbox (see Sections 4 and 9) during initialization and finalization of NOE assignments, whereas the generate_batch() method calls the NOESY module at the start of each different stage of the RASREC protocol (Fig. 1, green). The assign() interface used by both these points of entry marks the start of the NOE assignment algorithm and is available in NoesyModule.impl.hh (Fig. 5). Here, methods within the CrossPeakList class are called specifically to carry out score calculations during the assignment of NOESY cross-peaks. AutoNOE is designed to accept multiple unassigned cross-peak lists. Each input list is encapsulated within the CrossPeakList class that provides functionality to read/write from iostream/files, update relevant scores, calibrate peak lists, and generate centroid and full-atom restraints for use in RASREC.

Each cross-peak in the list is associated with a set of properties, which are held by the CrossPeak class. These include proton spins mapped with chemical shift frequencies, peak intensity, cross-peak information (for each proton dimension such as, atom name and tolerance, which are stored as instances of the CrossPeakInfo class), and input resonance list. The CrossPeak class implements methods to find and assign atom labels to spins based on chemical shift frequencies and their tolerances. The Spin class contained within CrossPeak is used to encapsulate resonance frequency and spin assignment identifier which can be used to index into the resonance table. For the purposes of faster access to the assignments and residues (instead of referring from CrossPeak every time), PeakAssignment and PeakAssignmentResidueMap classes are used. PeakAssignment class stores score terms, spin identifiers, and resonance labels for each assignment of the cross-peak. PeakAssignmentResidueMap stores a mapping of the PeakAssignment list to the target residue list.

Calibration of peak intensities is handled separately by the PeakCalibrator class that exhibits simple association with CrossPeak. PeakCalibrator is inherited by StructureDependentCalibrator and StructureIndependent Calibrator which are invoked depending on the availability of decoys (this relies on the stage of the protocol) during calibration (see Section 4).

To store input chemical shifts and their atom labels, ResonanceList and Resonance classes are used. ResonanceList is the container class that keeps a list of resonances (each defined by Resonance class), identifiers for each resonance, and the corresponding residue number in the protein sequence. In particular, ResonanceList can be used to obtain tuples of chemical shift frequencies and respective atoms. Resonance encapsulates chemical shift id (or label; this is just an integer), frequency, tolerance, atom, residue, and atom-specific calibration flag (see Section 4.3). It is inherited by ProtonResonance, LabelResonance, and FloatingResonance. Depending on the type of atom encountered while reading chemical shift file, ProtonResonance (for protons) or LabelResonance (for other atoms associated with protons) is instantiated. In addition, ifthe given chemical shift frequencies lie in the folded region of the NMR spectra, then FoldResonance is utilized so that necessary corrections are applied while mapping the frequencies to the NOE cross-peaks.

6.1. NMR

NMR spectroscopy provides a means to obtain local and global structural information of proteins, which is later converted to distance or torsion angle restraints used to drive the structure modeling calculations. Local restraints are typically provided by NMR chemical shifts and NOEs (Anglister, Srivastava, & Naider, 2016; Güntert, 2009), whereas global information is obtained from RDC (Blackledge, 2005), PCS (pseudocontact shift) and PRE (paramagnetic relaxation enhancement) (Clore & Iwahara, 2009) restraints.

NMR chemical shifts are highly informative parameters that report on protein backbone and side chain conformations and are utilized heavily by NMR-based structure determination programs to model near-native structures (Berjanskii & Wishart, 2017). CS-Rosetta exploits NMR chemical shifts in two ways: first, to derive torsion angle restraints using TALOS-N (Shen & Bax, 2013), that are applied during the fragment selection process; second, to measure their consistency with SPARTA-predicted (Shen & Bax, 2007) chemical shift values obtained from simulated (intermediate or final) structural models. This fitness measure (deviation of experimental and predicted chemical shifts) given by chi-squared term is used to adjust the score function. During structure modeling in CS-Rosetta, the contribution of chemical shifts to the score function can be observed under chem_shifts column in the output file containing decoys. The corresponding weights for this term during centroid or full-atom modes can be customized using -cenpool_chemicalshift_weight and -fapool_ chemicalshift_weight tags present in the flag files.

Chemical shifts are generally complemented with long-range distance restraints provided by NOEs (distance between atom-pairs that are close in space) (Güntert, 2003). While unassigned NOE peak lists can be used by AutoNOE algorithm to perform NOE assignment and structure calculation iteratively (see Section 4.3), it requires higher levels of resonance assignment completeness (> 70% of resonances must be assigned) (Evangelidis et al., 2018). Therefore, when it is not feasible to obtain near-complete resonance assignments, structure determination algorithms (Abrelax and RASREC) can use all available chemical shifts together with NOE restraints assigned manually to model structures. NOE restraints are employed primarily to penalize the models that disagree with their distance bounds during structure calculation runs (Table 1). The amount of penalty applied is given by the atom_pair_constraint score term in the output file containing decoys. The weights used to compute atom_pair_constraint (represented by the same name in the patch files) can be tuned up or down as per custom requirements during different stages of calculations in centroid or full-atom modes. The NOE distance restraints typically start with the keywords, AmbiguousNMRDistance or AtomPair (or NamedAtomPair), and further contain definitions of atom-pairs, and penalty function type and its parameters (Table 1).

Whereas chemical shifts and NOEs are sufficient to model lower-molecular-weight complexes (<25kDa), additional data providing information about orientation of local and global regions is generally necessary to determine structures of large proteins or oligomeric systems. Such information is available from RDC measurements, which are useful in describing the orientation of the amino acids in proteins with respect to the global alignment frame (or the external magnetic field) (Donald & Martin, 2009). In CS-Rosetta, backbone RDC measurements are used in the score functions to evaluate their fitness to the structural models generated during centroid and full-atom stages (Raman et al., 2012). Corresponding scores in the decoys file are represented using rdc and rdc_filt (rdc_filt = rdc score term * rdc weight term) terms. Whereas a bias on RDC scores is indicated by the rdc term in the patch files.

The presence of paramagnetic metal ions in proteins induce significant changes in chemical shift values resulting in long-range structural restraints known as PCS (Schmitz, Vernon, Otting, Baker, & Huber, 2012). PCS measurements can be combined with Abrelax to model protein structures with increased accuracy (Schmitz et al., 2012; Yagi et al., 2013). To perform PCS-based structure calculations for a target with one or multiple paramagnetic tags, preliminary models are determined using only the chemical shift-derived fragments. These models are used to tune target-specific PCS score weights for each paramagnetic tag, which are applied during low-resolution mode of the structure calculation runs using -abinitio::stageXX_patch option where XX = {1, 2, 3a, 3b, 4}. Since the PCS measurements are recorded primarily for backbone atoms, they are applied only during low-resolution modeling phase. Following low-resolution phase, models are refined using standard refinement protocols within Rosetta. These models are filtered using a combined score produced by PCS data (during low-resolution stage) and full-atom refinement.

PRE measurements can be used to obtain distances between C^α atoms of PRE-affected residues and C^β atoms of the residues to which a paramagnetic probe is attached. While NOEs generally provide distance restraints between protons that are <6 Å apart in space, PREs can be used similarly as bounded restraints (Table 1) providing information about atoms that are spaced between ~20–25 Å (Tamaki et al., 2016). Moreover, the weights of these distance restraints can be tuned up or down using standard patch file containing atom_pair_constraint term. Alternatively, CS-Rosetta can be combined with solvent accessible data acquired through PREmeasurements (sPRE-CS-Rosetta) (Hartlmüller, Göbl, & Madl, 2016). Such data as demonstrated by Hartlmüller et al. (2016) is used within scoring functions to validate modeled structures during low-resolution and high-resolution modes. The experimental data file containing sPRE values is input in a format similar to the chemical shift file format in TALOS. Score weight represented by spre term in a patch file is added to the various ab initio stages given by -abinitio::stageXX_patch option in the flag files, where XX = {1, 2, 3a, 3b, 4}. Simultaneously, more parameters used in scoring, such as grid resolution (used during prediction of sPRE values from models), nonstereospecificity (using sPRE values obtained by averaging over stereospecific atoms) in full-atom mode, distances between (i) paramagnetic agent and the nuclei and (ii) the outer boundary of paramagnetic grid and the protein, can be adjusted for custom targets.

6.2. Homologous Proteins

Restraints derived from homologous or evolutionarily related proteins provide complementary information alongside sparse experimental data. CS-Rosetta offers a mechanism to generate such restraints and further supports data obtained from external tools or webservers to model accurate structures (Nerli et al., 2018).

CS-HM-Rosetta (or Chemical Shift-Homology Modeling-Rosetta) (Thompson et al., 2012) is a protocol whose primary purpose is to generate restraints from homologous proteins that can be combined with NMR chemical shifts within main structure modeling methods (Abrelax and RASREC). Here, the protocol first performs multiple sequence alignment of a query protein to all the proteins in the PDB using HHpred webserver (Söding, Biegert, & Lupas, 2005) or HHsearch (Söding, 2005) stand-alone software. The resulting alignment profile is used by CS-HM-Rosetta to measure distances between C^α atoms of residues that are consistent in multiple template proteins and are validated using criteria including, local and global sequence similarities and buried surface (Thompson & Baker, 2011). This validation information is combined to produce interatomic distance restraints as a mixture of Gaussian probability densities. The restraints are given higher weights if they are observed between same pair of atoms across multiple templates. The combined Gaussian probability densities are represented using a SOGFUNC function in restraint definition (Table 1).

More recently, sequence covariance information has been used to derive structural restraints (termed evolutionary coupling (EC) restraints) to guide the search for native conformation (Balakrishnan, Kamisetty, Carbonell, Lee, & Langmead, 2011; de Juan, Pazos, & Valencia, 2013; Jones, Buchan, Cozzetto, & Pontil, 2012; Kamisetty, Ovchinnikov, & Baker, 2013; Marks, Hopf, & Sander, 2012; Ovchinnikov et al., 2017; Sułkowska, Morcos, Weigt, Hwa, & Onuchic, 2012). These restraints are generated by utilizing genomic data and applying statistical methods such as, pseudolikelihood maximization (Ekeberg, Lövkvist, Lan, Weigt, & Aurell, 2013) and/or direct information (Marks et al., 2011) to measure evolutionarily coupled residues that are closer in space. Currently, there are many webservers and stand-alone software tools available to predict such restraints (Jones et al., 2012; Kamisetty et al., 2013; Marks et al., 2012). These distance restraints can be used alone or combined with sparse NMR data and CS-Rosetta algorithms to model high-resolution structures. Toward this end, CS-Rosetta expects data to be in a distance restraint file format (Table 1). For instance, if the restraints are generated using EVFold (Marks et al., 2011, 2012), then Rosetta provides scripts (see Section 9) to convert them into Rosetta compatible constraint files (with Sigmoid function; Table 1) that can be used directly within structure calculation runs as demonstrated previously (Braun, Leman, & Lange, 2015).