Folding and misfolding of the papillomavirus E6 interacting peptide E6ap

Implicit Solvent LD Simulation Methods

The LD simulations take the total system energy U _total = U _b + U _bend + U _tors + U _{imp - tors} + U _ch(ε(r)) + U _vdW + U _solv(σ) as the sum of the usual types of interaction potentials between the protein atoms, whereas the solvent contributions are modeled by using a spatially dependent dielectric constant to screen charge–charge interactions U _ch(ε(r)) and a solvation potential U _solv(σ). The bonding interactions U _b, bond–bond bending interactions U _bend, and improper torsional energies U _{imp - tors} are modeled by harmonic potentials, the regular torsional potentials U _tors by standard periodic functions, and the van der Waals interactions by Lennard–Jones 6–12 potentials. The Coulomb interactions U _ch(ε) = Σ_{i > j} q _i q _j/ε(r _ij)r _ij are expressed in terms of atomic partial charges q _i and a Ramstein–Lavery style (38) distance-dependent dielectric constant ε(r). The microscopic solvation potential is modeled by using the Ooi–Scheraga solvent accessible surface area potential (39), , where σ_i is the accessible surface area of a hypersurface bisecting the first solvent shell surrounding protein atom i, and the empirical surface free energy parameters g _i depend on the atom type. Because the g _i are free energy parameters, the U _total generates a temperature-dependent free energy that contains contributions from solvent reorientation within a mean-field approximation.

The LD simulations use the velocity Verlet algorithm (40) with a time step of Δt = 2.5 fs for integrating the equations of motion for the protein atom positions and velocities. The computations are performed by using a modified version of the tinker 3.9 molecular design package (41) and the standard amber94 force field (42). The frictional forces and corresponding random forces acting on the protein atoms are computed by using the Pastor–Karplus accessible surface area model (43). The solvent accessible surface areas σ_i for the friction coefficients are calculated from the exposed surface area of the solute atom by using a probe of zero radius. The smaller probe size for friction coefficients is used to cancel effectively the results of (more expensive to calculate) hydrodynamic interactions. The accessible surface areas, atomic friction coefficients, and solvation potentials are updated every 100 dynamical steps (0.25 ps), because tests show that this approximation incurs negligible error because significant conformational variations occur on a much longer time scale (29). No cutoff is applied on the forces between nonbonded atoms.

All calculations are performed on a single central processing unit (CPU) (Pentium 4, 1.6 GHz/256 Mb), with an output of ≈30 ns per CPU day for E6ap whose sequence is IPESSELTLQELLGEERR. Before starting each simulation, the initial conformation is energy minimized with a steepest descent algorithm to remove close contacts that would otherwise be present on cooling the initial configuration to 300 K (see Simulations of E6ap Dynamics). A constant temperature is maintained by using a Berendsen-type (44) thermal bath coupling throughout the simulation, and the data are saved every 4,000 steps (10 ps) for subsequent analysis.

Simulations of E6ap Dynamics

Before beginning the folding simulations of E6ap, we first test its stability by conducting LD simulations at 300 K that start with a representative NMR structure in the presence of 40% TFE. An ensemble of 24 NMR solution structures is obtained from the Protein Data Bank (PDB ID code 1EQX). Of this group of structures that all satisfy the NMR nuclear Overhauser effect constraints within a defined tolerance, a single representative native structure is selected by using the nmrclust module (45) from the olderado (On-Line Database Ensemble Representatives And DOmains). Among the 24 NMR structures, structure 16 is chosen to be the representative structure for comparisons. We have run two independent 50- and 100-ns trajectories beginning from the same NMR structure but with different random seeds for the initial velocity distribution and random forces.

