Forbidden space and feasible space.

We define the conformation space by the internal dofs of the protein, namely the torsion angles that are allowed to change throughout the motion pathway (see below). The conformation space is divided into forbidden and feasible regions (referred here as C-forbid and C-feasible, respectively; for illustration see Text S1). The forbidden regions correspond to all the conformations that involve high energy values, namely energy score above a threshold parameter, whereas the feasible regions comprise the low-energy conformations.

RRT algorithm with branch termination.

In algorithmic robotics literature, RRT is often biased by manipulation of node sampling, e.g. sampling in certain regions of interest. Here we take a different approach and rely instead on terminating branches that do not improve a certain predicate (note that branch termination has been previously used in a rather different context only, namely for avoiding imminent collisions under kinematic-dynamic constraints by sensing the immediate environment [45], and to deal with moving obstacles [42]). The input to the algorithm is an initial conformation and a set of partial information predicates (detailed pseudo-code and a full list of parameters are provided in Text S1). The tree is grown iteratively in small incremental moves to guarantee the smoothness of the motion. At each iteration, a new conformation, q_rand, is randomly sampled from the feasible space C-feasible. The nearest neighbor in the tree is then expanded towards q _{rand ,} by linear interpolation of the degrees of freedom from q_near to q_rand. Each path in the RRT tree can be considered a fine discretization of a continuous motion pathway in the feasible conformation space. The simulation terminates when: (i) the number of nodes in the tree is larger than N, a parameter for maximal tree size, or (ii) the tree could not be expanded for k consecutive iterations. The partial information predicates are used to choose a motion path that leads from the initial conformation to the conformation with the best predicate score.

The partial information bias is introduced by a filtering step. The filtering step is applied only in every other iteration, to allow an escape from local minima traps. In the filtering step, the branch that grows towards q_rand is terminated if it does not improve the partial information predicate after m consecutive interpolation steps (typically m = 2, again to allow an escape from local minima traps). The branch is terminated even if it leads to energetically feasible conformations (Figure 1b). The aim of this filtering step is to avoid expensive energy calculations in undesired directions. We note that existing branches are not pruned, only the growth of the current branch is terminated.

The effect of avoiding local minima with respect to partial information.

Figure 1 shows a toy example where branch termination is applied in all iterations. As desired, this narrows down the search to relevant regions of the conformational space (Figure 1b, grey circle). If we compare to Figure 1a, where branch termination is not employed, we see that the overall coverage of unexplored regions is compromised, but the target is reached faster (Figure 1c and 1d). By applying the global filtering step in every other iteration, we gain a bias towards partial information predicates, but still benefit from the rapid sampling of the unbiased RRT algorithm.

Side-chain optimization and local minimization.

In full atom-mode (see below), side-chains of generated conformations along the pathways are locally refined by the Rosetta Rotamer-Trial procedure to alleviate local steric clashes and optimize the interaction of side-chains with neighboring residues [63]. In addition it is advisable for full-atom runs to include short gradient-descent minimization with respect to all torsion degrees of freedom: very slight rotations of torsion angles (∼0.1°–0.2°) can alleviate local steric clashes and reduce the energy score substantially. To restrict the change in φ/ψ values, we added a heavy Gaussian penalty for deviating from the initial backbone torsions (Rosetta energy constraint CST_PHI_PSI with weight = 250.0). Note however that minimization is time consuming (3–8 times slowdown). Local minimization is optional and can be invoked by turning on a run-time flag.

Energy function: Full-atom vs. centroid mode.

We experiment with both the Rosetta full-atom energy function (Rosetta score12 [50]), which was shown useful for discerning native structures at atomic detail, and the coarser Rosetta centroid-mode energy function, a united-atom representation where side-chains are represented as centroid spheres (Rosetta score4 [49]). The latter allows rapid calculations at the expense of atomic detail, and has been used in a wide range of applications in Rosetta to speed up the conformational search by optimizing coarse features prior to atomic level optimization.

