Stability of trimeric DENV envelope protein at low and neutral pH: An insight from MD study

Biochimica et Biophysica Acta 1834 (2013) 53–64

Contents lists available at SciVerse ScienceDirect

Biochimica et Biophysica Acta journal homepage: www.elsevier.com/locate/bbapap

Stability of trimeric DENV envelope protein at low and neutral pH: An insight from MD study Kshatresh Dutta Dubey, Amit Kumar Chaubey, Rajendra Prasad Ojha ⁎ Biophysics Unit, Department of Physics, DDU Gorakhpur University, 273009, India

a r t i c l e

i n f o

Article history: Received 27 March 2012 Received in revised form 13 August 2012 Accepted 14 August 2012 Available online 23 August 2012 Keywords: Trimeric interaction MM–PBSA Thermodynamical parameter B-factor Ionic concentration Dengue envelope protein

a b s t r a c t Change in pH plays a crucial role in the stability and function of the dengue envelope (DENV) protein during conformational transition from dimeric (pre-fusion state) to trimeric form (post-fusion state). In the present study we have performed various molecular dynamics (MD) simulations of the trimeric DENV protein at different pH and ionic concentrations. We have used total binding energy to justify the stability of the complex using the MMPBSA method. We found a remarkable increase in the stability of the complex at neutral pH (pH ~ 7) due to the increment of sodium ions. However, at very low pH (pH ~ 4), the total energy of the complex becomes high enough to destabilize the complex. At a specific pH, almost at a range of 6, the stability of the complex is significantly better than the stability of the trimer at neutral pH, which connotes that the trimer is most stable at this pH (pH ~ 6). © 2012 Elsevier B.V. All rights reserved.

1. Introduction Dengue is a mosquito-borne viral disease that has become a major public health concern worldwide in recent years [1]. However, there is no vaccine or treatment other than vector control and supportive medical care. Three structural proteins (C, M, and E) and a lipid bilayer package the positive-strand RNA genomes of flaviviruses [2]. The core nucleocapsid protein, C, binds directly to genomic RNA, while the major envelope glycoprotein, E, and the membrane protein, M, form the outer protein shell [3]. Membrane fusion is the primary molecular event during the entry of the dengue virus into the host cell. The envelope protein E mediates the fusion process in a pH dependent manner. Therefore, envelope protein is the most fertile field for inhibition of dengue infections. The structure of the dengue envelope protein has been determined in pre- and post‐fusion separately [4–6]. Every unit of envelope proteins has three domains, based largely on beta sheets. One of these domains, (domain II) an elongated finger like structure, bears a loop at its tip with a hydrophobic sequence which mediated the fusion process. More details of the structure of envelope proteins have been already stated elsewhere [7–12], which reveal a homo-dimeric structure before fusion and a homo-trimeric structure after fusion.

⁎ Corresponding author. Tel.: +91 551 2202167. E-mail address: [email protected] (R.P. Ojha). 1570-9639/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.bbapap.2012.08.014

Since, change in pH plays a crucial role in the transition from pre-fusion to post fusion state, it is necessary to study the effect of pH on an envelope protein. Although X-ray crystallographic study has utmost importance in the structural information of envelope proteins we have to rely upon dynamical methods, say MD simulations, to explore the conformational mechanism of these proteins. For several years now, computational techniques are employed as an effective tool for approaching biochemical themes, such as protein– protein interactions, protein–drug interactions, drug discovery, and binding site search, among others [13–17]. In order to obtain an accurate value of binding affinity of such biomolecular complexes, free energy methods are readily used. Free energy perturbation (FEP) and thermodynamic integration (TI) are rigorously used for reliable prediction of binding affinity. However, these methods are extremely time consuming [18], therefore, MMPBSA [19] and LIE [20] methods are widely used. These methods make use of potential energy terms and provide detailed insight into the structure–activity relationship whose detail may be found elsewhere [21–26]. Several applications of in-silico drug designing for dengue infections have been applied [27–29]. However, some molecular dynamics simulations have also been performed to study the dynamical characteristics of envelope proteins [30–32], but these studies were performed on a monomeric complex. Since an envelope protein exhibits a multimeric structure (dimer before fusion and trimer after fusion), therefore MD simulations of a multimeric complex will be more interesting. In our previous publication [33], we have investigated the protein–protein interactions between dimeric units of the dengue envelope protein

54

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

at normal and low pH using MD simulations, separately. We have shown that the interactions between dimeric units decrease significantly with a decrease of pH. However, a complete study of the role of pH on envelope proteins demands a molecular dynamical study of a trimeric complex in different pH and ionic concentrations. In the present study, we have performed molecular dynamics simulations at different pH (3, 5 and 7) and different ionic strengths for the dengue envelope protein in a post-fusion trimeric form to study the thermo-dynamical parameters in the stability of this complex, for the first time. 2. Method The initial structure was imported from the RCSB protein data bank (www.pdb.org) [PDB ID:1OK8] [4]. The structure was available as a biological assembly of three units of envelope protein (as shown in Supplementary Fig. 1Sa). Each unit of envelope protein was composed of 394 amino acids. The imported structure has missing hydrogen atoms and some heavy atoms like nitrogen and oxygen which were added by the Leap module of Amber [34] using the latest Amber force field (ff10) parameters [35–37]. We prepared five different systems at different pH and different numbers of ions. For low pH (pH ~ 4), we protonated all histidines, aspartic acids and glutamic acids, while for pH ≈ 6, we protonated histidines only [38–41]. We modified all HIE (histidine protonation state at neutral pH according to Amber library) to HIP, and GLU was modified to GLH. Three different models were generated for simulations at neutral pH conditions. For the first model, we used crystallographic chloride ion and no additional chloride ions were added (shown in Supplementary Fig. 1Sb); for the second model, 5 chloride ions were inserted while for the third model 18 chloride ions were added. To compensate for the charge of each system, we added 4, 8 and 21 sodium ions, respectively. For low pH, we prepared two different systems at pH ~ 4 and pH ~ 6, separately. We added 170 and 29 chloride ions to counterbalance the positive charges of the system at pH ~ 4 and pH ~ 6, respectively. After model generation, all systems were neutralized separately with Na+ ions. All models were immersed in a truncated octahedral shell of TIP3P water [42]. A time step of 2 fs and a direct-space nonbonded cutoff of 12 Å were used. After the protein preparation, all systems were minimized to remove the steric clashes. The systems were then gradually heated from 10 to 300 K over a period of 200 ps and then maintained in the isothermal–isobaric ensemble (NPT) at a target temperature of 300 K and a target pressure of 1 bar using a Langevin thermostat [43,44] and a Berendsen barostat with a collision frequency of 2 ps and a pressure relaxation time of 1 ps, respectively. The hydrogen bonds were constrained using SHAKE [45]. We have used the velocity-Verlet algorithm (default algorithm for the Amber MD package) for MD simulations. Particle mesh Ewald (PME) was used to treat long range electrostatic using default parameters [46]. After getting the systems at our target temperature and pressure of 300 K and 1 bar, respectively, the simulations were continued up to 3 ns for equilibration with the same parameters as described above, for all five models separately. For the analysis of the systems, production dynamics was initiated and continued for another 14 ns and maintained in the isothermal–isobaric ensemble at the target temperature of 300 K and target pressure of 1 bar using the same Langevin thermostat and Berendsen barostat as in the equilibration process. The structures in the trajectories were collected at 10-ps intervals. The analysis of trajectories was performed with the Ptraj module of Amber11. We have used VMD 1.6.7 [47], Chimera-1.3 [48] and Maestro [49] graphical programs for visualization purposes. For the binding free energy calculations, we used the standard MM–PB/SA method [19,50]. Before the MM–PBSA analysis, all water molecules and the sodium ions were excluded from the trajectory. The dielectric constant used for the solute and surrounding solvent