The similarity of the protein structure in a simulation to its native conformation is characterized by the rms deviation (rmsd), which measures the minimized rms distance of each backbone C_α atom from its native structure. Because the structure of E6ap in water is currently unavailable from experiment, the NMR structure in the mixed solvent is chosen as the “native” conformation solely to have a fixed reference for computing the rmsd. The backbone rmsd between the instantaneous and NMR native protein conformations, for both the whole peptide and the central segments (residues 4–13), increases initially, reaches a constant after ≈2 ns, and remains close to that level throughout the simulation. The average of the saturated limit of the rmsd for the central segments (1.5 Å) is much smaller and exhibits less fluctuations than that for the whole peptide (≈3.8 Å), which suggests that the average structure of the central segments is more stable, whereas the two ends are very mobile. Nevertheless, these values for the rmsd indicate the expected departure of the conformation in aqueous solution from the “native” NMR structure in 40% TFE solution.

Representative conformations adopted by E6ap are shown in Fig. 1 b–e. Compared to the initial NMR structure in 40% TFE solution (in Fig. 1 a), the E6ap peptide retains a fluctuating α-helix structure in the central portion (Ser-4–Leu-13), where the essential leucine patch is well conserved throughout the two trajectories. However, both the N and C termini are dynamically very flexible. The equilibrium fluctuations of the helical content occur quite rapidly, normally within 1 ns, which suggests a highly labile helical dynamics in E6ap. As shown in Fig. 1 b–e, each terminus can extend the helix structure either by forming a backbone hydrogen bond or by flipping back to interact with the central portion through side-chain interactions and therefore failing to produce a discernible secondary structure. Also, it appears that the back-and-forth flippings of the N and C termini are quite independent.

To study the E6ap folding mechanisms and the possible folding pathways and intermediates, LD folding simulations are conducted with E6ap starting from completely denatured states that are obtained by first running LD simulations at a high temperature (750 K) for about 1 ns each. After the 1-ns denaturation, the structures are quenched and free energy minimized at 300 K to serve as the initial conformations for the LD folding simulations. (Note that this method of denaturation is more extreme than any accessible in experiments.) Six independent folding trajectories, with time durations of 104, 105, 143, 179, 281, and 1,004 ns, respectively, have been generated starting from different initial states. The longest 1-μs trajectory is used to construct the relative free energy plots in Free Energy Landscape and Folding Path.

The denatured conformations are fairly extended, and all secondary-structure and native contacts are removed after heating at 750 K for 1 ns. The essential leucine hydrophobic core (Leu-9, Leu-12, and Leu-13) may either be preserved or destroyed after heating. At the end of the simulations, all six folding trajectories exhibit a loosely fluctuating central helix and two flexible tails, similar to the trajectories that begin from the NMR structure in 40% TFE. However, the folding paths for the simulations that begin with configurations that lack a leucine core are qualitatively different from the other simulations, as discussed below.