We assume here that the Rosetta energy function is relevant for our current task (see Discussion). The Rosetta energy function was optimized for native structures, but it includes physical van der Waals terms and solvation models, as well as a statistical hydrogen bonding term that was shown to correlate with quantum mechanical calculations [64]. However, we do note that PathRover is in principle not restricted to any specific energy function, and can be used in conjunction with other energy potentials as well.

Rosetta optimizations.

In addition to the full-atom and centroid-mode energy functions, the presented framework takes advantage of the elaborate infrastructure of Rosetta, including manipulation of molecular dofs, rapid side-chain optimization for fixed backbones (using the “rotamer trial” procedure described in [63]), energy minimization and caching of energy calculations. The framework can also use the gamut of other Rosetta features, such as closed loop sampling and sophisticated manipulation of backbone torsions, e.g. backbone fragment libraries or backrub motions [65]. These features are out of the scope of the current study, and will be explored in future work.

Running times.

All runs in this study were conducted on an AMD Opteron 275 2.2 Ghz/1 MB processor. Unless otherwise specified, in the Results section, centroid mode runs take the order of a few seconds to minutes each and Full-atom runs take 2–8 hours without energy minimization, and roughly 10–60 hours when the energy minimization flag was employed (see above), for growing trees of 30,000–100,000 conformations each. The number of dofs in different runs was between 8 and 198 backbone torsions (and all side-chains, see Table 2).

CesT domain swapping.

CesT is a type III secretion chaperone in Enteropathogenic E. coli that binds numerous effector proteins. In CesT the neighboring chains within the crystal lattice are domain swapped [66] (Figure 2a). The N-terminal domain (residues 1–33) and C-terminal domain (residues 38–134) from neighboring chains pack to form a monomer-like globular unit, the “pseudo-monomer” (Figure 2b). The pseudo-monomer can be well aligned to monomers in the other homologues (Figure S1), suggesting that there is a monomeric form of CesT that resembles the pseudo-monomer. Packing of non-swapped monomers against each other is mostly identical to their packing in the pseudo-monomers, as the sequence of domains is identical [56]. It is however not known whether the swapped conformation of CesT is a crystallographic artifact or whether it is the physiologically active peptide-binding form [66], and it is interesting to examine the possibility for domain-swapping motion of this protein.

Using the pseudo-monomer as a partial information predicate.

We first examine whether we can model a hinge motion in which a chain of CesT moves to the pseudo-monomer conformation where its two domains are interacting. We start from the swapped conformation (where the domains are farther apart), using chain A in pdb-id 1k3e [66]. We allowed backbone mobility in residues 34–37, a loop region that separates the two domains and is also predicted to be a hinge region by normal mode analysis (Figure 2a, spheres; see also Table 2 for additional evidence that these residues form a hinge). The algorithm was biased to minimize the RMSD between the initial conformation and the pseudo-monomer conformation, using the residues in the interface between the N-terminal domain and the C-terminal domain (Table 3). We applied RRT with branch termination, as described in Methods, and simulated the motion in both centroid mode and full-atom mode (where all side chain atoms are included, see Methods). Each centroid mode run was repeated 50 times, taking a few seconds only to complete. Full-atom mode simulations were performed with the energy minimization flag turned on, to relieve local steric clashes. Each such simulation took a few hours on a single processor, and was repeated 15 times.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Predicates used for guiding the domain swapping motion of CesT.

https://doi.org/10.1371/journal.pcbi.1000295.t003

We analyzed the runs that best minimized the predicate. Both in centroid mode and full-atom mode, a collision-free path towards the pseudo-monomer conformation was found. The initial conformation deviates by 16 Å RMSD from the pseudo-monomer, and the final conformation deviates by 0.8 Å RMSD in centroid mode, and 1.34 Å in full-atom (0.76 Å for a partial alignment without residues N24-I33). These runs provide a proof of concept that the biased RRT algorithm successfully employs bias for guiding the motion. The suggested motion is shown in Video S1 (full-atom mode) and Video S2 (centroid side-chains mode).

