Research Article

Atomic-Level Characterization of the Structural Dynamics of Proteins

See allHide authors and affiliations

Science  15 Oct 2010:
Vol. 330, Issue 6002, pp. 341-346
DOI: 10.1126/science.1187409


Molecular dynamics (MD) simulations are widely used to study protein motions at an atomic level of detail, but they have been limited to time scales shorter than those of many biologically critical conformational changes. We examined two fundamental processes in protein dynamics—protein folding and conformational change within the folded state—by means of extremely long all-atom MD simulations conducted on a special-purpose machine. Equilibrium simulations of a WW protein domain captured multiple folding and unfolding events that consistently follow a well-defined folding pathway; separate simulations of the protein’s constituent substructures shed light on possible determinants of this pathway. A 1-millisecond simulation of the folded protein BPTI reveals a small number of structurally distinct conformational states whose reversible interconversion is slower than local relaxations within those states by a factor of more than 1000.

Many biological processes involve functionally important changes in the three-dimensional structures of proteins. Conformational changes associated with protein folding (1), signal transduction (2), the catalytic cycles of enzymes (3), and the operation of molecular machines and motor proteins (4) often involve transitions among two or more structurally distinct states. These states are often characterized as “basins” separated by barriers on an “energy landscape” (5).

Substantial progress has been made, using both experimental (1, 6) and computational (7, 8) techniques, in characterizing conformational basins and the ways that proteins move within and among them. It has proven difficult, however, to structurally characterize sparsely populated or disordered states and to elucidate the “basin-hopping” mechanisms involved in the interconversion of various states.

All-atom molecular dynamics (MD) simulations are designed to provide a high-resolution view of the motions of biological macromolecules (9), producing continuous trajectories with the potential to connect static structural snapshots generated from experimental data. Computational constraints, however, have limited such simulations to ~1 μs of simulated biological time. [The longest previously published all-atom MD simulation of a protein, at 10 μs, required more than 3 months of supercomputer time (10).] This has limited the usefulness of MD, as many biological processes involve conformational changes that take place on time scales between 10 μs and 1 ms.

To access such time scales, we designed and constructed a special-purpose machine, called Anton (11), that greatly accelerates the execution of such simulations, producing continuous trajectories as much as 1 ms in length. This has allowed new insight into two fundamental processes in protein dynamics: protein folding and the interconversion among distinct structural states of a folded protein.

Specifically, we have been able to formulate a detailed description of the folding of a WW domain (12) as well as the folded-state dynamics of bovine pancreatic trypsin inhibitor (BPTI), a workhorse in the study of protein dynamics [and the subject of the first protein MD simulation (13) and of pioneering computational studies of protein folding (14)]. Our choice of biological systems was motivated in part by the fact that a considerable amount of experimental data is available for both of these proteins, providing us with various ways to test the reliability of our simulations.

Folding of a WW domain. WW domains are small, independently folding protein domains that bind to proline-rich sequences. The topology of WW domains is characterized by two β hairpins, which form a three-stranded β sheet (15). Mutational analyses of the folding of WW domains show that the rate-limiting step in the folding reaction involves the formation of the first hairpin (1618). This information facilitated the original design of the fastest-folding WW domain reported to date, FiP35 (12), which folds in 14 μs.

FiP35 has several features that make it attractive as a model system for use in computational studies of protein folding. A great deal of experimental data is available for this system (12, 16, 17), and previous attempts to characterize its folding mechanism through explicit-solvent simulations have been largely unsuccessful (10, 19). It has been suggested that WW domains such as FiP35 fold through multiple distinct routes that differ in the order of formation of the individual hairpins (2022), but this has not been conclusively demonstrated. It has also been speculated that the mutations involved in the design of FiP35 may have shifted the rate-limiting step relative to that of the Pin1 WW domain, which formed the basis for its design. Finally, it has been suggested (12) that the fast folding of FiP35 may be close to the downhill regime, with only a small (<3kBT) free-energy barrier. Many of the unresolved questions surrounding the folding of FiP35 raise issues central to the process of protein folding more generally, and we believed they might be amenable to investigation using very long atomistic simulations.