To elucidate the dynamics of the fluctuating central helices in pure aqueous solution, an analysis of the average helical content is performed by using a time average computed from the 1-μs trajectory. A residue is defined as being in the α-helix state at time t if the (ϕ, Ψ) angles fall within an angular range ±30° about (-57°, -48°), and the residue forms an (i, i + 4)-type backbone hydrogen bond with other residues. The second criterion is used to exclude cases where the instantaneous (ϕ, Ψ) angles form a right-handed α-helical structure without having structural support from hydrogen bonds. To compute the average helical content from the LD trajectory, it is convenient to define a variable ξ_μ(t _i) that is unity if residue μ is in a helical state at the ith time step t _i and that vanishes otherwise. The average of ξ_μ(t _i) over t _i and μ gives where n is the number of residues, T is the number of time frames in the trajectory, and h̄ _μ is the average helical content of the μth residue. Fig. 2 presents h̄ _μ as a function of residue μ by using the last 500 ns of the 1-μs trajectory (to assure that E6ap is well equilibrated). The flexible nature of the helical structure of E6ap in pure aqueous solution is quite evident from Fig. 2, which exhibits residues Leu-9 and Leu-12 as a center of helical stability, Glu-11 as a center of minimal helical stability, residues Pro-2 through Glu-6 as weakly helical, and the C-terminal residues Gly-14 through Arg-17 as virtually devoid of helical content. Using this stringent definition for the α-helix, the average helical content h̄ emerges as 22.3% and is lower than the helical contents deduced from other popular definitions [31.5% with the Kabsch–Sander definition (46), 39.5% with a criterion depending only on Ψ/ϕ (26)]. The amber94 force field is known to overemphasize helical propensities, and the simulation temperature of 300 K is higher than the experimental condition of 288 K, where the helical propensity would be lower. Therefore, combining the experimental CD helical measurement (20%) and the simulated value (22–40%), we conclude that the E6ap adopts a fluctuating helical structure in pure water. Moreover, computations of the correlation time τ_μ for survival of helical content in the μth residue support the fluctuational description of a loose central helix. Fig. 2 displays τ_μ in the nanosecond range but with no correlation between τ_μ and h̄ _μ. For example, Thr-8 has the third largest h̄ _μ but one of the smallest τ_μ, whereas Glu-3, Glu-6, and Leu-7 exhibit the largest τ_μ but only moderate h̄ _μ, suggesting that the residues participate in a stable helix that forms rarely but persists relatively long. The high helical stability of Leu-9 and Leu-12 is vital for the p53-binding activity of E6ap, because this stability maintains the structural integrity of the leucine-triplet patch that is required for binding to E6. Despite the low average helicity of E6ap in pure aqueous solution, this LD simulation displays the structure of the binding site as well preserved with relatively high local helix propensity.

Further analysis of the dynamics proceeds by viewing the time evolution of the radius of gyration (R _g) and the backbone rmsd for both the whole peptide and for the central segments (residues 4–13). In Fig. 3a (for one of the five folding trajectories), the rmsd curves initially drop sharply during the first few nanoseconds and then fluctuate around the same constant values (3.7 Å for the whole peptide and 1.5 Å for the central segments), as in the simulations beginning from the NMR structure in 40% TFE solution. The time course of the R _g exhibits behavior similar to the rmsd and saturates around an average of 9 Å. However, the rmsd curves in Fig. 3b (misfolding trajectory) initially increase to much larger constant values (7 Å for the whole peptide and 4 Å for the central segment) and remain there for >80 ns. Then, at ≈85 ns, they rapidly drop (in 3 ns) to the equilibrium values as in the stability simulations. The curves for R _g in Fig. 3b initially decrease to a smaller value (≈8 Å) and then suddenly increase back to the equilibrium value (≈9 Å) at around 85 ns. Note that except for the first few nanoseconds, the fluctuations in R _g are anticorrelated with those in the rmsd, which means the native state is not the most compact state, in contrast with the situation generally found for large proteins.

The folding process can be described more quantitatively by considering the time course of the fraction of the native contacts Q. Two residues are considered to be in contact when the distance between their α-carbons is within 6.5 Å, and they are separated by at least two residues in the sequence (to remove interferences from bonding constraints) (47). This conventional definition of the fraction of native contacts suffers from some shortcomings. By definition, Q can equal only discrete rational values n/[total number of native contacts], where n is an integer from zero to the total number of native contacts. The time course of Q along a single trajectory is very noisy when the total number of native contacts is small (only 38 for E6ap), and the information conveyed is partially obscured by the large fluctuations. Hence, we introduce a new parameter P as a measure of the average probability for the presence of native contacts. Instead of a sharp cutoff 6.5 Å for the contact distance as in the calculation of Q, the native contact probability for residues i and j is defined by a Gaussian distribution with mean and with a spread proportional to the interparticle distance in the native structure. Thus, the average total probability of finding native contacts at time t is defined as where r _ij(t) is the α-carbon distance (in angstroms) at time t between residue i and residue j (|i - j| ≥ 2), and N is the total number of residue pairs treated. α is an empirical parameter that controls the penalty for deviations from the mean separation distance in the native state, and a value of 14.0 for α is obtained by minimizing the difference between P and Q. As shown below, P is continuous, less scattered, and provides much clearer information than Q for the structural details arising from long-distance (>6.5 Å) contacts.