Biasing the motion with SigE, a homologue of CesT.

SigE (pdb-id 1k3s [66]) is one of a few distant homologues of CesT. The N-terminal and C-terminal domains of SigE are similar to those of CesT (RMSD of 1.7 Å and 2.5 Å respectively), and the pseudo-monomer of CesT can be aligned to the SigE monomer (Figure S1). However, the two structures share a very low sequence identity (18%), and the SigE monomer deviates by 3.96 Å from the pseudo-monomer (using the sequence alignment from [66]).

In order to investigate whether this distant homologue can indeed guide the motion towards the correct conformation, we devised a set of geometric predicates that use SigE as a reference for guiding the motion of CesT, such as the distance between specific atoms or the orientation between specific secondary structures (see Table 3 and Figure 3). For each predicate, we conducted 50 independent runs and analyzed the run that best minimized the predicate. We worked in centroid (united-atom) mode, as this example mainly serves to illustrate the effect of various predicates on the simulations. Running time was a few seconds for each simulation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Use of partial information in simulations of CesT domain swapping (in centroid mode representation).

(A) The final conformation along the motion pathway of CesT (cyan) is shown for five different examples of predicates (see Table 3). We show the best scoring run with respect to the specified predicate (out of 50 independent runs). The orientation of the N-terminal domain of the native pseudo-monomer is shown in red for comparison. (B) The RMS distance of the final structure, for simulations with different predicates, is plotted relative to SigE (the homologue that was used to guide the motion, blue) and the native pseudo-monomer of CesT (black). Even though the homologue was the reference for biasing the motion, the simulations reached the correct conformation with a better level of accuracy for several predicates. (C) Biasing the motion by combining several distance constraints (see Table 3 for details about the constraints): the results are shown as the fraction among 50 independent simulations that reached given RMSD thresholds (to the native pseudo-monomer).

https://doi.org/10.1371/journal.pcbi.1000295.g003

In Figure 3a and 3b we examine, for 5 of the predicates in Table 3, the RMSD distance of the final conformation from (1) SigE – the reference homologue protein, and (2) the native CesT pseudo-monomer (Figure 2b). Remarkably, in four cases, the final structure was more similar to the pseudo-monomer than to SigE, even though the reference for guiding the motion was SigE (the SigE monomer deviates by 3.96 Å from the pseudo-monomer). This suggests that the energy function prevented over-bias of CesT motion towards the structural features of SigE. This fact is particularly surprising since the simulations were conducted in centroid mode, without the atomic details of the side-chains, and demonstrates the effectiveness of a simplified and rapid energy function in this case. At such level of predictive precision (<2 Å), many side-chains can be already modeled quite accurately. Although our aim is motion prediction and not homology modeling, it is promising that the near-native conformation is recovered using very simple predicates.

Relative orientation of secondary structures.

In Sig-E and in the CesT pseudo-monomer (but not in the initial conformation), α-helices H1 and H3 each lie in a different domain of the protein (Figure 2b). In addition, β-strand B0 of the N-terminal domain is paired to β-strand B1 of the C-terminal domain (Figure 2b). In Table 3, we list the set of predicates that we formulated over the relative orientation of these secondary structures in the two domains. An interesting predicate is the Helix Line-Fit predicate, which combines three measures. We used least mean square line fitting (LMSLs, Table 1) to approximate the main axes of α-helices H1 and H3. We then measured the distance and angle between the fitted lines, as well as the distances between the centers of mass of each helix. This is a useful measure when relying on a homologue protein, since it is much less sensitive to alignment shifts that are characteristic for α-helices. In Figure 3b we observe that the final conformation was very close to the native pseudo-monomer (1.87 Å) but not to SigE, which is the reference for the biasing predicate. Looking closely at this example (Table S1), we saw that the line angle and line distance predicates were perfectly matched to their values in SigE, whereas the center of mass distance between the helices did not reach its value in SigE (15 Å) but rather reached its value in the CesT pseudo-monomer. Hence, the simulation was not over-biased by partial information constraints, and took into account the specific features of the simulated molecule.