Folding FiP35 and villin to experimental resolution using the same force field. We first ran a simulation at 337 K, a temperature at which FiP35 should be predominantly folded. The protein was initially configured in a fully extended state and was observed to fold to a stable conformation with a backbone root-mean-squared deviation (RMSD) of ~1 Å from the crystal structure (15) (Fig. 1B). A previous attempt to fold FiP35 computationally, using a 10-μs explicit solvent simulation (10), did not converge to the native state, and subsequent work (19) provided evidence that this was attributable to deficiencies in the force field. Our successful folding simulation was based on a modified version of the Amber ff99SB force field (23).

Fig. 1

Folding proteins at x-ray resolution, showing comparison of x-ray structures (blue) (15, 24) and last frame of MD simulation (red): (A) simulation of villin at 300 K, (B) simulation of FiP35 at 337 K. Simulations were initiated from completely extended structures. Villin and FiP35 folded to their native states after 68 μs and 38 μs, respectively, and simulations were continued for an additional 20 μs after the folding event to verify the stability of the native fold.

One potential concern in simulating the folding of an all-β protein is the possibility that the force field used might tend to overstabilize β sheets, thus “assisting” the folding process in a nonphysical manner. The force field we used, however, also folded a variant of the villin headpiece C-terminal fragment, a small all-α protein (24), to within an RMSD of ~1 Å from the crystal structure. (These results cannot be taken as evidence that the same force field would necessarily succeed in folding other proteins.)

Reversible folding and unfolding in an equilibrium simulation of FiP35. In the hope of observing a number of folding and unfolding events under equilibrium conditions, we then ran two independent 100-μs MD simulations of FiP35 at a temperature (395 K) approximating the protein’s in silico melting temperature. [The in silico melting temperature was estimated using a replica exchange metadynamics simulation (25) and is ~40 K higher than the experimental melting temperature.] In each of the two equilibrium simulations, FiP35 underwent multiple folding and unfolding transitions, for a total of 15 events in the two simulations (Fig. 2A). We found the folding time, as calculated from the average waiting time in the unfolded state, to be 10 ± 3 μs, in relatively close agreement with the experimental folding time [14 μs (12)]. The population of the folded state in the simulations was 60%.

Fig. 2

Reversible folding simulation of FiP35. (A) RMSD time series of two independent 100-μs simulations of FiP35 initiated from an extended state. RMSD with respect to the x-ray structure (15) was calculated for the Cα atoms of residues 4 to 32. A total of eight folding and seven unfolding events can be observed within the two simulations. (B) Representative sequence of events leading to folding (left) and unfolding (right). RMSD to the x-ray structure was calculated for different regions of the protein, namely the tip of hairpin 1 (residues 12 to 18, blue), the entire hairpin 1 (residues 8 to 22, green), hairpin 2 (residues 19 to 30, orange), and the full protein (residues 2 to 33, red). Analyses of the individual transitions reveal a consistent sequence of events leading to folding (upper panel): The tip of hairpin 1 forms first, followed by hairpin 1, followed by hairpin 2, followed by the rest of the protein. The reverse sequence of events is observed in unfolding transitions. (C) Representative members of the transition-state ensemble reveal that the first, but not second, hairpin is structured in the TSE for folding. The transition-state ensemble was identified from equilibrium simulations using a previously described procedure (29), and representative structures were superimposed using Theseus software. (D) Commitment probability (Pfold) distribution of the TSE as calculated from 101 structures, using four simulations for each structure. The observed distribution of Pfold (red) is compared with the binomial distribution expected for a true TSE (black). (E) Comparison of experimental and calculated ϕ values. Two sets of ϕ values (red and green) were calculated independently from the two simulations and are compared to the experimental values (black) obtained for wild-type Pin1 WW domain (17). Error bars in (D) and (E) correspond to ±1 SEM.

In our simulations, a well-defined sequence of events leads from the disordered, unfolded state to the native state: In all folding transitions, formation of the tip of the first hairpin is followed by formation of the entire first hairpin, then by formation of the second hairpin, and finally by consolidation of the hydrophobic core (Fig. 2B, top). Unfolding transitions follow the reverse pattern. Thus, the folding mechanism of FiP35 appears to be dominated by a single pathway, with any flux through alternative pathways being small. This is in contrast with recent studies of WW domains (2022)—using simulations shorter than the folding time—which found structurally heterogeneous dynamics with multiple parallel pathways connecting folded and unfolded states.