Fig. 3 c and d displays the fraction of native contacts Q and the average probability of native contacts P as a function of time for a typical folding trajectory and for the misfolding trajectory. The behavior of the parameter P agrees very well with Q but contains significantly reduced fluctuations, so attention is focused on P for the remainder of the discussion. During the folding trajectory, P quickly grows from 0.2 to 0.7 within the first 10 ns and then fluctuates between 0.6 and 0.8 ns for the rest of the trajectory, which accords well with the initial decrease and the subsequent dynamics of the rmsd in Fig. 3a. P for the stability trajectories (not shown) exhibits similar fluctuations between 0.6 and 0.8 ns. On the other hand, P for the misfolding trajectory experiences an initial rapid decrease to 0.3 ns, then slowly but steadily grows back to ≈0.7 ns, and at ≈90 ns jumps in a similar although less dramatic behavior than the rmsd but with several nanoseconds lag. Because the time evolution of P more or less follows the side chain dynamics and because the rmsd tracks only the backbone movement, the difference in their time courses suggests that the side chain movement is slow and continuous, whereas alterations in the backbone structure are more discrete and occur on a shorter time scale.

The behavior of both the rmsd and P implies that dramatic changes occur in the E6ap conformation at ≈85 ns in the misfolding trajectory. Typical spatial configurations for E6ap before 85 ns (the misfolded structure) and the NMR structure are depicted in Fig. 4, with the core residues Leu-9, Leu-12, and Leu-13 displayed in atomic detail. Although the NMR structure in 40% TFE solution has a single central helix with the leucine core on the same side of surface, the misfolded structure is composed of two smaller helices located at separate arms of a hairpin-like structure. Also evident in Fig. 4 is the absence of the leucine hydrophobic core in the misfolded helical hairpin, with the three core leucine residues situated at different sides of the backbone chain and especially with Leu-9 not involved in either of the two helices. The role of the leucine triplet is better illustrated in terms of solvent accessible areas. The misfolded state has a particularly large surface area of 322.5 Ų for Leu9-12-13, nearly as large as the 323.8 Ų for the unfolded state, and much greater than the 298.1 Ų for the native structure. This surface burial of the leucine triplet and the favorable Lennard–Jones interactions on clustering not only stabilize the native state but also guide the misfolded structure toward it.

Free Energy Landscape and Folding Path

The misfolding trajectories display the peptide as initially trapped in a misfolded state before folding back into the somewhat loose native state. To elucidate the nature of the intermediate misfolded state more fully, we conduct a long trajectory (≈1 μs), beginning with the same initial conformation as in the previous misfolding simulation but with a different random number seed for the initial velocity distribution and random forces. The time course of the 1-μs trajectory exhibits the same initial trapping (rmsd_4–13 ≈ 3.5 Å) in the misfolded state, although the transition time varies slightly, with a later return to the equilibrium state (rmsd_4–13 ≈ 1.5 Å). Notably, the two misfolding trajectories both display a conversion time from the misfolded to folded states at ≈80–100 ns. Subsequent to the initial trapping before 90 nsec, the rmsd for the long trajectory increases again toward the misfolded state value before finally settling down to that characteristic of the native state central portion value at t ≈ 350 ns. The protein conformation at t ≈ 300 ns appears to be similar to that at ≈70 ns, as displayed in Fig. 4 b, which indicates that the peptide revisits the same misfolded state reversibly.

A full understanding of the folding mechanism requires a global overview of the free energy landscape. However, it is difficult to construct the free energy surface for several reasons. First, the free energy depends on a large number of variables. To deduce meaningful information, the multidimensional surface is generally projected onto some suitable folding coordinates (9, 15). Second, to compute the free energy, a complete sampling of conformational space is necessary. Although umbrella-sampling (48, 49) and multicanonical or entropy-sampling techniques (5, 7) appear very promising, a large number of simulations are required just to map out the structure of the free energy surface near the native structure basin.