β-sheet formation.

The combination of helix and sheet predicates (Table 3 and Figure 3) was sufficient to direct the motion to the native conformation. How does β-sheet formation affect α-helix orientation and vice versa? We observed that when the motion was guided by partial information on the orientation of the α-helices alone, the β-strands B0 and B1 still came close together and the final conformation exhibited similar structure to the native pseudo-monomer. In contrast, the helices did not move to the correct orientation when the only partial information provided was β-sheet formation. We note however, that this may also be an artifact of the smaller number of atoms involved in the strand-pairing predicate.

Combinations of atomic distance constraints.

Not surprisingly, a constraint on the distance between a single pair of atoms is not deterministic enough for guiding the motion (Figure 3a and 3b–left bar). The structural alignment between the final conformation of CesT and the pseudo-monomer is rather poor. It is clear that the distance between a single pair of atoms should be combined with other partial information or atomic detail constraints, in order to derive a more reliable target conformation and motion pathway. Therefore, we have examined what combination of distance constraints suffices for biasing the motion. Combinations of two or three distance constraints (Table 3) were used to guide the motion. In Figure 3c, we plot the percentage of 50 independent simulations that reached the native pseudo-monomer conformation up to various degrees of similarity (in RMSD). We observe that 2 or 3 constraints are still not enough to guide the motion in all simulations, but they lead to a much higher percentage of runs that reach the native conformation. This suggests that the combination of just a few distance constraints is an effective way of constraining motion-planning simulations.

Ribose binding protein (RBP): Ligand-binding-induced hinge movement: Incorporating loop closure constraints with simple predicates.

The problem of fixing remote structural segments that are connected by a flexible loop is known in the literature as the protein loop closure problem [67]. It might require complex loop closure calculations or interpolation of internal coordinates motion from Normal-Mode Analysis [38],[68]. Previous attempts have been made for ad-hoc solutions to this problem during RRT simulations [68],[69], as well as in the broader context of structural modeling [67],[70],[71].

The Ribose Binding Protein (RBP) belongs to a family of ligand-binding proteins that comprise two domains, connected by a hinge. Upon binding of the ligand in a cleft between the two domains, the domains approach each other to close the cleft (Figure 4a). However, unlike CesT, in RBP each domain is discontinuous with respect to the sequence. The hinge that connects the two domains is made of three separate stretches of sequence (Figure 4b and Table 2). Consequently, the hinge torsion angles must change in a coordinated way, to prevent the two domains to disintegrate. Although violating the integrity of domains is energetically unfavorable, a lot of running time may be consumed on sampling non-favored conformations that disintegrate the domains. Indeed, when we simulated the motion of RBP without any external constraints, the domains wobbled and partially disintegrated during the motion, with high energy fluctuations (results not shown). Although this type of motion cannot be negated completely, domain disintegration during relatively fast substrate binding motion contrasts basic biological intuition.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Ribose Binding Protein (RBP) architecture.

Domain A (cyan) is connected to Domain B (blue) by a hinge (red). (A) The RBP structure in its open and closed form, pdb-ids 1urp [83] and 2dri [84] respectively. (B) The architecture of RBP–each domain consists of discontinuous segments of residues. The two domains are connected by three hinges that must move in a coordinated way to maintain domain integrity. Domain boundaries are rough estimates. (C) The sequence of RBP showing the discontinuous domains and the secondary structures. This illustration was taken from 1urp Protein Data Bank entry at http://www.rcsb.org/pdb/ [5]; and domain assignments are from [85].

https://doi.org/10.1371/journal.pcbi.1000295.g004