To explain why FiP35 folds along a single, dominant route, we performed reversible folding simulations of peptides corresponding to the two hairpins of FiP35. Both the first and second hairpin, in isolation, fold to the structure found in the full WW domain, with similar time constants (1.2 ± 0.1 μs and 0.86 ± 0.06 μs for hairpins 1 and 2, respectively). The two hairpins, however, differ substantially in their stabilities, with the population of the folded state being 25% and 4% for the first and second hairpins, respectively. Thus, the order in which the individual hairpins form during folding of the full WW domain mirrors their intrinsic thermodynamic stabilities, with the slower but more stable hairpin 1 forming first. We suspect that the alternative pathway in which the second hairpin forms first may not be substantially utilized because of the relatively large difference in stability between the two hairpins (~1.5 kcal mol−1). Indeed, experiments have shown that variations in substructure stabilities can cause substantial shifts in folding pathways (26, 27), as have computational studies in the context of the diffusion-collision model for protein folding (28).

Our simulations show that both hairpins also form with similar rate constants in the context of the full WW domain. The first hairpin forms with a time constant of 5 ± 1 μs, slower than in isolation by a factor of 4; such a slowdown is greater than expected for a simple hydrodynamic drag. To determine its origin, we performed simulations in which we systematically added residues in silico, three at a time, to the two ends of the hairpin (25). Adding 12 C-terminal residues resulted in a slowdown by a factor of 2. There was no change in the folding rate upon adding the first three N-terminal residues (residues 1 to 3), but the addition of residues 4 to 6 slowed hairpin formation by a factor of 2. Thus, approximately half of the observed slowdown is determined by the addition of just a few residues.

Characterizing the transition state for folding. MD simulations can in principle provide detailed insight into all steps of the protein-folding pathway. Experiments, however, are generally limited to the characterization of stable states and rate-limiting transition states. To facilitate comparison with experiment, we determined the transition state in our simulations. We used a previously described variational approach (29) to find an optimized reaction coordinate that separates the transition state from the stable folded and unfolded basins (25). To validate the resulting transition-state ensemble (TSE, Fig. 2C), we calculated the commitment probability (Pfold), the probability that simulations initiated from a given structure will fold before unfolding. For each of 101 structures sampled at random from the TSE, we ran four simulations until either the folded or unfolded state was reached, requiring in total an additional 151 μs of simulation.

The distribution of Pfold obtained with this approach peaked at 0.5, closely resembling the binomial distribution expected for an ideal TSE (Fig. 2D), thus validating the reaction coordinate and the TSE. The average commitment time observed in the simulations initiated from the transition-state region is relatively long (0.36 μs and 0.40 μs for folding and unfolding, respectively), suggestive of a diffusional process over a flat free-energy barrier. Inspection of the TSE shows that formation of the first hairpin, but not the second, is part of the rate-limiting step. We note also that the time constant of formation of the first hairpin determined above (5 ± 1 μs) is half of that for folding (10 ± 3 μs), consistent with our finding that formation of this hairpin is rate-limiting in folding. (The factor of 2 difference in time constants arises from the 50% probability of folding from the TSE.)

The most important experimental strategy for inferring the structural properties of the TSE is the protein engineering method, in which folding and unfolding rates are measured for a series of mutants. The results are typically presented as ϕ values (30). A ϕ value of ~1 suggests that the interactions formed by a residue in the native state are also present in the TSE, whereas a value close to zero indicates that the native-state interactions are not present in the TSE. Commonly, ϕ values are calculated from simulation by approximating the mutational free-energy changes from the fraction of native side-chain contacts lost upon mutation (31). We used this approach to calculate ϕ values for side chains and a free-energy perturbation approach (25) for the backbone, and compared the results to experimental measurements in a related WW domain (Fig. 2E) (17). The values obtained confirm the observation that the first hairpin is essentially fully formed in the TSE, whereas the second hairpin only makes a fraction of the contacts found in the native state.

Although the agreement between simulation and experiment (Fig. 2E) is encouraging, it may also be somewhat fortuitous, as the number of atomic contacts is only a rough approximation for the free energy (32). We thus also calculated ϕ values directly from the folding and unfolding rates obtained from reversible folding simulations of mutant proteins (33, 34) and compared these results to ϕ values calculated by applying the contact approximation to the TSE for FiP35.