was 1 and 80, respectively. During the analysis of the MM–PBSA trajectory, snapshots were gathered at the interval of 100 ps between the 4th and 12th nanoseconds. Similar applications of these methods can be found in our earlier publications [51–54]. 3. Results Every unit of trimer in post-fusion and pre-fusion structures has an equal number of amino acid residues. The fusion peptide at domain II in post-fusion has almost the same conformation but different orientations relative to pre-fusion. Three units of trimer form an extended cavity that runs along the three fold axis. In the crystal structure, a chloride ion is found at the tip of domain II of the trimeric protein, which is ligated by amide nitrogen of Lys 110 of each unit. We performed a molecular dynamics study of a trimeric envelope protein at normal and low pH to account the conformational flexibility during simulations. The solvent accessible surface area for all five complexes is reported in Table 1. The details of SASA for every complex at different snapshots can be found in the supplementary materials. From the table, we see that the surface area for the trimeric complex at pH ~ 4 is maximum, while it is least when the complex is at neutral pH with 5 chloride ions. The surface area at pH ~ 6 is greater than the surface area at neutral pH. It reveals that the trimeric complex is more solvent exposed at low pH. This result is also validated by an experimental study where it was found that the surface buried per monomer at neutral pH (before fusion) is almost half of the surface buried in a trimer (post-fusion) [5]. 3.1. MD at neutral pH After visualization of the trajectories (Video VS1 available in supplementary data) from production dynamics we see that there are no major changes in the conformation of the complex at neutral pH with reference to the initial structure. However, a large fluctuation in a peptide loop, oriented towards domain III is observed. Fig. 1 shows a single unit of the trimeric complex with a flexible loop colored in green. This loop corresponds to residues 145–155 of the first unit (similarly residues 539–549 of the second unit and residues 933–943 of unit III). The movement of ions within the trimeric complex is shown in Fig. 2. From the figure, we see that Na + 1185 (green lines) shows the least and almost constant deviation during the entire simulation. It interacts electrostatically with δ oxygen (OD2) of Asp 192 of the first unit and Asp 582 of the second unit and becomes a life partner of the complex. Na + 1185 creates a link that joins both units of the trimer by strong electrostatic interactions. These interactions, probably, give a stability to the trimeric complex (binding of Na + into the cavity of the trimer is shown in Supplementary Fig. 2S). The RMS (root mean square) deviations and thermal fluctuations (B-factor) are shown in Fig. 3a and b. We observe that the thermal fluctuations of all three units are almost symmetrical. The B-factor

Table 1 Mean solvent accessible surface area (SASA) at different pH and ionic concentrations. Complex

Mean SASA

System System System System System

66,253.77 65,807.13 66,317.83 68,338.03 66,601.21

1 2 3 4 5

The unit of area is Å2. Here, system 1 = normal pH and 1 chloride ion, system 2 = normal pH and 5 chloride ions, system 3 = normal pH and 18 chloride ions, system 4 = low pH (pH~4) and system 5 = low pH (pH~6). SASA for every system at different snapshots can be found in supplementary materials.

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

55

Fig. 1. Structure of a single unit of the trimeric complex. The loop that shows much flexibility during simulations is shown in green color.

of residues 145–155 of the first unit (similar for other units) is largest among other residues. It corresponds to a loop whose flexibility was noticed during visualization of the trajectory. This loop is composed of Gly, Lys, Ser, Thr, Asp and Glu. The structures from the trajectories reveal solvent exposure of these residues and there is no other possibility of main chain hindrance; which reduces stiffness of the loop. The RMS deviation shows that unit I and unit III have small deviations while the RMS deviation of unit II is slightly greater than units I and III. From Fig. 2 (and also video data in supplementary materials), we see that Na+ binds with units I and III that might be giving the stability to these units, hence the stability of units I and III is slightly greater than unit II. 3.2. MM–PBSA calculations We performed an MM–PBSA study to investigate the stability of the trimeric complex at neutral pH. From Table 2, we see that electrostatic energy (−21,528.35 kcal/mol) strongly favors the stability of

