Structural Ensemble of CD4 Cytoplasmic Tail (402–419) Reveals a

Article pubs.acs.org/JPCB

Structural Ensemble of CD4 Cytoplasmic Tail (402−419) Reveals a Nearly Flat Free-Energy Landscape with Local α‑Helical Order in Aqueous Solution Navjeet Ahalawat,† Simran Arora,†,‡ and Rajesh K. Murarka*,† †

Department of Chemistry, Indian Institute of Science Education and Research Bhopal, Bhopal By-pass Road, Bhauri, Bhopal 462066, Madhya Pradesh, India S Supporting Information *

ABSTRACT: The human cluster determinant 4 (CD4), expressed primarily on the surface of T helper cells, serves as a coreceptor in T-cell receptor recognition of MHC II antigen complexes. Besides its cellular functions, CD4 serves as a primary receptor of human immunodeficiency virus (HIV) type 1. The cytoplasmic tail of CD4 (residues 402−419) is known to be involved in direct interaction with the HIV-1 proteins Vpu and Nef. These two viral accessory proteins (Vpu and Nef) downregulate CD4 in HIV-1 infected cells by multiple strategies and make the body susceptible to all forms of infections. In this work, we carried out extensive replica exchange molecular dynamics simulations in explicit water with three popular protein force fields Amber ff99SB, Amber ff99SB*-ILDN, and CHARMM36 to characterize the equilibrium conformational ensemble of CD4-tail (402−419) and further validated the simulated ensembles with known NMR data. We found that ff99SB*-ILDN gives a better description of the structural ensemble of this peptide compared with ff99SB and CHARMM36. The peptide adopts multiple distinct conformations with varying degree of residual secondary structures. In particular, we observed 28, 7, and 5% average α-helical, β-strand, and 310-helix content, respectively, for ff99SB*-ILDN. The peptide chain shows the tendency of helix formation in a cooperative manner, seeding at residues 407−410, and subsequently extending toward both ends of the chain. Furthermore, we constructed Markov state model (MSM) from large-scale molecular dynamics simulations to study the dynamics of transitions between different metastable states explored by this peptide. The mean first passage times computed from MSM indicate rapid interconversion of these states, and the time scales of transitions range from several nanoseconds to hundreds of microseconds. Our results show good agreement with experimental data and could help to understand the key molecular mechanisms of T-cell activation and HIV-mediated receptor interference.

■

presented by MHC-II.7,8 Thus, by down-regulating CD4 in lymphocytes HIV makes the body susceptible to all forms of infections, as the immune system can no longer recognize antigens.9 In other words, HIV makes the body immunodeficient. CD4 is known to have 370 amino acids in its well-structured extracellular region,10 a transmembrane helix spanning 25 amino acids, and a cytoplasmic tail of 38 residues1,10−12 which plays an important role in both activating immune response and CD4 down-regulation. The membrane proximal residues 402− 419 are necessary for down-regulation of CD4 by the HIV accessory proteins Nef and Vpu.5,6,13−17 Furthermore, the presence of an α helix in this region was identified to be important for the binding of both Nef and Vpu.15−19 Kim et al.20 demonstrated that although Lck binds via Zn2+ at the CxCP motif (residues 420−423) in the cytoplasmic tail of

INTRODUCTION Human cluster determinant 4 (CD4) is a type I transmembrane glycoprotein and consists of three regions: an extracellular domain, a short transmembrane region, and a cytoplasmic region at the C-terminal end. CD4 is expressed on a subset of T-cells, primarily helper T cells.1 T-cells are a type of lymphocyte that is either CD4+ (alarm cells) or CD8+ (killer cells). By binding to major histocompatibility complex class II (MHC-II) molecules, CD4 serves as a coreceptor during antigen recognition by the T-cell receptor (TCR). The extracellular domain of CD4 interacts with the complex made of MHC-II, antigen-peptide and TCR, following which the cytoplasmic tail of CD4 binds noncovalently to the lymphocyte-specific protein kinase (Lck). This association is required for efficient antigen-induced T-cell activation and further TCR signaling.2 Besides its cellular functions, CD4 serves as a primary receptor for human immunodeficiency virus (HIV) type 1.3,4 Two regulatory proteins of HIV-1, viral protein U (Vpu) and negative factor (Nef), interact directly to the cytoplasmic part of CD4 and downregulates CD4 from the cell surface by multiple strategies.5,6 CD4 activates the immune response of our body upon recognition of the antigens © XXXX American Chemical Society

Special Issue: Biman Bagchi Festschrift Received: March 31, 2015 Revised: June 22, 2015

A

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