In order to enforce loop closure, we have simply added a partial information predicate that penalizes disintegration of the domains (in terms of RMSD to the native domain). In the resulting pathway (Video S3), the two domains are kept in one piece throughout the motion, and only a small β-strand at the C-terminus of the protein (residues 266–269) deforms during the motion. The simulation of the motion in centroid mode is performed within a few minutes time. This demonstrates the flexibility of the partial information framework to efficiently address diverse settings, without the need for explicit ad-hoc calculations.

The central hinge: Allowing rotation in residues 48–55.

It has been suggested that residues 48–55 between the two repeat domains of Cyanovirin-N form a hinge region for domain swapping [72],[73]. Additional support comes from structural conservation patterns, difference in torsion angle values between alternative structures and Gaussian Network Models for detecting hinge regions (Table 2). We refer to this region as the central hinge of Cyanovirin-N.

Using the central hinge set of dofs (Table 2), we first used the simplified centroid mode representation to generate a low-energy motion pathway within minutes. Considering the experimentally determined high energy barrier for this motion, it is rather surprising that such a pathway could indeed be easily created. We thought this might be an artifact of the simplified representation of the structure in centroid mode: the barrier might be apparent only at a higher resolution level. We therefore proceeded to a full-atom representation:, When all side-chains atoms and hydrogen atoms were modeled explicitly, it was impossible to unlock the intertwined monomer, unless the energy threshold was substantially raised to 10⁵ Rosetta Score-12 units (which allows for extreme steric clashes). The domains did not unpack even in a long run of RRT, consisting of 300,000 conformations and taking a few days to run. The protein moved by no more than 1.5 Å from the initial conformation, over 13 Å away from the swapped target conformation. Video S4 demonstrates how the side-chains of one domain are tightly locked within the other domain, and the motion is confined within a steric “cage”.

The effect of local energy minimization.

Could very slight “breathing” motions of other degrees of freedom allow the domain-swapping of the protein? Local energy minimization involves slight changes (∼0.1–0.2) in all backbone and side-chain degrees of freedom and as such might suffice to alleviate local steric clashes of the sensitive full atom energy function (see Methods section for details). Indeed, in a new simulation with freedom of motion in the central hinge together with local energy minimization (“Central-Hinge (M)” in Table 2) we were able to generate a clash-free motion pathway in full-atom. Local minimization slowed down the rate of generating new nodes (3–8 fold), but allowed the initial unpacking of the domains after less than an hour run. It is striking to observe how minute structural breathing in all degrees of freedom can alleviate steric clashes and allow motion that was not possible when only the central hinge is mobile.

The effect of adding secondary hinges.

We now examine if the introduction of several additional dofs may provide the simulation with sufficient “breathing” flexibility to allow the large scale motion of the hinge, even without energy minimization. We inspected the structure and located two symmetry-related loops at residues 36–40 and at residues 87–91 that connect the two β-sheets in each domain. These loops appear as “weak links” in the protein chain between the two sheets (see Table 2 for more reasoning behind this choice). The addition of flexibility in the two loop regions (“Secondary-Hinges” in Table 2) allowed our simulations to find low-energy clash-free motions in full-atom mode (Video S5), without energy minimization in all dofs.

Analyzing hinge residues by restricted sampling of all DOFs.

We saw that restricted local energy minimization of all dofs, as well as the introduction of secondary hinges both enabled the domain swapping motion. In order to analyze the motion of all residues, we conducted another simulation where the central hinge (residues 48–55) is free to move, and all other dofs can also rotate by up to 30° from their initial value, allowing for more extensive “breathing” motion in all dofs (“Partially-Restricted (M)” in Table 2). The total number of dofs in this simulation is 198, with 16 dofs that are completely free to move. Using this large number of dofs (and including local energy minimization), the simulations took 3–4 days. A movie of such a simulation (see Video S6) shows the breathing motion of all dofs, and in addition suggests that side-chains L1 and W49 (marked in red) act as “gate-keepers” that interfere in the unpacking of the two domains. It would be interesting to examine the role of these residues experimentally and in-silico, although this is out of the scope of this work.