the complex. At neutral pH the complex forms several salt bridges with all three units of the trimer, which give large electrostatic interactions that stabilize the complex. The table showing the details of the salt bridges, is reported in Table 1S of supplementary materials. The van der Waals interaction energy is positive (8670.97 kcal/mol) which does not signify a favorable contribution in the stability of the complex. The total internal energy of the system is very high. The internal energy of a system is an implication of the kinetic energy or the molecular motion of the system, hence, the complex will have high molecular agitation during the simulation which does not favor the stability of the complex. Furthermore, we also notice that nonpolar contribution is positive (403.45 kcal/mol) which does not favor the stability of the trimeric complex. However, polar energy signifies a favorable interaction and it dominates the non-polar interactions such that total solvation becomes significantly negative which favors the formation of the trimeric complex. During simulations we found several water molecules forming hydrogen bonds with different parts of the complex as clear from the trajectories. These interactions provide stability to the complex and therefore, a favorable solvation contribution is observed during the MM–PBSA study. EELE + POL and EVDW + NONPOL are the compensation of electrostatic energy with polar contribution and van der Waals terms with the non-polar energy term, respectively. Furthermore, the total energy ETOT is also negative, which energetically favors the stability of the complex.

3.3. MD simulations of the trimer in low concentration (5 Cl− and 8 Na +) of ions

Fig. 2. The motion of ions in the trimeric complex at neutral pH. Different ions are shown in different colors which are indicated in the legend box. Ions closer to the complex, are shown here. Here, the unit of displacement on the y-axis is Å2 (mean square displacement).

We have performed molecular dynamics simulations in the presence of 5 chloride ions and 8 Na+ ions to investigate their role on the stability of the trimeric complex. The side and depth views of the initial structure are shown in Supplementary Fig. 3S. We see that chloride ions are arranged in the channel formed by the units of the trimer. Na+ ions are distributed outside of the trimeric complex. The root mean square deviations and thermal fluctuation are shown in Fig. 4a and b. It shows the stability of the trimeric complex in the context of RMS deviations and B‐factors. We see that the root mean square deviation of unit I and unit II is significantly small in the presence of Na+ and chloride ions while a large deviation in unit III is observed. The B-factor is shown in Fig. 4b. We notice

56

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

Fig. 3. (a) Root mean square deviation and (b) thermal fluctuation (B-factor) in Å2 residue-wise during simulation at neutral pH. The RMSD of the first, second and third units are shown in black, red and green colors, respectively. The first unit (residues 1–394) is shown at the top, the second unit (395–788) is shown in the middle while the third unit (789–1182) is shown at the bottom of the figure.

Table 2 Various energy contributions in the stability of the trimeric complex at neutral pH as calculated by MM–PBSA method. All energy parameters are calculated in kcal/mol. Here, EELE: Coulombic energy; EVDW: van der Waals energy; EINT = internal energy, EGAS = EELE + EVDW + EINT; ENONPOL: non polar solvation free energy; EPOL: polar solvation free energy; ESOL = ENONPOL + EPOL; EELE +POL = sum of electrostatic and polar energy terms; EVDW+NONPOL = sum of van der Waal energy and non-polar contributions; ETOT = EGAS + EPBSOL. Contributions

Mean value

Std dev

EELE EVDW EINT EGAS ENONPOL EPOL EELE+POL EVDW+NONPOL ESOL ETOT

−21,528.35 8670.97 21,775.49 8918.10 403.45 −12,408.13 −33,936.48 9074.42 −12,004.68 −3086.58

9556.15 3943.84 9665.05 4057.64 179.07 5508.83 7532.49 5329.80 1629.73

a sudden rise in the thermal fluctuation for unit III, which reveals instability in the trimer. From the trajectories (Video data VS2 in supplementary materials) we found three sodium ions, Na + 1188, Na + 1190 and Na + 1191, initially bound in the cavity formed by three units of trimer. The motion of these ions during simulations is shown by the diffusion plot in Fig. 5. We notice that Na + 1188, 1190 and 1191 bind with the complex during 10 ns of simulation because their diffusion in reference of the trimeric complex is very small. However, Na + 1188 has the highest resident time and it may be treated as a life partner of the trimeric complex. The binding of these ions gives stability to the complex at the initial stage of the simulations, shown in RMS deviations. As the simulation proceeds, the binding of Na + 1190 and 1191 fades away. However, Na + 1188, lives permanently with the complex. Since, Na + 1188 binds with the aspartate of units I and II and becomes a life partner of units I and II, these two units exhibit a stable fluctuation during the entire simulations. Due to exclusion of two ions, Na + 1190 and Na + 1191 from the cavity, the stability of the complex decreases. We monitored the stability of the complex thermodynamically using the MM–PBSA method as well. The various energy contributions

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

57

Fig. 4. (a) The root mean square fluctuations and (b) B-factor during the entire simulations in the presence of five chloride ions at neutral pH. The unit of the B-factor is Å2. All units have the same color and representations as Fig. 3.

during stability calculation are shown in Table 3. We see that the energy contributions due to electrostatic interactions are highly negative, which strongly favors the stability of the trimeric complex. However, van der Waals interaction energy is unfavorable for the stability of the complex, as it has a positive value. Moreover, total internal energy is also positive, which disfavors the stability of the complex. We see that nonpolar interaction of the solvent molecule does not favor the stability of the complex, while low energy values of polar interactions dominate the nonpolar interaction and as a whole solvent interaction gives stability to the complex. We see that, the sum of EELE, EVAN, and EINT strongly disfavors the stability of the complex, however, the total solvation energy strongly favors the stability of the trimeric complex. Furthermore, the total energy of the complex is negative which supports the overall stability of the complex. Hence, it may be concluded that the effect of solvation plays an important role for the stability of the trimeric complex. 3.4. MD simulations at high (18 Cl − and 21 Na +) ionic concentrations We performed a separated MD simulation by increasing the number of chloride and sodium ions using the same parameters as

described for previous simulations. The initial structure of the complex with ions is shown in Fig. 4S. All chloride ions are put into the channel while sodium ions are distributed around the complex. The motion of these ions during the entire simulation is shown in Fig. 6. From the diffusion of ions, we see that several ions bind with the trimer and show least diffusion with respect to the complex. We found that a Cl − 1196 also binds with the complex and lives as a life partner of the complex. We see that several water molecules are in good contact with sodium ions outside and inside, as well. The root mean square deviations and B-factor are also depicted in Fig. 7a and b. From Fig. 7, we see that the increment in the chloride and sodium ions leads to a stable conformation of the trimeric complex which is validated by small fluctuations of all units of the trimer. In the previous simulation, we notice that the deviation for unit III was very high which was responsible for the instability of the complex. By MD movie (Supplementary material Video data VS3), we see that several sodium ions and chloride ions interact with each unit of the trimer. These interactions provide stability to the complex. We observe that the complex has very small root mean square deviation (b3 Å), which elucidates a well equilibrated and stabilized system.