We chose to study the effects of six mutations that we expected to have different effects on the folding and unfolding rates (Fig. 2E). Ser13 → Ala, located in the tip of the first hairpin, and Arg11 → Ala, located in the central β strand, are expected to have high ϕ values; Tyr19 → Leu and Phe21 → Leu, both located in the central β strand, are expected to have intermediate ϕ values; and Leu4 → Ala and Trp8 → Phe, located in the hydrophobic core, are expected to have low ϕ values.

All six mutants folded reversibly to the native state, albeit with different rates and stabilities (Table 1). For most mutants, the changes in stability upon mutation were in good agreement with experimental data. Most of the ϕ values calculated from the folding and unfolding rates were in reasonable agreement with the magnitude calculated from the TSE of FiP35 (Table 1), although our results support the notion that individual ϕ values are best interpreted qualitatively (35, 36).

Table 1

Computational ϕ-value analysis of FiP35. In columns 2 and 3, the ϕ values for six selected mutants, calculated from the folding and unfolding rates obtained from reversible folding simulations, are compared with the values estimated from a contact approximation. In columns 4 and 5, the calculated free-energy changes upon mutation are compared to the values measured experimentally at the melting temperature of the hPin1 WW domain (49).

View this table:

A notable exception is the Arg11 → Ala mutation, whose low ϕ value calculated from the folding kinetics (0.2) is substantially smaller than the high value expected from the contact approximation (0.8). Although the reasons for this discrepancy remain unclear, our results overall support the use of experimentally derived ϕ values to infer the structural properties of the TSE, with the caveat that the contact approximation may fail in individual cases (35); all the experimental ϕ values should thus be considered simultaneously when inferring the overall structural properties for a TSE (31). Finally, we note that all mutant proteins fold via the same overall pathway as FiP35, although some of the mutations appear to cause a noticeable Hammond-like shift in the structure of the TSE (fig. S4).

It is perhaps worth noting that these simulations may also provide a computational gold standard for future studies exploring the accuracy and efficiency of methods for the prediction of mutational free-energy differences and folding rates.

FiP35 folds across a small free-energy barrier. We determined the free-energy profile and position-dependent diffusion constant along the optimized reaction coordinate (25). We found that the free-energy barrier for folding is small (1.6 kcal mol−1 or ~2kBT), consistent with the suggestion that FiP35 is an incipient downhill folder (12). The transition-state region is broad and flat (Fig. 3A), helping to explain the long commitment time observed in the Pfold analysis. Langevin simulations on the free-energy profile (Fig. 3B) approximate well the folding dynamics observed in the MD simulations, and we are thus able to use this kinetic model to simulate a temperature-jump experiment (25). In addition to the slow phase associated with folding, we observed a fast “molecular” phase whose amplitude and time constant depend on both the size of the temperature jump and the spectroscopic probe used (25). Such features are spectroscopic indications of protein folding across a low free-energy barrier, and they support the notion that experimental studies of fast-folding proteins might be used to probe directly the spectroscopic properties of the TSE (37).

Fig. 3

Folding kinetics across a low energy barrier. (A) Free-energy profile along an optimized reaction coordinate. The profile exhibits two minima, centered at 0.1 and 0.7, corresponding to the folded and unfolded basins, respectively. The folding and unfolding free-energy barriers are 2kBT and 3.5kBT, respectively. (B) Langevin simulation of WW folding in a one-dimensional model. The simulation was based on the one-dimensional free-energy profile in (A) and a position-dependent diffusion coefficient, both derived from the MD simulation data.

It has been argued that the fast molecular phase provides an estimate of the time scale for transition paths in folding of FiP35 (37). The value obtained in these experiments (≤0.7 μs) is in agreement with theoretical estimates (0.3 μs) as well as with the upper bound (200 μs) obtained directly through single-molecule experiments (6). These values also agree with the average transition path time observed in our equilibrium simulations (0.4 ± 0.1 μs). Thus, a range of different techniques (simulation, theory, ensemble, and single-molecule experiments) provide independent evidence for transition path times for protein folding on the order of 1 μs.

Native-state dynamics of BPTI. Dynamic changes in protein structure typically occur not only during but also after the folding process. The 58-residue protein BPTI was the subject of the first nuclear magnetic resonance (NMR) experiments of the internal motions of proteins (38). NMR studies showed that on time scales ranging from nanoseconds to milliseconds, several internal water molecules exchange with the bulk (39, 40), a number of aromatic rings rotate (38, 41), and a disulfide bridge isomerizes (42, 43). We used a 1-ms MD simulation at a temperature of 300 K to reproduce and interpret the kinetics of folded BPTI.