We should note that the S–S bonds between two adjacent β-strands, from C8 to C22, and from C58 to C73, were not modeled in the simulation due to technical limitations. However, we note that both of these bonds connect adjacent β-strands, and atomic distances between these pairs of residues are close to constant during all the simulation, suggesting that S–S bonds will not play a critical role.

Consistency of the hinge regions between runs.

In order to examine the residues involved in the motion, we have scored each residue for how often it serves as a hinge during the motion, at different resolutions of motion (see Methods and Figure S2). Low resolution hinges are involved in strong hinge motions, and high-resolution hinges are involved in milder ones. We conducted three independent simulations, and compared the consistency of the detected hinges at different resolutions of hinge motion (parameter ρ). We used Pearson's linear correlation (with values ranging between 1 for full-linear correlation, and 0 for no correlation). The correlation over all residues is very high and rises with decreasing resolution, such that the most prominent hinge motions are consistent between runs (Figure 5b, blue line; correlations range from 0.77 at ρ = 0.5 Å, to 0.96 for ρ = 4 Å). Since the central hinge is inherently biased by the run parameters, we also analyzed the correlation when excluding the central hinge region (green line). In this case the correlation is lower but still significant. As expected, the correlation is low at the highest resolution (ρ = 0.5 Å), where small, flickering, movements are measured. For a detailed inspection, we plot in Figure 5c the hinge scores in the three simulations at different resolutions (1.5 Å and 3.5 Å), and also show a Cyanovirin-N structure in cartoon representation colored based on the residue hinge score. Finally, a plot of the weighted average of hinge scores at different resolutions is shown in Figure 5d. Low-resolution hinges are involved in larger hinge motions and are assigned a higher weight, so that pronounced hinge regions are inspected:

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Hinge regions in independent simulations of Cyanovirin-N domain swapping.

All runs use the Partially-Restricted (M) set of dofs (Table 2), where the central hinge is allowed free motion, and all other residues can rotate by ±30°. (A) In each run, each residue is scored by how often it tends to be in hinge regions. Hinges are extracted by structural comparison between conformations along the motion sequence. They connect sub-domains that remain rigid during the simulation. Resolution parameter ρ states the RMSD threshold used for rigid alignments. Mild hinges appear only at higher resolutions (low value of ρ), and salient hinges appear at low resolutions (see Figure S2 for a detailed protocol). (B) Pearson's Correlation between the hinge scores of three independent simulations, for different values of ρ. In blue, the correlation over all residues, including the central hinge (residues 48–55). In green, the correlation when excluding the central hinge. (C) Hinge scores for each residue in three independent simulations, for ρ = 1.5 Å and ρ = 3.5 Å. The y-axis denotes how often each residue appears in hinge regions (see Figure S2 for more details). Secondary structures (according to DSSP [86]) are plotted along the x-axis. Observe that milder hinges disappear at lower resolution (3.5 Å). R and R* are the average Pearson's correlations between runs with and without the central hinge region, respectively. For each plot, the crystal structure of Cyanovirin-N is colored according to the corresponding hinge score, with warm colors indicating higher scores. (D) The weighted average of the hinge scores for different values of ρ (see Results). Since higher resolutions contain milder hinges, they were assigned a lower weight.

https://doi.org/10.1371/journal.pcbi.1000295.g005

Not surprisingly, a large peak at the central hinge (marked in dashed lines in the plots) dominates in all figures. Also interesting are hinges in other regions, where restricted motion of ±30° was allowed. Flexibility in the tail region is apparent, although quite trivial. Interestingly, the hinges do not appear in random location, but rather are consistent between independent runs. For instance, residues 27–30 form an apparent hinge at resolution ρ = 1.5 Å, for all independent runs. At resolution of ρ = 3.5 Å the consistency seems less remarkable at first, but a focus on residues 10–50 shows very similar hinge patterns between runs (Figure 5c, inset).