58

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64 Table 3 Various energy contributions during stability calculation of the trimeric complex in the presence of few chloride ions (3 Cl−). All energy parameters are calculated in kcal/mol. All parameters have the usual meaning as described in Table 1. Contributions

Mean

Std dev

EELE EVDW EINT EGAS ENONPOL EPOL EELE+POL EVDW+NONPOL ESOL ETOT

−25,802.45 9992.70 25,731.98 9922.23 473.81 −14,339.73 −40,142.18 10,466.51 −13,865.92 −3943.68

278.42 810.26 79.20 1056.19 3.77 239.01 258.71 235.88 891.45

3.5. MD at pH ~ 4

Fig. 5. The diffusion of ions in the trimeric complex at neutral pH with 5 chloride ions. Ions closer to the complex are shown in different colors whose information is indicated in the legend box. Here, the unit of diffusion is Å2 (mean square displacement).

The calculation of the residue-wise B-factor of the entire system was also performed to enumerate the residue fluctuations during simulation. Fig. 7b elucidates that thermal fluctuations (B-factor) of residues for all units are almost similar. However, we see that some residues (Gly104 of unit II, Asp548 for unit II and Asp942 for unit III) have the largest fluctuations during simulations. The large B-factor of a residue shows large flexibility in that residue. Hence, these residues have the highest fluctuation rate during simulations. We used the standard MM–PBSA method to investigate the stability of the trimeric complex in a high concentration of ions. The results from the stability calculation are given in Table 4. From the table we see that electrostatic energy also favors the stability of the complex. The electrostatic energy has the largest contribution among all other energy parameters. Hence, for the stability of the trimer, it is the most important interaction. The positive van der Waals interaction energy does not favor the stability of the complex. We notice that the internal energy is also very high. The internal energy is a significant molecular motion; hence it arises from the atomic fluctuations during simulations. The non polar interactions do not contribute in the stability of the trimer, as it is positive. However, polar contributions have a significantly high value which strongly favors the stability of the complex. As we see the favorable contributions due to electrostatic interaction, are dominated by unfavorable contributions due to internal and van der Waals energy, hence polar contributions are the major factor which stabilizes the complex. We observe that the total energy, is negative, and it shows significantly strong contributions. The solvation energy is the major contributor to the favorable total energy. We know that polar energy arises due to the interaction of the solvent with the complex system; hence water molecule plays a crucial role for the stability of the complex. This speculation is firmly validated by analysis of the trajectories. We see that several water molecules are in interaction with the complex at various places during the entire simulation. The interactions of these water molecules give a favorable solvation energy as is clear from Table 4. We may conclude that the increment in the concentration of Na + leads to a stable conformation of the trimeric complex. Since excess of Na + ions produces an environment of more positive charges analogous to low pH, we have performed a separate simulation at low pH.

To study the conformational behavior of the trimeric complex we performed a separate MD simulation at low pH (pH ≈ 4). The initial structure of the complex at low pH in the presence of chloride ions is shown in Fig. 5S. Due to protonation of all histidines, glutamic acids and aspartic acids, there was a significant increment in the charge of the complex, therefore 174 chloride ions are added in the ensemble to neutralize the system. We have inserted some ions in the channel formed by three units of the trimer to investigate the effect of these ions in the cavity. During the simulation, we saw that the chloride ion continuously moves throughout the cavity (Video data 4VS in supplementary materials). The motions of these ions are depicted in Fig. 8. We found that Cl − 1197 binds permanently with the complex and may be treated as a life partner of the system. Other chloride ions show very large diffusions. Three different snapshots representing the motion of ions are shown in Fig. 9(a–c). Furthermore, we see that the area of the channel increases significantly during simulations and the distance between each unit of trimer also increases. The increase in separation between the units occurs due to decrement in interaction between the units. The RMS variations and residue wise fluctuations during simulation are shown in Fig. 10a and b. We see that root mean square deviations for domain I and domain II are relatively smaller at the initial stage of simulations (almost till 4 ns of production dynamics). The root mean square deviation for unit III is large relative to other units. This deviation increases abruptly during 5–8 ns of simulation. However, it decreases and almost becomes steady during the rest of the simulations. The RMS deviations reveal that all units go away from each other during simulation, which weakens the interaction between the units. Due to the increased separation, the interaction between the units drops off and destabilizes the trimeric system. A similar pattern is followed by the B-factor scheme during simulations. We notice that domain III has the largest thermal fluctuations during simulations; however, domain I and domain II exhibit small thermal fluctuations. In previous simulations the trimeric complex was more stable at neutral pH. It may occur due to the presence of Na+ ions in the complex at neutral pH. At neutral pH, we observed that a sodium ion binds with domain I and domain II that provide stability to the trimeric complex. The increment in the positive charge due to a decrease of pH, disables the binding of Na+ with domain I and domain II which destabilizes the complex. After the visualization of the movie of the production dynamics and Fig. 6 (diffusion of ions) we notice that several Cl − propagate through the channel formed by the trimeric interface. Due to the propagation of chloride ions, the electrostatic interactions vary significantly, which leads to major changes in the conformation of the trimer. Therefore, we predict that low pH is not a sufficient condition for the stability of the trimer and there should be some other parameters which are also necessary

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

59

Fig. 6. The motion of ions in the trimeric complex at neutral pH and with 18 chloride ions. Ions closer to the complex are shown in different colors whose information is indicated in the legend box. Here, the y-axis shows the mean square displacements of ions measured in Å2.