that are otherwise kinetically inaccessible to conventional MD simulations. Temperature replica exchange molecular dynamics (REMD) simulation, which is an enhanced sampling technique that efficiently samples the protein conformational space,42,43 is used to identify the structural ensemble of this peptide to alleviate the problems associated with poor sampling in conventional MD simulations. Furthermore, the transferability of the force fields is another major issue for studying the IDPs as the modern protein force fields are mainly optimized for well-folded globular proteins. The simulated ensemble thus generated can be biased toward the force field used to model the system, and different force fields may lead to very different results. It is, therefore, important to verify robustness of the simulation results obtained using different biomolecular force fields against the experimental data for the system under investigation. Here we have used three popular protein force fields, Amber ff99SB,44 Amber ff99SB*-ILDN, 45,46 and CHARMM36,47 to model the peptide in explicit water and validated the simulated ensembles with known solution NMR data. We have constructed Markov state model (MSM), a kinetic network model that can predict long time scale dynamics from many short simulations, to study the dynamics of this peptide. For this purpose, extensive unbiased MD simulations (for a total of ∼25 μs) in explicit water were carried out starting from different initial conformations obtained from REMD simulations. We found that the different metastable states identified for this peptide interconvert rapidly with mean first passage times (MFPTs) ranging from several nanoseconds to hundreds of microseconds. Our study reveals that the freeenergy landscape of the peptide is almost flat with transiently populated residual secondary structures. The findings of the present study could provide important insights into conformational preferences of this peptide in different physical environments or with different binding partners under physiological conditions. Furthermore, it could help us to understand the key molecular mechanisms of T-cell activation and HIV-mediated receptor interference.

CD4, residues 406 to 415 form a helix while residues beyond 423 are unstructured in the CD4-Lck-Zn2+ complex. The structural study carried out by Wray et al.10 using circular dichroism (CD) and NMR spectroscopy for the complete cytoplasmic domain of CD4 (containing 38 amino acids) reported that the peptide is unstructured in aqueous solution and forms a stable α-helix that extends from residues 402 to 417 only in the presence of 50% membrane mimetic trifluoroethanol (TFE). Similarly, Kim et al.20 found that in aqueous solution the isolated cytoplasmic tail of CD4 is largely unstructured, that is, intrinsically disordered. Moreover, NMR studies carried out for transmembrane and cytoplasmic domains of CD4 in the presence of dodecylphosphocholine (DPC) micelles found that membrane proximal residues 404− 413 in the cytoplasmic tail form an α helix that is reasonably stable even at 45 °C.21,22 In contrast, the structure of isolated CD4 peptide (residues 403−419) studied in aqueous solution by Willbold and Rosch23 reported that the peptide has the tendency of forming an α-helix between residues 403−412, although the average α-helical content was found to be ∼25% for this peptide. The structural plasticity offers a variety of functional advantages to intrinsically disordered proteins/peptides (IDPs), viz., flexibility to interact with different partners by the same peptide, specific but low affinity binding, and ability to modulate their functions by post-translational modifications.24 IDPs are therefore a central player in key cellular processes including cell signaling, regulation, and transcription and are also linked to several human diseases.25,26 Characterization of the conformational substates of the IDPs under physiological conditions is important to understand their functions at the molecular level.24,27−31 Standard biophysical techniques are, in general, of limited use for studying IDPs due to a lack of stable conformation, and further the secondary structures are weakly populated, transient, and confined to small fraction of the protein sequence. IDPs do not fold into a single free-energy minimum state unlike globular proteins; rather, they rapidly interconvert between an ensemble of conformational substates with high temporal fluctuations in the backbone dihedral angles and atomic positions.26,27,32 Despite the advancement of experimental methods such as NMR, electron spins resonance (ESR) and small-angle X-ray scattering (SAXS), characterizing the native state ensemble of IDPs remains a challenge. Molecular dynamics (MD) simulations are increasingly becoming a powerful tool to complement experiments by providing information on both structure and dynamics at atomic resolution.30,31,33−41 Previous experimental studies confirmed that the residues 402−419 of the cytoplasmic tail of CD4 are necessary and sufficient for the down-regulation of CD4 by Nef and Vpu.5,6,13−17 Despite its biological significance, our understanding at the molecular level of the binding mechanism and function of CD4 tail is very limited. It is essential to systematically characterize the conformational substates explored by this peptide in solution and determine the nature of the underlying free-energy landscape that can provide a microscopic description of its internal dynamics to have an in-depth understanding of its functioning mechanism. In the present study, we aim at studying the equilibrium structural ensemble of the peptide in aqueous solution using atomistic MD simulations; however, the complexity of the conformational space of natively disordered proteins/peptides demands high computational cost to sample functionally relevant states

■