Energy analysis.

What is the contribution of different energy terms to the motion? For analyzing the estimated energy landscape of the motion, we used the Rosetta energy score with dampened van der Waals potential [75], to reduce the dominant effect of fluctuative contributions of slight steric clashes (in the simulations themselves we use the classic van der Waals potential, so any steric clashes are heavily penalized). In Figure 6a we observe the energetic barrier that results from unpacking of the two domains, both disrupting favorable attractive forces (vdW–blue, hydrogen bonds - green) and causing increased repulsion due to the motion in a cluttered environment (repulsive force, red line). The solvation term decreases as polar residues are exposed (the hydrophobic effect is evident in the much higher increase in the attractive vdW force).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. The contribution of different energy terms in the domain-swapping simulation of Cyanovirin-N.

We used the Rosetta soft-repulsive energy score for this analysis, in order to dampen repulsive fluctuations that are due to mild sterical clashes (see text, note that the simulation itself was conducted with the Rosetta score-12 energy function, with a classical vdW potential). Note that the y-axis shows the total energy score, and the specific energy terms are shifted by a constant number of units, for convenient comparison with the other terms. (A) Energy plot for simulation with restricted dofs (central hinge is free to move, and all other dofs can fluctuate by ±30°). (B) Energy comparison for the aligned sections of the restricted (left) and free (right) simulations. In the free simulation, all dofs are free to move. The reaction coordinates of the two simulation were aligned by a string matching algorithm [39] based on structural similarity between conformations (see text, Figure S3 and Video S8).

https://doi.org/10.1371/journal.pcbi.1000295.g006

The repulsive forces subside after the domains have separated (steps 30–69 on the x-axis). Now that the two domains are unlocked, they may be free to sample many conformations without significant clashes. Indeed, there are several alternative conformations of domain swapped Cyanovirin-N [73].

Comparison of restricted and free simulations.

Is the biasing for the central hinge indeed justified? In order to answer this question, we performed a simulation where all residues (except for the first and the last residue) were completely flexible and no bias was introduced (“Free” in Table 2). For 198 free dofs, the simulation took 5–6 days to generate a tree of 100,000 conformations. For the free simulation, the domains unpacked from each other substantially, but did not manage to reach the target conformation within the limitations of the running time (Video S7), probably due to the large number of non-restricted dofs.

We compared the free simulations to the “Partially-restricted” simulations described above and in Table 2. Each simulation results in a motion pathway that comprises a sequence of conformations. We aligned the simulated motion pathways based on RMSD between corresponding conformations (using a path alignment scheme where corresponding frames in the two simulations are aligned by a string matching algorithm, similar to sequence alignment methods, see [39]). The movie of the aligned motion of the restricted and free simulations (Video S8) demonstrates that the two simulations are very similar, and Cα RMSD between aligned conformations stays within 2–3 Å throughout most of the motion, growing to 3.5–4 Å only towards the end (Figure S3). Remarkably, comparison of the hinge scores in the restricted and free simulation (Figure 7), shows that the central hinge is the most prominent hinge at low resolution (left panel, ρ = 4 Å), which means it is involved in the largest scale motion in the free simulation. Milder hinge motions at resolution ρ = 1.5 Å are less correlated, although the hinges at residues 27–30 are still markedly observed in both simulations. Note that since the free simulation spanned a part of the domain swapping motion, the alignment is partial, comprising half of the partially-restricted simulation, and the entire free simulation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. Robustness to restriction of dofs in Cyanovirin-N simulations.

The hinge score for each residue is plotted for (A) restricted simulations, where only the central hinge (residues 48–55) is free to fully rotate, and the other residues are restricted to ±30° deviation from the initial conformation, and (B) a free simulation (magenta) where all backbone degrees of freedom are free to rotate (see Table 2). Hinge scores are plotted for resolution parameter values ρ = 4 Å and ρ = 1.5 Å. The central hinge is the most salient feature in the free simulation, and therefore it appears even in low-resolution plots. Milder hinges are less robust to the restriction of dofs (see text for more details).