for the stability of the complex in addition to pH. Furthermore, we also performed standard MM–PBSA calculations to study the stability of the trimeric complex. The energy contributions during MM–PBSA calculations are shown in Table 5. From the table, we see that electrostatic energy (71,407.12 kcal/mol) does not show favorable interactions within the trimeric complex. Furthermore, van der Waals interaction does not favor stability of the complex as well. However, we notice that there are favorable polar and solvation interactions in the complex. From the snapshots of the complex during simulation, we see that the solvent accessible area of the trimeric complex increases during simulations. The table for the solvent accessible surface area (SASA) for all simulated systems is shown in Table 1. Due to increase in solvation area, more solvents are being interacted by the complex and hence there is good solvent interaction as inferred from the data. Therefore, we see that simulation at very low pH (≈ 4) is not meaningful regarding the stability of the trimer. Hence we have performed a separate simulation at different ionic concentrations and pH conditions using identical simulation parameters. 3.6. Molecular dynamics simulations at pH ≈ 6 We have protonated all histidines to produce a condition whose pH is in the range of 6, while aspartic acids and glutamic acids are left unaltered. The initial structure prior to equilibration is shown in Fig. 6S. We arranged all chloride ions into the cavity formed by contacting the interface of the trimer. Some chloride ions are also arranged outside the trimer as shown in Fig. 6S. Since the addition of protons on histidines increases the charge of the system, no further sodium ion is required to neutralize the system. The motion of ions into the complex is shown in Fig. 11. We found that several chloride ions bind with the trimeric complex and become a life partner of the complex. These ions have utmost importance in the stability of the trimeric complex. The root mean square deviations and the thermal fluctuations are shown in Fig. 12a and b. From the figure, we see that during the equilibration, all units of trimer are almost stable while deviation increases during production dynamics. Furthermore, the second unit exhibits a stable RMSD relative to other units while the first and third units show significantly

more RMS deviation relative to the second unit. However, for longer simulation, unit I also shows a stable RMSD value, but the third unit has a significantly larger RMS deviation during the entire simulations. If we compare these results with the very low pH simulations, we notice that in the later case (pH ≈ 6) the trimeric complex is stable relative to the previous simulation at pH ≈ 4. The thermal fluctuation is helpful in determining the flexibility of residues in the complex. We see that the thermal fluctuation of the third unit is largest while it is least for the second unit. We see that, the thermal fluctuation schemes of all three units are entirely different. For the first unit, we see that the B-factor for residues 145–150 has the highest root mean square fluctuation; hence these residues are the most flexible part of unit I. These residues belong to a loop of domain I which is composed mostly of the glutamates and histidines as already shown in Fig. 1. Furthermore, residues 300–394 are also more flexible relative to other residues of unit I, which belong to the domain III of the first unit. The result is also validated by previous studies where the flexibility of domain III is already stated [55,56]. Comparing domain-wise fluctuations, we found that domain I and domain II are relatively less flexible than domain III. Domain III has almost similar fluctuations in unit I and unit II. The fluctuations of domain I are entirely different for unit I and unit II. Domain I was almost stable in unit I while it shows a significant fluctuation in unit II. Residue 1078 of unit III shows the largest fluctuation during simulations. This residue is part of domain I which connects domain III and comprises a loop. The fluctuations of different parts of the trimer can also be explored with the help of RMS deviation. We see that, the RMS deviations of the third unit, as shown in green lines in Fig. 10a, show more deviations during simulations. Hence, the B-factor for this unit is also high. Similarly, the RMS deviation for the second unit is less than the RMS deviation of the first unit; hence thermal fluctuation is relatively small for the second unit. The above results elucidate that at pH ~ 6 the trimeric complex has stable conformations, which is validated by the B-factor as well as the root mean square deviations. We have used the standard MM–PBSA method to investigate the stability of the complex thermo-chemically. The various energy contributions as the result of stability calculation are shown in

60

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

Fig. 7. (a) Root mean square deviations and (b) residue wise thermal fluctuations, during the entire simulation in the presence of 18 chloride ions at neutral pH. Colors have the same meaning as described earlier.

Table 6. We see that electrostatic contribution (−22,921.31 kcal/mol) significantly favors the stability of complex. The van der Waals energy is again positive (10,172.83 kcal/mol) which does not favor the stability of the complex. The total internal energy of the complex is also positive (22,921.66 kcal/mol) and it dominates favorable electrostatic interactions. The sum of van der Waals interactions and internal energy makes the total bonded interactions (EGAS) positive. The positive non Table 4 Various energy terms in the stability of the trimeric complex in the presence of more chloride ions (18 Cl−) at neutral pH. All energy parameters are in kcal/mol and follow the same explanations as in Table 1. Contributions

Mean

Std dev

EELE EVDW EINT EGAS ENONPOL EPOL EELE+POL EVDW+NONPOL ESOL ETOT

−25,344.84 9939.18 25,640.32 10,234.67 477.40 −14,797.11 −40,141.95 10,416.58 −14,319.71 −4085.05

189.14 536.11 108.79 663.33 2.77 173.82 181.48 171.90 573.93

polar energy does not favor the stability of the complex. The non polar interaction arises due to the cavitations of the solvent molecule into the complex. However, the polar contributions are significantly large which strongly favors the stability of the complex. Since the polar interaction energy is very large relative to the non polar interactions, the total solvation energy favors the stability of the trimeric complex. In addition the total energy (sum of all energy terms calculated by PB methods) is also negative, which is a strong evidence of the stability of the complex. 4. Discussion Three simulations at neutral pH reveal that the trimer is in a stable form at neutral pH and high concentration of ions. Therefore, some sodium ions are necessary for the stability of complex and they interact with the complex for a longer simulation time. Since sodium ions do not interact with the third unit of the trimer, its deviations are in some extent larger than the first and second units. In Section 3.3 where the numbers of chloride and sodium ions are significantly small, the deviation from the initial structure is very high. In this case, unit III has an abnormal deviation at larger simulations and it dissociates from other units which destabilize the complex. However,

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

61

Fig. 8. The motion of ions in the trimeric complex at very low pH (pH ~ 4). Protonation of His, Glu and Asp produced extra positive charges, therefore, only chloride ions were added to neutralize the complex. Ions closer to the complex are shown in different colors whose information is indicated in the legend box. The motion of ions on the y-axis is represented by mean square displacement measured in Å2.