Our simulation of BPTI transitioned reversibly among a small number of structurally distinct long-lived states (Fig. 4A) of lifetimes ranging from 6 to 26 μs. The two most populated states in the simulation, which together accounted for 82% of the trajectory, are supported by experimental data. The average structure of the “crystallographic” state, which was occupied 27% of the time in our simulation, had an RMSD (between nonsymmetric heavy atoms of residues 4 to 54) of only 0.8 Å from the crystal structure (44). The most populated state (occupied 56% of the time) exhibits a left-handed conformation of the disulfide bridge between Cys14 and Cys38, which has been previously suggested by NMR experiments (42, 43). The imbalance in populations between simulation and experiment represents an error in the conformational free-energy difference of just a few times the thermal energy, which falls within the expected accuracy range of the force field representation.

Fig. 4

Native-state dynamics of BPTI. (A) All-residue backbone RMSD from the crystal structure with PDB ID 5PTI (44). Each point shows the median value in a window of 50 ns. The color of the data points denotes cluster membership. [See (25) for PDB files representative of each of the five most-visited states in the 1-ms simulation of BPTI, together with a more detailed analysis of the local structural features that best discriminate each state from the others.] (B) Dynamical content of the P2 internal correlation functions (25) and its decomposition into side-chain and backbone contributions. The peak near 28 ps results from the relaxation of the side chains within a conformation, whereas the peak around 10 μs results from the relaxation of the backbone during jumps between states. The cartoon shows the fractional contribution of each residue to the decay of the average internal correlation function between lag times of 10 ns and 20 μs. (C) Crystal structure of BPTI, highlighting the aromatics that rotate slowly in purple and those that rotate quickly in orange, with representative structures from each of the four additional conformations observed in the simulation. (D) Survival probability distributions for each of the four internal water molecules of BPTI. The arrow at 14 μs marks the lifetime of the slowest waters, as determined from a double-exponential fit of the tail of the W122 survival times.

In addition to the two states anticipated by previous experiments, the systematic analysis of the dynamics provides evidence for the existence of at least three additional states. In some, the exposed surface area of BPTI was large, as the conformation displayed an additional water-filled cavity or even a pore; in others, the small N-terminal helix unfolded.

Separation of time scales in protein dynamics. As part of our analysis of the native-state dynamics of BPTI, we performed dynamical content calculations, derived from the autocorrelation of bond orientations over a very wide (and previously inaccessible) range of time lags, to quantify the structural relaxation occurring on different time scales and the persistence of various structural features (25). These analyses revealed a distinct separation of time scales: Hopping among conformational basins occurs on time scales on the order of 10 μs, whereas motions within the individual basins occur on a time scale that is faster by several orders of magnitude. Substantially less dynamical content is measured within the wide gap that lies between these two time scales (Fig. 4B).

This separation of time scales, a hallmark of barrier-crossing dynamics (5), allowed us to partition the trajectory into distinct long-lived states (color-coded in Fig. 4A) using a new kinetic clustering scheme, which was designed to retain important aspects of the long–time scale behavior observed in the MD simulation (25).

We found that the fast relaxations, which extend up to the 10-ns time scale, originate primarily from side-chain motions, whereas the slow relaxations—corresponding to hops between well-separated basins—originate almost exclusively from backbone motions (Fig. 4B). Apart from low-amplitude vibrations with frequencies several orders of magnitude higher than those of conformational state transitions, bond orientations within the backbone exhibit little variation within a given state; motion is “frozen out” not only at the level of the protein’s global structure, but at a local level within the polypeptide chain. Side chains, on the other hand, experience large fluctuations on the nanosecond time scale; these fluctuations, however, are nearly identical in all conformational states, leading to the absence of a long–time scale signature in the side-chain dynamical content. These findings provide a unifying model, encompassing an extremely wide range of time scales, that integrates the previously established picture of “liquid-like” side chains attached to a “solid backbone” (45, 46) with the common observation of distinct conformational states (2, 3).

The transition path time for conformational transitions was generally at least several hundred nanoseconds (25). This time—which might, if anything, tend to increase with protein size—serves as a lower bound for the lifetime of individual conformational states. Even for a protein as small as BPTI, this lower bound is more than an order of magnitude larger than the time scale of intrabasin side-chain motions. This observation suggests that a distinct separation of time scales, with a region of reduced dynamical activity in the interval between those time regimes, is likely to be a common feature of the dynamics of folded proteins.