https://doi.org/10.1371/journal.pcbi.1000295.g007

In Figure 6b we also compare the energy of the aligned simulations. In both simulations we observe the energetic barrier that results from unpacking of the two domains, both disrupting favorable attractive forces (blue) and causing increased repulsion due to the motion in a cluttered environment (repulsive force, red line). In both simulations the repulsive forces subside as the domains continue to separate and the attractive forces increase.

It is worth noting that motion-planning with such a large number of dofs (198) is not a trivial task, and that both simulations converged only when we used local energy minimization. Although energy minimization may often increase the running time, it may allow to deal with a larger number of flexible dofs, since unlikely random perturbations that lead to clashes can be countered by local minimization.

Experimental validation and analysis of simulations.

One of the big challenges of computational biology is the interface between computational and experimental observations. While full-atom experimental motion pathways of high resolution are still not in sight, significant progress has been recently made in experimental research of transient conformations. Distance constraints from FRET experiments, Paramagnetic-Resonance Enhancement, Residual Dipolar Coupling and other spectroscopic methods for assessing molecular dynamics can be used for (1) constraining simulations of molecular motion using measurable constraints, and for (2) validating motion pathways of suggested simulations, by comparing the measured distance constraints to simulated predictions. Our vision is that innovative experimental measurements of limited scope can focus and enhance computational techniques, effectively allowing researchers to generate realistic motion pathways that incorporate as much external information as possible within the currently suggested framework of PathRover. Particularly, this can allow for the design of experiments that target specific states within a motion pathway based on in silico predictions of large-scale motions. The predicted motions can be also used to suggest mutations, such as our suggested mutations in L1 and W49 for Cyanovirin-N (see residues marked in red in Videos S6 and S7).

Computational observations are most meaningful when they are well-defined in a way that poses them as clear, lab-testable hypotheses. To that end, it is not sufficient to rely on raw simulations results. In our work, we have therefore devoted significant effort for developing analysis and visualization tools for extracting physical features from simulated motion, including the protocol for analyzing hinge regions in simulated motion, as well as the visualization and space-time alignment of multiple motion paths. We believe that developing novel analysis and visualization tools is an important direction of future research, which is just as important as the simulations themselves, as it can provide the missing link between experimental and computational observations.

Applications to other types of molecular motion.

In this work PathRover was applied to motions of domain swapping and substrate binding. However, different types of molecular motions might have different characteristics with respect to the number of torsion angles that are involved in the motion, the scale of the motion, the role of side-chains, etc. One challenging class of molecular motions involves allosteric protein motions [80]. In this case, a large number of torsion angles are often involved in the motion, but each of them changes by rather small increments, and partial information might constrain the overall nature of the motion. Another interesting type of motion involves more than one molecule, such as docking of a protein or a flexible peptide onto another protein. Motion-planning techniques have been used for small-molecule docking [40], but to the best of our knowledge not for docking of two globular proteins or for protein-peptide docking. Of particular interest within this framework are cases where partial information can provide details about the approximate location of the interface, and conformational backbone flexibility of the monomer needs to be modeled efficiently [71],[81].

Analyzing multiple motion pathways.

One of the advantages of RRT-based techniques is their relative speed. A large body of motion pathways can be created at atomic level that includes side-chain atom positions. A large number of pathways provide further insights about the connectivity of the conformational space under a wide range of settings. In contrast, it is difficult to generate a large number of pathways using, e.g., MD simulations, due to slower running times. In Video S8 we showed an alignment between two motion pathways. We recently proposed a method to compare, cluster and merge multiple motion pathways from independent runs of the RRT algorithm. The merged pathways have lower energy or shorter length than all input pathways [39]. It would be interesting to examine the clusters of pathways that are generated with different types of partial information.