when we increase the concentration of sodium ions, the deviation becomes significantly small relative to the previous one. From the thermochemical data, we see that the increment in the concentration of chloride ions, increases the electrostatic interactions but it also lessens the van der Waals interactions. We notice that, electrostatic interaction of the trimeric complex is least at neutral pH and highest in the presence of 5 Cl − ions. However, the RMS deviation elucidates that the trimeric complex will have the least structural deviations at neutral pH without Cl− ions. Hence, we may conclude that the electrostatic interaction does not affect significantly due to the variation of salt concentrations. We notice that the van der Waals interaction for neutral pH, which is the most stable system according to the RMS deviations, has the least positive energy while, it is the highest for simulations with 5 chloride ions. Therefore, van der Waals interaction has reasonable effect on the stability of the trimeric complex. By low pH simulations we see that, change in pH strongly affects the electrostatic interaction. Electrostatic energy at pH = 4, has a high positive value, which disfavors the stability of the complex. The protonation of all histidines, aspartate and glutamate gives a large amount of positive charges on each unit. The repulsion between similar charge groups dominates all favorable electrostatic interactions and as a result electrostatic interaction becomes highly positive. Furthermore, we notice that no other interactions are favorable except the solvation energy. We notice a large conformational rearrangement in each unit of trimer which gives the entire system a staggered shape. This conformational rearrangement increases the solvent exposed surface area of the entire system which is directly responsible for the solvation interaction energy. We notice that several water molecules interacted with the trimeric complex, which are also responsible for favorable solvation energy. From thermodynamical respect, we see that the trimer will not be stable at very low pH (pH = 4). The protonation of histidines only, produces a pH of around 6–6.5. Therefore, the complex at this pH is significantly stable as is obvious in the root mean square and thermal fluctuation scheme. We see that the total energy (GTOT) is significantly negative and it strongly favors the stability of the complex. Comparing results from the binding energy we see that the increase of concentration of sodium ions continuously increases the stability of the complex. At neutral pH,

Fig. 9. (a–c) Snapshots during simulations at low pH (~4) at different time steps (4 ns, 8 ns and 12 ns respectively). The separation between units is shown by an arrow.

when concentrations of sodium ions were least, the total energy of the complex is − 3086.58 kcal/mol, which is lowest among all simulations. When we increase the concentration of ions at neutral pH (for Section 3.3), we get the total energy of the system as −3943.68 kcal/mol. Furthermore, at a higher concentration of sodium ions at neutral pH (Section 3.4), the total energy has been slightly increased. These results show a transition towards more stable conformation due to the increase of ionic concentrations. Total energy is significantly better for pH = 6 than neutral pH. 5. Conclusion From the above results, we see that the role of salts is crucial for the stability of the trimeric complex. The increment in the

62

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

Fig. 10. (a) The root mean square deviation and (b) thermal fluctuation (B-factor) residue-wise during simulation at low pH (~4).

concentration of sodium ions increases the stability of the complex in the neutral pH, which indicates a nucleophillic character of the receptor. At low pH, there may be excess protons which may interact with the nucleophillic receptor to stabilize the complex relative to

Table 5 Different energy contributions in the stability of the trimeric complex at very low pH (pH ~ 4). All energy values are in kcal/mol and follows earlier explanations. Contributions

Mean

Std dev

EELE EVDW EINT EGAS ENONPOL EPOL EELE+POL EVDW+NONPOL ESOL ETOT

71,407.12 10,416.58 26,612.62 108,436.31 492.37 −106,617.15 −35,210.95 10,908.85 −106,124.78 2311.54

199.78 879.42 135.61 1056.99 2.23 195.66 197.72 196.26 973.90

Fig. 11. The motion of ions in the trimeric complex at pH~6. Ions closer to the complex are shown in different colors whose information is indicated in the legend box. The motion of ions on the y-axis is represented by mean square displacement and is measured in Å2.

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

63

Fig. 12. (a) The root mean square deviation and (b) residue wise thermal fluctuation (B-factor) during simulation at pH ~ 6. Colors have the same meaning as described earlier.

neutral pH. There are few chloride ions which bind the trimer at low pH and may be responsible for the increased electrostatic interaction. Therefore, we conclude that pH alone is not a factor but the role of ions is also very crucial for the stability and mechanism of the dengue envelope protein. In our recent study we have presented the stability of a dimer at low and high pH separately, and found that a decrease in pH drastically changes the interaction between dimeric partners [33]. In the present study, we see that the stability of the trimeric complex at pH six is better than the stability of the complex at neutral pH. We have also investigated that if we decrease pH to below six, the trimer becomes unstable and there is no favorable interaction between the trimeric partners. Therefore, we speculate that the trimer will be in

the most stable form at pH 6. We also speculate that the protonation of histidine is crucial for the stability of the complex. We found that electrostatic interaction is significantly affected by change in pH, however, van der Waals interaction plays a crucial role in binding of the trimeric partner. Supplementary data to this article can be found online at http:// dx.doi.org/10.1016/j.bbapap.2012.08.014. Acknowledgements The authors acknowledge the Department of Science and Technology for the FIST scheme. Major parts of the computational works were

64

K.D. Dubey et al. / Biochimica et Biophysica Acta 1834 (2013) 53–64

Table 6 Energy contributions in the stability of the trimeric complex at pH ~ 6. All energy parameters are measured in kcal/mol. Contributions

Mean

Std dev

EELE EVDW EINT EGAS ENONPOL EPOL EELE+POL EVDW+NONPOL ESOL ETOT

−22,921.31 10,172.83 25,966.25 13,217.77 479.53 −17,665.96 −40,587.27 10,652.36 −17,186.43 −3968.66

54.60 777.31 80.19 829.29 2.12 65.71 60.30 65.11 856.16