We have computed the relative free energy for the folding free energy landscape of the E6ap peptide. The relative free energy is taken to be proportional to the logarithm of the probability for finding the peptide with specific values of the chosen folding coordinates (21–23) during the long 1-μs trajectory. The folding coordinates are the rmsd_4–13 and R _g, whose instantaneous values during the folding process are projected onto a 2D contour map. Even though all the folding simulations yield the same free energy basin as in the simulations that begin from the NMR structure, only the 1-μs folding trajectory is long enough that the system crosses the barrier between folded and misfolded free energy basins more than once. Hence, the following relative free energy plots are understood as being only qualitative.

Fig. 5 displays the contour plot of the free energy surface as a function of the rmsd and R _g for the 1-μs trajectory. The large and flat native basin at rmsd_4–13 ≈ 1.5 Å and R _g ≈ 9 Å for the equilibrium simulation (not shown) indicates a flexible and dynamical native state. The folding trajectory produces a single energetic descent toward the native basin. Two free energy minima are clearly evident in the free energy surface for the 1-μs misfolding simulation (Fig. 5). In addition to the same minimum at rmsd ≈ 1.5 Å and R _g ≈ 9 Å, another minimum exists at rmsd_4–13 ≈ 3.5 Å and R _g ≈ 7.5 Å, thereby demonstrating the existence of a metastable state with relatively lower R _g and larger rmsd.

Further inspection of the energy components at the misfolded and folded free energy basins provides better insight into the underlying forces during the folding-misfolding transition. The total potential energy for the long 1-μs simulation displays an overall descent throughout the simulation, varying sharply at the critical time of the misfolding–folding transition. The Coulomb and solvation energies, however, exhibit quite different behaviors. The solvation free energy is positively correlated whereas the Coulomb energy negatively correlated with respect to the behavior of rmsd. The instantaneous peptide conformation in the 1-μs trajectory is classified as being in the misfolded or native free energy basins if its R _g and rmsd values fall into the ranges described Table 1. Table 1 also presents average values of the total energy, Coulomb energy, and solvation energy for E6ap in the misfolded and native free energy basins. The average total energy of the folded state is ≈7.3 kcal/mol lower than that of the misfolded state. The average solvation free energy is 2.4 kcal/mol lower, whereas the average Coulomb energy is 6.0 kcal/mol higher for misfolded than for the native state. The lower solvation energy is as expected because of the more compact structure of the misfolded state (see Fig. 4 b). The excess electrostatic energy in the folded state is mainly due to the formation of more backbone hydrogen bonds in the core α-helix, and this energy contribution finally drives the peptide back to the folded structure.

Discussion

We have examined the dynamics of the folding of the E6ap peptide by using implicit solvent LD simulations. Our results indicate that: (i) E6ap adopts a flexible helical structure in the LD simulation, and the leucine hydrophobic core lies on the same face of the helix, both of which agree well with NMR experiments; (ii) the helical content derived from the long time simulation is 22–40%, which overestimates, as expected, the 20% from CD measurements and suggests the reliability of this description; (iii) the central leucine core is not only crucial for its biological function but also important for the stability of native structure; (iv) starting from different initial configurations, there are multiple folding pathways leading to the native state, and E6ap can behave either as a two- or three-state folder; (ν) the free energy surface contains a significant barrier between the native and the misfolded state basins; and (vi) the Coulomb energy gained in forming backbone hydrogen bonds is the main driving force for the peptide to fold back into its native structure. The richness observed in the diversity of folding pathways for this simple 18-residue peptide suggests that perhaps even greater richness is to be found from simulations for the folding of larger and more complex proteins through computations that are enormously facilitated by the use of accurate, validated, implicit solvent models and by continually increasing computational power.