METHODS System Preparation and REMD Simulations. The NMR structure from protein data bank (PDB ID 1WBR) was selected as the starting structure, and the missing N-terminal residue (Arg 402) was added with the help of pymol software (https:// www.pymol.org/). The resulting peptide was capped with acetyl group on N-terminus and methyl amide on C-terminus. The peptide Ace-RQAERMSQIKRLLSEKKT-NME was modeled using three popular atomistic force fields, Amber ff99SB,44 Amber ff99SB*-ILDN,45,46 and CHARMM3647 to check the force-field dependence of the structural ensemble generated from our simulations. All simulations in this study were performed using Gromacs 4.6.3 suite.48,49 The peptide was initially minimized for 10 000 steps using steepest descent algorithm; subsequently, the minimized conformation of the peptide was solvated in a dodecahedron box using TIP3P50 water model with 5158 water molecules. To neutralize the overall charge, 4 Cl− ions were added in the box. The minimum distance between the atoms of peptide and the box edge was 12 Å. The particle mesh Ewald method51 was used for long-range electrostatic interactions with a real space cutoff of 10 Å. LINCS52 method was used to constrain all of the bonds with H atoms, and periodic boundary conditions were employed in all of our simulations. The bonds and the angles of the water molecules were constrained using the SETTLE53 algorithm. B

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Equilibration phase consists of the following steps: first, the peptide was only allowed to relax by applying restraints on positions of all the solvent molecules and ions; subsequently, solvents were allowed to relax by restraining the positions of all the peptide heavy atoms. Finally, all of the restraints were removed, and the full system was equilibrated at 300 K using modified Berendsen thermostat54 (v-rescaling method) with a coupling constant of 0.5 ps for 2 ns. The system was further equilibrated in NPT ensemble at 300 K and 1 atm for 5 ns using Berendsen barostat55 with a coupling constant of 2 ps. The average density obtained at the end of the NPT simulation was used for all subsequent NVT simulations. The system was simulated in NVT ensemble at 300 K for another 2 ns so as to release the strain of the barostat. For each of the three force fields, a total of 64 replicas were used to carry out the REMD simulations in the NVT ensemble. The temperatures of different replicas were chosen so as to get an average acceptance ratio of ∼0.2.43,56 This resulted in following temperature distribution: 300, 302.5, 305, 307.5, 310, 312.5, 315, 318, 321, 324, 327, 330, 333, 336, 339, 342, 345, 348, 351, 354, 357, 360, 363.5, 367, 370.5, 374, 378, 382, 386, 390, 394, 398, 402, 406, 410, 414, 418, 422, 426, 430, 434, 438, 442, 446, 450, 455, 460, 465, 470, 475, 480, 485, 490, 495, 500, 505, 510, 515, 521, 527, 533, 539, 545, and 551. Each replica was initially equilibrated at the given temperature in the NVT environment for 5 ns. Subsequently, the simulations were run for 200 ns per replica with a time step of 2 fs, and an exchange between adjacent replicas was attempted at every 2 ps. The coordinates were saved at every 2 ps for data analysis. Dihedral Principal Component Analysis. Dihedral angles are the main degrees of freedom responsible for backbone flexibility and the formation of secondary structural elements. We performed dihedral principal component analysis (dPCA) of trajectories sampled by REMD to construct the free-energy landscape. In the dPCA method, Cartesian coordinates are replaced by dihedral angles; specifically it considers the dihedral angles (ϕi, ψi) of the peptide backbone, where i = 1, ..., Np and Np is the number of peptide groups. Instead of using the dihedral angles directly in the PCA and to avoid potential problems due to periodicity of dihedral angles, Mu et al.57 introduced the variables

energy of the conformation. To weigh each structure with its potential energy to get a proper average value, a variant of the weighted-histogram analysis method (WHAM) suitable for parallel tempering simulations, PTWHAM,58 was employed. MSM Construction and Validation. We divided all of the REMD conformations at 300 K into 500 clusters using the kcenters clustering algorithm59 and then randomly chose one conformation from each cluster as initial conformations for subsequent unbiased MD simulations. The length of each simulated trajectory was more than 40 ns, giving a total of ∼25 μs simulation time. All of the MD simulations were performed at 300 K with the same setup as used for REMD simulations. For the construction of MSM, simulation trajectories were saved every 2 ps. To provide a suitable basis for the discretization of the state space, we first reduced the dimensionality of the data using the time-structure-based independent component analysis (tICA) method.60 Each conformation was represented as a single vector whose elements were defined as the distances between all 18 Cα pairs plus the backbone dihedral angles of the CD4 peptide. This high-dimensional simulation data were then projected onto five dimensions (tICs) with a 2 ns lag time, which identified the slowest varying variables. The k-means clustering algorithm was used to discretize the data into 1000 microstates. We also checked the consistency of the results by increasing the number of clusters and projected dimensions. The transition between different states was described by the memoryless Master equation,61−66 Ṗ (t) = KP(t), where P(t) is a column vector of state probabilities at time t, and K is the transition rate matrix with off-diagonal elements Kji ≥ 0 (the rate of transitions from state i to state j) and diagonal elements Kii = −∑j≠iKji. The formal solution of this equation for a time step τ is given by P(t + τ) = exp(τK)P(t), and the transition probability matrix T(τ) is defined as T(τ) = exp(τK). To obtain T(τ), the number of transitions between the microstates at an interval of a certain lag time τ was counted and the transition count matrix was then symmetrized and normalized by column. We first identified the lag time τ at which the model is Markovian by observing the implied time scales at different lag times. At a specific lag time τ, the implied time scales can be calculated as −τ τk = ln[λk (τ )]

x4i − 3 = cos(ϕi), x4i − 2 = sin(ϕi)

where τk is the relaxation time scale of the kth eigenvector of T(τ), τ is the lag time, and λk is the kth eigenvalue of T(τ). If the model is Markovian at a certain lag time, the relaxation time scales should not change when using longer lag times. As the lag time τ defines the time resolution of the model, the minimal lag time at which the model is Markovian was chosen to build the MSM. This MSM was then used to calculate the equilibrium state populations and the kinetic properties. To identify the metastable states in CD4 peptide, we further lumped the 1000 microstates that are kinetically close into 100 macrostates using the PCCA+ lumping algorithm.67,68 The MSM construction was performed using MSMBuilder software.69−71 Mean First Passage Time. MFPT is defined as the average transition time going from the initial state to the final state.72 The MFPT from the initial state i to the final state f (MFPTif) was determined as

x4i − 1 = cos(ψi ), x4i = sin(ψi )

The covariance matrix C is constructed using sines and cosines of (ϕ, ψ) protein backbone dihedral angles. By diagonalizing the covariance matrix, one obtains a set of orthogonal principal directions or components (dPCs) and corresponding eigenvalues. After finding the PCs, we used projections of the trajectory along two main dPCs to map out the free-energy landscape in these collective coordinates. The free-energy surface in terms of two dPCs was calculated as ΔF = − kBT[ln P(dPC1, dPC2) − ln Pmax(dPC1, dPC2)]

where P(dPC1,dPC2) is the probability distribution and Pmax(dPC1, dPC2) is subtracted to ensure that ΔF = 0 for the lowest free-energy minimum. We performed root-meansquare deviation (RMSD)-based clustering to find the representative structures for each dominant local minimum. The average property calculated depends on the probability of each conformation. The probability is proportional to the Boltzmann factor, which, in turn, depends on the potential

MFPTif =

∑ T(τ)ij (τ + MFPTjf ) j

C

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B where τ is the lag time and T(τ) is the corresponding transition probability matrix, MFPTjf is the mean first passage time of the state j to final state f, and the boundary condition is MFPTff = 0.

Before using the trajectories for analysis, it is important to make sure that the simulations are converged. By evaluating the fluctuations in certain structural properties with time, one could verify whether the trajectories represent the equilibrium state of the system. We calculated the fluctuations in root-mean-square deviations (RMSDs) and the radius of gyration (Rg) for trajectories at different temperatures from 300−551 K. Figure S3 in the SI shows the RMSD plots as a function of time, where the reference structure is taken as the starting structure obtained from PDB, and Figure S4 in the SI shows Rg as a function of time. It can be clearly seen that both RMSD and Rg fluctuate randomly about their equilibrium values, and no net drifts are observed in these quantities throughout the simulated trajectories. We also observed an increase in average Rg values as the temperature is increased (Figure S5 in the SI), indicating more expanded structures of the peptide at higher temperatures. Furthermore, the convergence of REMD trajectories was verified by calculating the potential of mean force (PMF) along Rg and backbone RMSD (using PTWHAM) in two different time intervals 0−150 and 0−200 ns. As shown in Figures S6 and S7 in the SI, there is little difference in the PMF of both Rg and RMSD at 300 K, suggesting that an adequate sampling of the conformational space of the peptide has indeed been achieved, and the simulation data can reliably be used to compute the equilibrium properties of the system. We therefore decided to use the full-length trajectories for all further analysis in this study. REMD Simulations Reproduced the NMR Chemical Shifts and Scalar Coupling Constant 3JHNHα. To evaluate the ensembles produced by different force fields, we first examine the Hα and HN chemical shifts of each amino acid residue in the CD4 tail, for which solution NMR data are available.23 Furthermore, to avoid the dependency of results on the engine used for prediction of chemical shift, the Hα and HN chemical shifts of each amino acid residue were calculated for all of the simulated conformations at 300 K using the

■

RESULTS AND DISCUSSION Sampling of Phase Space using REMD and Convergence of Trajectories. The starting configuration for each Table 1. Uncertainty (χ2δ) between Experimental and Calculated Chemical Shift and Root-Mean-Square Difference (RMSD) between Experimental and Calculated Average 3JHNHα Values for Different Force Fieldsa SHIFTX

a

χ2δ

SPARTA+ χ2δ

force fields

(HN)

(HA)

χ2δ (HN)

χ2δ (HA)

RMSD (3JHNHα)

ff99SB ff99SB*-ILDN CHARMM36

0.0277 0.0985 0.1062

0.242 0.160 0.361

0.0293 0.0909 0.0712

0.275 0.062 0.109

1.085 0.754 0.786

Experimental values are taken from ref 23.

replica was indexed, and the replica index was traced with time at different temperatures for all three force fields to confirm proper sampling by each trajectory. To demonstrate the performance of the REMD algorithm, Figure S1 in the SI shows a good exchange among replicas for three representative temperatures. For example, the first trajectory (i.e., 300 K) has received conformations from the other 63 replicas, assuring random walk in replica space. Figure S2 in the SI displays potential energy distribution for all of the replicas; one can clearly see that a finely tuned overlap between adjacent replicas took place. The average exchange probability calculated for all pairs of adjacent replicas is ∼0.25, further confirming a good exchange among adjacent replicas.

Figure 1. Comparison between experimental23 and calculated Hα chemical shifts of CD4-tail at 300 K: (A) SHIFTX, (B) SPARTA+, and (C) comparison between experimental23 and calculated 3JHNHα for ff99SB*-ILDN. Standard error is shown by the error bars. D

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 2. Free-energy landscape of CD4-tail at 300 K obtained from the REMD simulation (with ff99SB*-ILDN force field) as a function of the first two principal components dPC1 and dPC2. The color scale shows the free energy in kilocalories per mole. The representative structures at the local minimum of free energy basins (a−j) are displayed on the periphery of the image. N represents the N-terminal of the peptide. Secondary structures are colored in purple (α helix), yellow (β strand), blue (310 helix), tan (β bridge), cyan (turn), and white (coil).

Figure 3. Free-energy landscape of CD4-tail at 300 K obtained from the REMD simulation (with ff99SB force field) as a function of the first two principal components dPC1 and dPC2. The color scale shows the free energy in kilocalories per mole. The representative structures at the local minimum of free-energy basins (a−e) are displayed on the periphery of the image. N represents the N-terminal of the peptide. Secondary structures are colored in purple (α helix), yellow (β strand), blue (310 helix), tan (β bridge), cyan (turn), and white (coil).

program SHIFTX73 and SPARTA+74 (Tables S1 and S2 in the SI). In recent studies, Ball et al.34,75 introduced a measure χ2δ to quantify the agreement between the chemical shifts obtained from experiment and calculated from the simulated ensemble, defined as χδ2

1 = N

N

∑ i=1

uncertainty (RMSD, from experiment). Reported uncertainties for the SHIFTX calculator are σδ = 0.23 ppm for Hα and σδ = 0.49 ppm for HN and for SPARTA+ calculator σδ = 0.25 ppm for Hα and σδ = 0.49 ppm for HN. Calculated χ2δ values (Table 1) for each of the three force fields show very good agreement between experimental and theoretical proton chemical shifts. Thus, chemical shifts data alone are not sufficient to differentiate the structural ensembles generated using these three force fields. To evaluate which ensemble is more close to experimental data, we further calculated three-bond scalar coupling constant, 3JHNHα (Table 1 and Table S3), using the Karplus equation76

(δi ,cal − δi ,exp)2 σδ2

where δi,exp and δi,cal are the experimental chemical shift and calculated average chemical shift of residue i, respectively; N is the number of residues in the peptide; and σδ is the calculator E

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 4. Free-energy landscape of CD4-tail at 300 K obtained from the REMD simulation (with CHARMM36 force field) as a function of the first two principal components dPC1 and dPC2. The color scale shows the free energy in kilocalories per mole. The representative structures at the local minimum of free-energy basins (a−g) are displayed on the periphery of the image. N represents the N-terminal of the peptide. Secondary structures are colored in purple (α helix), yellow (β strand), blue (310 helix), tan (β bridge), cyan (turn), and white (coil).

Figure 5. Distributions of secondary structure of CD4-tail (A) α helix, (B) β sheet, (C) turn, and (D) 310 helix content at 300 K based on DSSP criteria as a function of residue index. Blue represents ff99SB, green for ff99SB*-ILDN, and red for CHARMM36 force field.

three-bond scalar coupling constant (3JHNHα) and Hα and HN chemical shifts together clearly demonstrates that ff99SB*ILDN performed better compared with other force fields for this peptide (Table 1 and Figure 1). In fact, several other groups have already been reported that ff99SB*-ILDN better reproduces the data obtained from NMR experiments for IDPs.45,78 CD4-Tail is Intrinsically Disordered having Transient Secondary Structures in Solution. Principal Component Analysis. A protein composed of N atoms has a 3N dimensional conformational space. PCA is a useful technique

J(ϕ) = A cos2(ϕ − 60) + B cos(ϕ − 60) + C

where ϕ represents the protein backbone dihedral angle, with the parameter set A = 6.51, B = −1.76, and C = 1.60 originally suggested by Vuister and Bax.77 Table 1 shows the 3JHNHα RMSD values for all three force fields. Compared with ff99SB, the other two force fields (ff99SB*-ILDN and CHARMM36) show a reasonably good agreement between experimental23 and calculated 3JHNHα values from the simulated ensemble. (Note that, the reported RMSD of 3JHNHα with these parameter values of A, B, and C is 0.73 Hz.77). However, overall, the analysis of F

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 6. Probabilities of observing total n helical residues (A−C) and n consecutive helical residues (D−F) as a function of residue index for all three force fields throughout the simulations at 300 K. (A,D) ff99SB, (B,E) ff99SB*-ILDN, and (C,F) CHARMM36. Color shows fractional content of a particular type of helix.

Figure 7. Structural representations of the top eight highest populated clusters from cluster analysis of ff99SB*-ILDN trajectory at 300 K based on backbone RMSD of peptide. The structures are taken from the center of each cluster. N represents the N-terminal of the peptide. Secondary structures are colored in purple (α helix), yellow (β strand), blue (310 helix), tan (β bridge), cyan (turn), and white (coil).

to reduce the dimensionality of this complex space to study the collective motions of the protein.79,80 PCA diagonalizes the covariance matrix of the atomic coordinates. The eigenvectors define the directions of the motion, whereas the corresponding eigenvalues represent amplitudes. The eigenvalues are rank ordered in decreasing values; however, PCA in Cartesian coordinates often suffers from not being able to separate the internal and overall motions of the protein. PCA based on (ϕ, ψ) backbone dihedral angles of a protein evolved as a useful alternative technique to overcome this problem.81,82 Here we employed dPCA to disentangle interesting conformational modes from the structural ensemble sampled at 300 K for this peptide. The contribution of first few PCs to the internal

motion of the peptide is shown in Table S4 in the SI. Contributions from the first two slowest modes are ∼20, ∼33, and ∼24% for ff99SB, ff99SB*-ILDN, and CHARMM36, respectively. Two-dimensional free-energy surfaces along the first two PC modes corresponding to largest eigenvalues are shown in Figure 2 for ff99SB*-ILDN, Figure 3 for ff99SB, and Figure 4 for CHARMM36. The energy landscape of globular proteins is generally characterized by a funnel-like shape with a well-defined free-energy minimum representing an ensemble of conformations rather than a single structure;83−86 however, we found that the free-energy landscapes of the CD4 peptide are almost flat with several energy minima of varying depths (Figures 2−4). We clustered the structures populated in G

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 8. Ramachandran free-energy surface for the 16 residues 403−418 (a-p) from the trajectory of ff99SB*-ILDN at 300 K.

Figure 10. Equilibrium populations of the 100 macrostates. Figure 9. (A) Implied time scales as a function of lag time for the 1000 microstate MSM. (B) Implied time scales as a function of lag time for the 100 macrostate MSM. In each case, relaxation of the 10 slowest eigenvectors of the transition probability matrix is shown.

smaller length (6−9 residues) and located mostly in the central part of the peptide. Structures in the other two local energy basins “d” and “e”, however, are enriched in random coils with partial α-helical structures. Free-energy surfaces derived for two other force fields (Figures 3 and 4) are relatively more shallow in nature with many major and minor basins compared with ff99SB*-ILDN. The free-energy landscape for CHARMM36 contains seven dominant basins, in which three of them, “a”, “c”, and “g”, represent random coil structures with turns, whereas basin “b” corresponds to β-hairpin structures; the other three basins “d”, “e”, and “f” consist of α-helical structures having different lengths; however, ff99SB force field shows a highly flat energy surface with several low-energy basins, where most of them represent either β-sheet structures or random coil structures, except the basin “e” that contain α-helical structures. Overall, the feature of the free-energy landscapes reveals the presence of distinct conformations having low free-energy

dominated basins, and the representative structures for each basin are shown in Figures 2−4. All three force fields sampled transient secondary structures such as α-helix, β-sheet, and 310helix with large population of random coils and turns. Out of five major low-energy basins that were identified for ff99SB*ILDN, named as “a”, “b”, “c”, “d”, and “e” (Figure 2), the most dominant basins are “a” and “c”; “a” consists of structures having largest α-helical region of length 12−16 residues near the N-terminal region, whereas “c” represents β-hairpin structures. Basin “b” that is close to basin “a” and separated by small energy barrier (∼kBT) is also populated by α-helical structures, but compared with basin “a” these helices are of H

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 11. Transition network analysis for the top 10 most populated metastable states of CD4 peptide. Each metastable state is represented by a circle, and the size of the circle is proportional to its corresponding equilibrium populations. The transition probabilities between metastable states are represented by the thickness of the lines connecting them; thicker line implies a higher transition probability. Dotted lines have been used to connect metastable state “S2” with states “S3”, “S5”, and “S6” because these transitions occur through other small-populated states (not shown). Representative conformations and the associated equilibrium populations of these metastable states are also shown.

is an algorithm that calculates the energy of the intrahydrogen bonds in protein backbone based on the distance between CO and NH groups of different residues in the sequence. It defines a hydrogen bond if the bond energy is less than −0.5 kcal/mol. It assigns α helix (H), 310 helix (G), π helix (I), β bridge (B), extended strand (E), turn (T), curvature (S), and loop (blank). In the case of overlap between two secondary structures for a residue, priority is given in the order of H, B, E, G, I, T, and S. We compared the secondary structure propensities (in particular, α helix, β sheet, 310 helix, and turn) of all residues for each of the three force fields in Figure 5. CHARMM36 shows that N-terminal residues have more propensities to form helix, whereas residues in the central portion have more helical tendency in the case of ff99SB*-ILDN. On the contrary, ff99SB shows a homogeneous distribution of low α-helical content. It should be noted that ff99SB*-ILDN shows greater helical propensities throughout the length of the peptide compared with the other two force fields. The average secondary structure content calculated for ff99SB, ff99SB*-ILDN, and CHARMM36 is found to be 8, 28, and 11% for α helix; 5, 5, and 1% for 310 helix; 12, 7, and 4% for β sheet; and 18, 19, and 10% for turn, respectively. Previous studies10,23 with the help of circular dichroism spectroscopy did not find any β strand in aqueous solution, but in membrane mimetic environment,

barriers, suggesting that the structural ensemble of CD4 tail is quite heterogeneous where the conformations can interconvert rapidly. Therefore, on the basis of the dPCA analysis it is clearly evident that this peptide is intrinsically disordered, having preferences to form some residual α-helical structure in aqueous solution under physiological conditions. Secondary Structure Analysis. Solution structure of CD4 peptides 403−41923 and 396−43310,20 of the cytoplasmic domain of CD4 has been previously studied by NMR and CD spectroscopy. Willbold and Rösch23 reported that this peptide (residues 403−419) has the tendency to form an α helix between Gln403 to Arg412 in aqueous solution. The estimated average α-helical content of the peptide was ∼25%, indicating that only a subset of conformations forms the helical structure.23 On the contrary, two independent studies10,20 carried out considering the entire cytoplasmic tail of CD4 (with 38 residues) reported it as an unstructured peptide under aqueous conditions. N-terminal regions of the cytoplasmic tail are proximal to the membrane, and previous studies21,22 have shown that membrane-mimicking environment induces helicity in this region ranging from residues 404−413. We computed the secondary structure content of CD4-tail (at 300 K) for each residue using the computer program DSSP87 (Dictionary of Secondary Structure of Proteins). DSSP I

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

starting from 403 to 418 (terminal residues exhibit no helical propensity throughout the simulation) when using ff99SB*ILDN (Figure 6E). Interestingly, in the case of ff99SB*-ILDN, the peptide chain shows tendency of helix formation in a cooperative manner, seeding at residues 407−410 and subsequently extending toward both ends of the chain (Figure 6E). We found that CHARMM36 has helix-forming tendency mostly toward the N-terminal; however, the number of such conformations is less as compared with ff99SB*-ILDN (Figure 6E,F). Clustering of Structural Ensemble. To identify the most populated clusters in the REMD ensemble (at 300 K), the structure-based clustering algorithm89 was employed. The conformations are sorted into clusters based on backbone root-mean-square distance (RMSD). The residues 404−414 or a subset of these residues are reported in literature to form a stable helix in membrane mimicking environment. We considered only these residues with a RMSD cutoff of 1 Å for clustering. Central structures of resulting top eight clusters are shown in Figure 7 for ff99SB*-ILDN, Figure S8 in the SI for ff99SB, and Figure S9 in the SI for CHARMM36. The populations of most populated clusters from ff99SB, ff99SB*ILDN, and CHARMM36 are 5, 13.5, and 11.8%, respectively; even the top eight clusters together represent ∼27, ∼43, and ∼40% of total structures in the ensembles, respectively. This indicates that the simulated ensembles consist of highly diverse conformations. Cluster 1 and cluster 3 from ff99SB*-ILDNsimulated ensemble represent α-helical structures, whereas cluster 2 consists of β-hairpin structures; other clusters (clusters 4−8) have mostly random-coil conformation with some residual secondary structures. Although the top eight clusters from CHARMM36 have α-helix and β-hairpin structures, their contributions are less compared with ff99SB*-ILDN (Figure S9 in the SI). On the contrary, the clusters from ff99SB are more or less uniformly populated, where the topmost cluster (cluster 1) contributed only 5%. Moreover, they are mostly dominated by β-hairpin structures; only cluster 8 (with a population of 2.5%) shows α-helical structures (Figure S8 in the SI). In fact, these results are consistent with our findings of secondary structure analysis and dPCA discussed in previous sections. Taken together, structural ensembles of this peptide show a very heterogeneous distribution of conformations, a typical characteristic feature of an IDP. Ramachandran Maps. The backbone flexibility of a peptide is measured by the regions of the Ramachandran space (ϕ−ψ space) that the peptide explore. The dihedral angles ϕ (between N−Cα) and ψ (between Cα-C′) have restricted space in the Ramachandran plot due to steric clashes. In addition, each secondary structure has a fingerprint region in the Ramachandran plot: α ∈ −160° < ϕ < −20° and −120° < ψ < 50°; β ∈ −180° < ϕ < −90° and 50° < ψ < 240° or 160° < ϕ < 180° and 110° < ψ < 180°; ppII ∈ −90° < ϕ < −20° and 50° < ψ < 240°; others, everything else.90 The local conformational sampling of residues can be represented by the free-energy maps of the Ramachandran plots. The freeenergy maps at 300 K for all residues (except terminal residues, 402 and 419) are shown in Figure 8 for ff99SB*-ILDN, Figure S10 in the SI for ff99SB, and Figure S11 in the SI for CHARMM36. The results of the free-energy maps for residues 403−418 that form the peptide show that these residues explore all of the allowed regions in the Ramachandran space with finite population of different secondary structures. This indicates that the statistical coil model91 could be a good

Figure 12. MFPTs between different metastable states. (A) MFPTs between top 10 macrostates and (B) distribution of MFPTs between all pairs of macrostates.

Wittlich et al.22 reported 12−15% β-strand conformations for the entire cytoplasmic tail of CD4. We observed some β content in our simulated ensembles from all three force fields. It has been reported that although ff99SB successfully produces values of thermodynamic and structural properties close to the experimental observations, the propensity of α-helix is somewhat underestimated.88 Our results also demonstrate that ff99SB underestimates the propensity to form α helix. The calculated average α-helical content for ff99SB*-ILDN is 28%, which is in excellent agreement with the experimental value of ∼25%.23 Thus, our secondary structure analysis results also suggest that ff99SB*-ILDN performs better compared with ff99SB and CHARMM36. We further studied the distribution of lengths of helix during simulations. Figure 6 shows populations of observing total n helical residues (A, B, and C) and n consecutive helical residues (D, E, and F) throughout the simulations at 300 K for all three force fields. It is clear that the population of having more than four helical residues in the peptide is very low for ff99SB (Figure 6A). In addition, the helix length calculated with respect to peptide sequence shows that ff99SB does not have any preference for a particular region of this peptide; the longest helix observed for ff99SB was 12 residues long (i.e., helix with 3 turns) starting from residue 406 to 418 (Figure 6D). Figure 6C displays that finding of total number of helical residues (even up to 12) is considerably more for CHARMM36 compared with ff99SB, whereas the peptide shows a significant population of total number of helical residues (up to 14) in the case of ff99SB*-ILDN (Figure 6B). Furthermore, our results show that the longest helix formed using CHARMM36 is 14 residues long starting from 405 to 418 (Figure 6F), whereas the largest possible helix found for this peptide is 16 residues long J

DOI: 10.1021/acs.jpcb.5b03092 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

“S2” to “S1”, “S4”, and “S7”−“S10” occur via “S3”, “S5”, or “S6” as intermediate states (Table S6 in the SI). Native states (most populated states) of well-folded proteins act as kinetic hubs, which means they are kinetically highly accessible and mediate “non-native to non-native” (one unfolded state to another unfolded state) transition by first folding and then unfolding rather than transitioning directly between two unfolded states.93−95 The most populated metastable states of the CD4 peptide do not show hub-like behavior, and transitions between different unfolded states do not prefer to be mediated via the folded states (Table S6 in the SI). To study the kinetics of transition between the metastable states, we calculated the MFPTs from microstate MSM. The MFPTs thus computed indicate that the time scales of transitions between different metastable states range from several nanoseconds to hundreds of microseconds. Figure 12A shows the MFPTs between the top 10 metastable states, and the distribution of MFPTs between all pairs of metastable states is plotted in Figure 12B. The fastest transition occurs at ∼48 ns and the slowest occurs at ∼298 μs. The fastest transition corresponds to the conversion from the state “S9” to the state “S1”, while the slowest transition is from the state “S2” to the state “S99”. From the pairwise MFPTs analysis, we found that the slowest transitions correspond to exchange between the βsheet structures and other metastable states. This is because it involves breaking of hydrogen bonds, reorganizing the conformation, and then forming the new contacts. MFPTs between other metastable states, except those involving the state “S2” are relatively faster (Figure 12A) and are in line with the relatively low energy barriers between free-energy minima observed in dPCA (Figure 2).

description for the CD4-tail (402−419) in aqueous solution. The calculated fractions from simulated ensembles of ff99SB, ff99SB*-ILDN, and CHARMM36 are 0.32, 0.54, and 0.22 for α, 0.28, 0.17, and 0.26 for β, and 0.24, 0.14, and 0.34 for ppII (polyproline-II), respectively. Table S5 in the SI shows the fractions of calculated secondary structures for each residue of this peptide. In the case of ff99SB*-ILDN, residues 403−416 explore frequently α regions during simulation, well reflected in the values of fractional content ≥0.5, whereas the probability to be in α regions of residues is small (