performed on the SCFBIO, Indian Institute of Technology, New Delhi (IITD). Authors are thankful to Prof. B. Jayram, IITD, for facilitating supercomputer access at IITD. KDD acknowledges the Council for Scientific and Industrial Research for the SRF. References [1] G.N. Malavige, S. Fernando, D.J. Fernando, S.L. Seneviratne, Dengue viral infections, Postgrad. Med. J. 80 (2004) 588–601. [2] B.D. Lindenbach, H.-J. Thiel, C.M. Rice, Flaviviridae: the viruses and their replication, in: D.M. Knipe, P.M. Howley (Eds.), Fields Virology, 5th ed., Lippincott Williams and Wilkins, Philadelphia, PA, 2007, pp. 1102–1152. [3] R.J. Kuhn, W. Zhang, M.G. Rossmann, S.V. Pletnev, J. Corver, E. Lenches, C.T. Jones, S. Mukhopadhyay, P.R. Chipman, E.G. Strauss, T.S. Baker, J.H. Strauss, Structure of dengue virus: implications for flavivirus organization, maturation, and fusion, Cell 108 (2002) 717–725. [4] Y. Modis, S. Ogata, D. Clements, S.C. Harrison, A ligand binding pocket in the dengue envelope protein, Proc. Natl. Acad. Sci. 100 (2003) 6986–6991. [5] Y. Modis, S. Ogata, D. Clements, S.C. Harrison, Structure of the dengue virus envelope protein after membrane fusion, Nature 427 (2004) 313–319. [6] V. Nayak, M. Dessau, K. Kucera, K. Anthony, M. Ledizet, Y. Modis, Crystal structure of dengue virus type 1 envelope protein in the postfusion conformation and its implications for membrane fusion, J. Virol. 83 (2009) 4338–4344. [7] S.C. Harrison, Viral membrane fusion, Nat. Struct. Mol. Biol. 15 (2008) 690–698. [8] S. Bressanelli, K. Stiasny, S.L. Allison, E.A. Stura, S. Duquerroy, J. Lescar, F.X. Heinz, F.A. Rey, Structure of a flavivirus envelope glycoprotein in its low-pH-induced membrane fusion conformation, EMBO J. 23 (2004) 728–738. [9] M. Kielian, Class II virus membrane fusion proteins, Virology 344 (2005) 38–47. [10] F.A. Rey, Dengue virus envelope glycoprotein structure: new insight into its interactions during viral entry, Proc. Natl. Acad. Sci. 100 (2003) 6899–6901. [11] Y. Modis, S. Ogata, D. Clements, S.C. Harrison, Variable surface epitopes in the crystal structure of dengue virus type 3 envelope glycoprotein, J. Virol. 79 (2005) 1223–1231. [12] M. Kielian, F.A. Rey, Virus membrane-fusion proteins: more than one way to make a hairpin, Nat. Rev. Microbiol. 4 (2006) 67–76. [13] G.E. Nybakken, T. Oliphant, S. Johnson, S. Burke, M.S. Diamond, D.H. Fremont, Structural basis of West Nile virus neutralization by a therapeutic antibody, Nature 437 (2005) 764–769. [14] J.W. Lee, J.J. Chu, M.L. Ng, Quantifying the specific binding between West Nile virus envelope domain III protein and the cellular receptor alphaVbeta3 integrin, J. Biol. Chem. 281 (2006) 1352–1360. [15] J.P. Dal Molin, M.A.A. da Silva, A. Caliri, Effect of local thermal fluctuations on folding kinetics: a study from the perspective of non extensive statistical mechanics, Phys. Rev. E 84 (2011). [16] V. Frecer, S. Miertus, Design, structure-based focusing and in silico screening of combinatorial library of peptidomimetic inhibitors of dengue virus ns2b-ns3 protease, J. Comput. Aided Mol. Des. 24 (2010) 195–212. [17] A.R. Leach, M.M. Hann, Molecular complexity and fragment-based drug discovery: ten years on, Curr. Opin. Chem. Biol. 15 (2011) 489–496. [18] X. Chen, A. Tropsha, Calculation of relative binding affinity of enzyme inhibitors using the generalized linear response methods, J. Chem. Theory Comput. 2 (2006) 1435–1443. [19] H. Gohlke, D.A. Case, Converging free energy estimates: MM–PB(GB)SA studies on the Protein–protein complex Ras–Raf, J. Comput. Chem. 25 (2003) 238–250. [20] J. Aqvist, C. Medina, J.E. Samuelson, A new method for prediction binding affinity in computer aided drug design, Protein Eng. 7 (1994) 385–391. [21] U. Bren, V. Martinek, J. Florian, Decomposition of the salvation free energies of deoxyribonucleoside triphosphate using the free energy perturbation method, J. Phys. Chem. B 110 (2006) 12782–12788. [22] M. Bren, J. Florian, J. Mavri, U. Bren, Do all pieces make a whole? Thiele cumulants and the free energy decompositions, Theor. Chem. Acc. 117 (2007) 535–540. [23] A. Predih, U. Bren, T. Solemajer, Binding free energy calculations of N-sulphonyl‐ glutamic acid inhibitors of MurD ligase, J. Mol. Model. 15 (2009) 983–996. [24] U. Bren, J. Lah, M. Bren, V. Martinek, J. Florian, DNA duplex stability: the role of preorganized electrostatics, J. Phys. Chem. B 114 (2010) 2876–2885.