Local probes that report on large-scale conformational change. Seven of the eight aromatic rings of BPTI rotated in the simulation (25) with rates that generally agree with experiment (38, 41). We found that the rings that rotated slowly in the simulation were all located in the portion of the protein that did not change during hops between basins (purple side chains in Fig. 4C). In contrast, the rings that rotated quickly in the simulation were located in the portion of the protein that changed in the different conformational states (orange side chains in Fig. 4C) (25).

In agreement with experiment (40), the simulation displayed three distinct behaviors for the four internal water molecules (Fig. 4D), with W111 exchanging very rapidly with the solvent, W112 and W113 exchanging on a slower time scale, and W122 exchanging even more slowly. The tail of the survival probability of W122 displays single-exponential decay with a lifetime of 14 μs, suggesting that the longest-lived waters have a single escape mechanism. The binding events of the eight longest-lived instances of W122 (each bound for more than 9 μs) all occurred while the simulation was in the “crystallographic” state. Accordingly, the lifetime of this state in the simulation (25 μs) is longer than the lifetime of W122 (14 μs). The lifetime of the longest-lived waters may thus provide an experimentally measurable lower bound on the lifetime of the protein conformation that most tightly binds to that water.

The water escape path lasted only a few picoseconds. Before the escape, a small pocket opened on the protein surface that allowed a single water molecule to form a hydrogen bond with W122; W122 took the place of that hydrogen-bonded water, and a transient empty cavity was created inside the protein. The wide contrast between the lifetime of W122 and the duration of its exit and entry pathways suggests that ligand entry and escape can proceed without the need for long-lived intermediate conformations; thus, binding and escape events can be considerably faster than the lifetime of the ligand.

Transition path times for conformational changes. Conformational changes such as those described here for FiP35 and BPTI differ from simpler chemical reactions in that the former occur on a rough free-energy landscape characterized by many local minima. The roughness of the landscape manifests itself in a substantial slowdown of conformational changes and in a kinetic preexponential factor, k0, that is many orders of magnitude larger than the value expected for elementary reactions in small molecules. Using the folding rate and barrier height we obtained for FiP35, we estimate the preexponential factor to be ~1 μs−1. A related quantity is the average transition-path time 〈tTP〉 (i.e., the time required for the actual conformational change to occur). This value is also expected to be sensitive to the landscape roughness. We obtained a 〈tTP〉 value of 0.4 μs for the folding reaction of FiP35 and 0.3 μs for the conformational changes observed in BPTI.

Indirect estimates for k0 and 〈tTP〉 have previously been proposed on the basis of theory and experiment (6, 12, 24, 37, 47, 48), but these fundamental quantities have proven difficult to measure. The microsecond time scales associated with these constants are not only in general agreement with previous estimates, but also close to the observed time scales for conformational changes in FiP35 and BPTI, thus strengthening the notion that landscape roughness is a major determinant of the rates for conformational transitions in biological macromolecules.

Conclusions. The specialized machine we developed has allowed us to perform continuous, all-atom molecular dynamics simulations of proteins in an explicitly represented solvent environment over periods as much as 100 times longer than was previously feasible. Comparison of the results of these simulations with experimental measurements provides evidence for the non-obvious finding that existing force fields are capable of realistically describing of the structure and dynamics of proteins over even these extended time scales. More generally, our findings suggest that very long molecular dynamics simulations can serve as a powerful tool for elucidating the atomic-level behavior of proteins on a biologically critical but previously inaccessible time scale.

Supporting Online Material

Materials and Methods

Figs. S1 to S14

Tables S1 to S3

PDB files S1 to S5

References and Notes

  1. See supporting material on Science Online.
  2. We are very grateful to all members of the Anton hardware and software teams, without whom this work would not have been possible. We thank G. Hummer for providing us with the software to calculate position-dependent diffusion constants, A. Pan for helpful suggestions, T. Tu for assisting with trajectory analysis, A. Philippsen for helping with the BPTI renderings, K. Mackenzie for monitoring and supporting the BPTI simulation, and R. Kastleman and J. McGrady for editorial assistance.

Stay Connected to Science

Navigate This Article