[25] U. Bren, V. Martinek, J. Florian, Free energy simulations of uncatalyzed DNA replication fidelity: structure and stability of T-G and dTTP-G terminal DNA mismatches flanked by a single dangling nucleotide, J. Chem. Phys. B 110 (2006) 10557–10566. [26] V. Martinek, U. Bren, M.F. Goodman, A. Warshel, J. Florian, DNA polymerase [beta] catalytic efficiency mirrors the Asn279-dCTP H-bonding strength, FEBS Lett. 581 (2007) 775–780. [27] T. Knehans, A. Schuller, D.N. Doan, K. Nacro, J. Hill, P. Guntert, M.S. Madhusudhan, T. Weil, S.G. Vasudevan, Structure-guided fragment-based in silico drug design of dengue protease inhibitors, J. Comput. Aided Mol. Des. 25 (2011) 263–274. [28] C. Hetenyi, D. van der Spoel, Toward prediction of functional protein pockets using blind docking and pocket search algorithms, Protein Sci. 20 (2011) 880–893. [29] R. Yennamalli, N. Subbarao, T. Kampmann, R.P. McGeary, P.R. Young, B. Kobe, Identification of novel target sites and an inhibitor of the dengue virus E protein, J. Comput. Aided Mol. Des. 23 (2009) 333–341. [30] Z. Zuo, O.W. Liew, G. Chen, P.C.J. Chong, S.H. Lee, K. Chen, H. Jiang, C.M. Puah, W. Zhu, Mechanism of ns2b-mediated activation of ns3pro in dengue virus: molecular dynamics simulations and bioassays, J. Virol. 83 (2009) 1060–1070. [31] K. Wichapong, S. Pianwanit, W. Sippl, S. Kokpol, Homology modeling and molecular dynamics simulations of dengue virus ns2b/ns3 protease: insight into molecular interaction, J. Mol. Recognit. 23 (2010) 283–300. [32] M.K. Prakash, A. Barducci, M. Parrinello, Probing the mechanism of pH-induced large‐scale conformational changes in dengue virus envelope protein using atomistic simulations, Biophys. J. 99 (2010) 588–594. [33] K.D. Dubey, A.K. Chaubey, R.P. Ojha, Role of pH on dimeric interactions of DENV envelope protein: an insight from molecular dynamic study, Biochim. Biophys. Acta 1814 (2011) 1796–1801. [34] D.A. Case, T.E. Cheatham, T. Daren, H. Gohlke, R. Luo, K.M.A. Onufriev, C. Simmerling, B. Wang, R.J. Woods, The Amber biomolecular simulation programs, J. Comput. Chem. 26 (2005) 1668–1688. [35] I.S. Joung, T.E. Cheatam III, Molecular dynamics simulations of the dynamic and energetic properties of alkali and halide ions using water-model-specific ion parameters, J. Phys. Chem. 113 (2009) 13279–13290. [36] P. Banas, D. Hollas, M. Zgarbova, P. Jurecka, M. Orozco, T.E. Cheatam III, J. Sponer, M. Otyepka, Performance of molecular mechanics force field for nucleic acids: stability of UUCG and GNRA hairpins, J. Chem. Theory Comput. 6 (2010) 3836–3849. [37] K. Lindorff-Larsen, S. Piana, K. Palmo, P. Maragakis, J. Klepesis, D.E. Shaw, Improved side-chain torsion potentials for the Amber ff99SB protein force field, Proteins 78 (2010) 1950–1958. [38] D.S. Mueller, T. Kampamann, R. Yennamalli, P.R. Young, B. Kobe, A.E. Mark, Histidine protonation and the activation of viral fusion proteins, Biochem. Soc. Trans. 36 (2008) 43–45. [39] T. Kampmann, D.S. Mueller, A.E. Mark, P.R. Young, B. Kobe, The role of histidine residues in low-pH-mediated viral membrane fusion, Structure 14 (2006) 1481–1487. [40] Z. Qin, Y. Zheng, M. Keilian, Role of conserved residues in the low pH dependence of the Semiliski Forest virus fusion, Proteins 83 (2009) 4670–4677. [41] F.A. Carneiro, F. Stauffer, C.S. Lima, M.A. Juliano, L. Juliano, A.T. Da Poian, Membrane fusion induced by vesicular stomatitis virus depends on histidine protonation, J. Biol. Chem. 278 (2003) 13789–13794. [42] W.L. Jorgenson, J. Chandrashekhar, J.D. Madura, R.W. Imprey, M. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys. 79 (1983) 926–935. [43] J.A. Izaguirre, D.P. Catarello, J.M. Wozanaik, R.D. Skeel, Langevin stabilization of molecular dynamics, J. Chem. Phys. 114 (2001) 2090–2098. [44] H.J.C. Berendsen, J.P.M. Postama, W.F. van-Gunsteren, A. DiNola, J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (1984) 3684–3690. [45] J.P. Ryckaert, G. Cicotti, H.J.C. Barendsen, Numerical integration of the Cartesian equation of motion of a system with constraints: molecular dynamics of n-alkanes, J. Comput. Phys. 23 (1977) 327–341. [46] T. Darden, D. York, L. Pedersen, Particle mesh Ewald: an N⋅log(N) method for Ewald sums in large systems, J. Chem. Phys. 98 (1993) 10089–10092. [47] W. Humphrey, A. Dalke, K. Schulten, VMD — virtual molecular dynamics, J. Mol. Graph. Model. 14 (1996) 33–38. [48] E.F. Pettersen, T.D. Goddard, C.C. Haung, G.S. Couch, D.M. Greenblatt, E.C. Meng, T.E. Ferrin, UCSF chimera — a visualization system for exploratory research and analysis, J. Comput. Chem. 25 (2004) 1605–1612. [49] Maestro 9.2, Schrodinger LLC, New York, USA, 2011. [50] F. Fogolari, A. Brigo, H. Molinari, Protocols for MM/PBSA molecular dynamics simulations of proteins, Biophys. J. 85 (2003) 159–166. [51] K.D. Dubey, R.P. Ojha, Binding free energy calculation with QM/MM hybrid method for Abl kinase inhibitor, J. Biol. Phys. 37 (2011) 69–78. [52] K.D. Dubey, R.P. Ojha, Conformational flexibility, binding energy, role of salt bridge and alanine-mutagenesis for c-Abl kinase complex, J. Mol. Model. (2011), http://dx.doi.org/10.1007/s00894-011-1199-9. [53] A.K. Chaubey, K.D. Dubey, R.P. Ojha, Stability and free energy calculations of LNA modified human Tel-24 quadruplex structure: a molecular dynamics study, J. Comput. Aided Mol. Des. (2011), http://dx.doi.org/10.1007/s10822-012-9548-z. [54] K.D. Dubey, R.P. Ojha, Conformational flexibility and binding energy profile of c‐Abl tyrosine kinase complexed with Imatinib: an insight from MD study, Mol. Simulat. 37 (2011) 1151–1163. [55] J.T. Roehrig, R.A. Bolin, R.G. Kelly, Monoclonal antibody mapping of the envelope glycoprotein of the dengue 2 virus, Jamaica, Virology 246 (1998) 317–328. [56] D.L. Gibbsons, M.C. Vaney, A. Roussel, A. Vigouroux, B. Reily, J. Lepault, M. Keilian, F.A. Rey, Conformational change and protein-protein interactions of the fusion of semilski forest virus, Nature 427 (2004) 320–325.