Molecular Dynamics Simulation of Barnase: Contribution of Noncovalent Intramolecular Interaction to Thermostability

Chen, Zhiguo; Fu, Yi; Xu, Wenbo; Li, Ming

doi:https://doi.org/10.1155/2013/504183

Mathematical Problems in Engineering

On this page

Abstract Introduction Results and Discussion Conclusion Acknowledgments References Copyright Related Articles

Special Issue

Stochastic Systems: Modeling, Optimization, and Applications

View this Special Issue

Research Article | Open Access

Volume 2013 | Article ID 504183 | https://doi.org/10.1155/2013/504183

Molecular Dynamics Simulation of Barnase: Contribution of Noncovalent Intramolecular Interaction to Thermostability

Zhiguo Chen,¹Yi Fu,¹Wenbo Xu,¹and Ming Li²

Academic Editor: Shuping He

Received10 Oct 2013

Accepted07 Nov 2013

Published17 Dec 2013

Abstract

Bacillus amyloliquefaciens ribonuclease Barnase (RNase Ba) is a 12 kD (kilodalton) small extracellular ribonuclease. It has broad application prospects in agriculture, clinical medicine, pharmaceutical, and so forth. In this work, the thermal stability of Barnase has been studied using molecular dynamics simulation at different temperatures. The present study focuses on the contribution of noncovalent intramolecular interaction to protein stability and how they affect the thermal stability of the enzyme. Profiles of root mean square deviation and root mean square fluctuation identify thermostable and thermosensitive regions of Barnase. Analyses of trajectories in terms of secondary structure content, intramolecular hydrogen bonds and salt bridge interactions indicate distinct differences in different temperature simulations. In the simulations, Four three-member salt bridge networks (Asp8-Arg110-Asp12, Arg83-Asp75-Arg87, Lys66-Asp93-Arg69, and Asp54-Lys27-Glu73) have been identified as critical salt bridges for thermostability which are maintained stably at higher temperature enhancing stability of three hydrophobic cores. The study may help enlighten our knowledge of protein structural properties, noncovalent interactions which can stabilize secondary peptide structures or promote folding, and also help understand their actions better. Such an understanding is required for designing efficient enzymes with characteristics for particular applications at desired working temperatures.

1. Introduction

The stability of thermophilic proteins has been viewed from different perspectives and there is yet no unified principle to understand this stability. Detailed knowledge of the thermal stability of proteins and the phenomenon of protein folding are essential not only to understand protein structure and function but also to design thermostable proteins for industrial applications [1–5]. Comparisons of homologous proteins from thermophiles and mesophiles have shown that related proteins can perform similar functions yet have very different stabilities [6–8]. Recent comparisons between the structures of mesophilic and hyperthermophilic proteins have identified a number of structural features that are believed to give rise to increased thermal stability, such as optimized electrostatic interactions [9–11], integrating disulfide bridges [12], compacting of the hydrophobic core [13], the greater number of salt bridges [14–17], and stabilization of the α-helices [18]. It is clear that many different factors and methodologies [19–21] can contribute to enzymatic temperature adaptation and that no single factor can be invoked to explain adaptation in general. Although considerable progress has been made in theoretical and experimental studies of the protein folding and thermal stability, our knowledge is still limited for fully understanding this subject especially on mathematical modeling [22, 23]. Therefore, further investigations of the protein folding and thermal stability mechanisms are crucial since it may provide relevant information on the evolutionary aspects involved and on the general mechanisms underlying protein stability. Moreover, it may help in the design of thermostable biomolecules that are functional at high temperature.

Barnase (E.C.3.1.2.7) is an extracellular 110 residue ribonuclease from Bacillus amyloliquefaciens. It has broad application prospects in agriculture, clinical medicine, pharmaceutical, and so forth, especially that its special biological function of resistance to virus and tumor attracts more researchers’ attention recently [24–27]. Barnase is a multidomain protein with three helices in the first half of the sequence, followed by a five-stranded antiparallel β-sheet. There are three hydrophobic cores in protein. It is well suited for studying protein folding and stability because it is one of the smallest globular proteins that does not contain any prosthetic groups, metal cofactors, or disulfide bonds which all contribute to protein stability. It has been well studied [28–31] whereas the role of electrostatic interactions in Barnase thermal stability has not yet been investigated in detail. The question of whether salt bridges stabilize the native state of proteins has yet not definite answer. This is most likely related to the fact that the location, geometry, and optimization of the electrostatics vary greatly in proteins. We aim to address this issue by MD (molecular dynamics) simulation, and the major focus of this study was to investigate the contributions of noncovalent intramolecular interactions to protein stability, and how changes in these interactions are correlated with the initial events associated with tertiary structure unfolding in response to thermal stress.

The small size of the protein, the absence of metal cofactors, and the absence of disulfide bonds make Barnase an ideal test case for studying the contribution of salt bridge to thermostability by computer simulation. In this paper we have performed molecular dynamics simulation of Barnase at five different temperatures, namely, 300 K, 400 K, 500 K, 550 K, and 600 K. By calculating the root mean square deviation (RMSD) and root mean square fluctuation (RMSF) values for backbone and Cα atoms, the thermal sensitive regions have been identified. The dynamic properties of Barnase at different temperatures have been compared in terms of secondary structure content, molecular flexibility, intramolecular hydrogen bonds, and salt bridge interactions. Our analysis specifically focuses on the contribution of salt bridge to the stability of the protein along the unfolding transition. The study provides insight into the structure-stability relationship of Barnase, which may help enlighten our knowledge of protein structural properties and noncovalent interactions that stabilize secondary peptide structures or promote folding. A better understanding of specific interactions change during the unfolding simulation also allows us to predict more precisely a specific site or region critical to increase Barnase’s thermostability.

2. Methodology

2.1. Molecular Dynamics Simulation

Molecular dynamics simulations [32] are based on classical mechanics such as Newton’s second equation of motion. If the force acting on each atom is known, the acceleration in a system can be acquired. The equations of motion can be determined via the acceleration, resulting in a trajectory of the positions, velocities, and accelerations of each particle in given system. The force on an atom with mass and position is determined from the potential of the system , as shown in (1). The potential energy can be used to derive the position of a particle as a function of time in (2). The equations of motion are got for each atom in the system. The forces acting on atoms in new positions can be calculated and the simulation will continue as many time steps as necessary:

2.2. The Potential Function

Empirical energy functions can only be used to approach the force fields for large biological systems composed of many atoms. Here the CHARMM (Chemistry at HARvard Macromolecular Mechanics) force field is selected, which has been developed and continuously revised to better match new experimental data for over 25 years.

The potential function is a sum of some interaction energies. The value of the potential is determined by summing the bonded terms with the nonbonded potential terms : The nonbonded energy terms specify interactions between nonbonded atoms and atoms that are more than three covalent bonds away from each other in a molecule. These energy terms are modeled by the van der Waals energy and electrostatic energy, namely, the first and second summation items in (4) separately: where , and represent the Lennard-Jones diameters of the and atoms, respectively, , with and being the Lennard-Jones well-depth of the and atoms, and is the effective charge on atom .

The nonbonded terms calculation in the potential function is the most time-consuming part of the MD simulation. Generally, the interactions between every pair of atoms should be calculated definitely, meaning that, for an N-atom system, number of operations would be required. To reduce the computation complexity, methods were developed to ignore the interactions between two atoms separated by a distance greater than given cutoff distance. For the van der Waals interactions, the potential is “shifted” off over a distance from to . Atoms farther from each other than the distance are supposed not to interact. For the electrostatic interactions, the Ewald summation [33] is used to separate the potential into a slowly decaying long-range component and a quickly varying short-range component.

2.3. Temperature and Pressure Control

The Langevin dynamics [34] is used for constant temperature control. Langevin dynamics is an effective method for controlling the kinetic energy of the system, thus controlling the temperature and/or pressure. This methodology uses the Langevin equation [35] for a single particle. The friction kernel in Langevin equation is taken to be space and time independent for each particle, and the influence of the environment on the internal force is denoted as an average sense. Thus, explicit hydrodynamic interactions can be ignored, and the internal force is introduced by a frictional term proportional to the velocity and a random force , which approximately simulate molecular collisions and viscosity in the realistic cellular circumstance: Here is the mass matrix, is the potential energy governing the solute, is the damping coefficient, and is a Gaussian white noise force vector that has mean zero.

From classic theories on Brownian motion, it can be seen that, although molecular collisions are random, the ensemble of these collisions produces a systematic effect. That is to say, the molecular random motions exist at thermal equilibrium as a fluctuation. Hence one can see that the frictional force is a correlate of the random force by the fluctuation/dissipation theorem [35]. This relation can be expressed by the -dependence of the covariance of : where is the Boltzmann constant, is the target temperature, and is the usual Dirac symbol.

The Dirac- function is in the form

The random force is chosen independently for each step. And the covariance matrix is diagonal as hydrodynamic interactions between particles have been discounted.

A physical value of for each particle can be selected according to Stokes’ law for a hydrodynamic particle with radius . Stokes’ law describes how the frictional resistance of a spherical particle in solution varies linearly with its radius. The practical force magnitude is times the particle’s velocity, where is the solvent viscosity. Stokes’ law is often applied to particles of molecular size. Therefore, the in Langevin equation can be expressed as where is the particle’s mass.

The damping coefficient controls not only the magnitude of the frictional force but also the variance of the random forces. It can ensure that the system converges to a Boltzmann distribution characterized by the temperature . The larger the value of , the greater the influence of the surrounding fluctuating force. Small value of implies inertial motion. In this study, the main objective is to control the temperature, so we need to use small value of . The temperature of the system is maintained via the relationship between and .

For constant pressure and temperature simulations in which Langevin dynamics are used to control temperature, the pressure can be controlled in NAMD with a modified Nose-Hoover method. This method entitled Nose-Hoover [36–38] adds a fictive degree of freedom to the physical system with “coordinate” parameter (effectively a scaling parameter [39]), mass , and thermodynamic friction coefficient . (This friction coefficient is relative to and the corresponding momentum .) Besides the effective coordinate, mass, and friction set associated with the fictive thermostat variable, a set associated with virtual pressure piston (barostat) is also adopted. The effective equations of motion for a 3-dimensional system can be expressed as where is an external volume variable, is the number of degrees of freedom in the system, is the external applied pressure, and is the internal pressure, defined as

The internal virial vir is proportional to the inner product of the each atom’s position vector with the corresponding force component acting on atom due to all particles :

The conserved quantity under these augmented equations of motion is

By this means, the magnitude of the system fluctuates under specified thermostat and barostats, and the system is driven to steady state at which the average internal pressure is equal to the external applied force .

3. Results and Discussion

All simulations were performed on a PC with a Pentium 4 2.8 GHz dual core processor running Windows operating system and using the molecular dynamics program NAMD [40] with CHARMM27 [41] force fields. In order to run MD simulation, we need to do the following things.(1)A Protein Data Bank (pdb) file which stores atomic coordinates and/or velocities for the system is needed. The coordinates for starting configurations Barnase was obtained from the Protein Data Bank (PDB entry codes 1RNB [42]), which consisted of 110 residues.(2)A Protein Structure File (psf) which contains all of the molecule-specific information needed to apply a particular force field to a molecular system is needed. Coordinates of the atoms that were missing in the crystallographic structure were reconstructed using the PSFGEN structure building utility, a module of NAMD.(3)The protein needs to be solvated and put inside water, to more closely resemble the cellular environment. The protein was solvated in a cubic box consisting of TIP3 water molecules [43] with periodic boundary conditions. The system was neutralized by adding ions (Cl−) at physiological concentration using VMD’s solvate and autoionize plugins [44].(4)We need a force field parameter file. A force field is a mathematical expression of the potential which atoms experience in the system. A CHARMM forcefield parameter file contains all of the numerical constants needed to evaluate forces and energies, given a PSF structure file and atomic coordinates. The CHARMM parameters are available for download from the website: http://www.charmm.org/.(5)Create the simulation script, in which we specified all the options that NAMD should adopt in running a simulation. NAMD parses its configuration file using the Tcl scripting language.(i)First we specify the files that contain the molecular structure and initial conditions. Setting the Tcl variable temperature makes it easy to change the target temperature for many options. The outputName prefix will be used to create all of the trajectory, output, and restart files generated by NAMD run.(ii)Next is the parameter file itself and the options that control the nonbonded potential functions. These are mostly specified by the CHARMM force field. In force-field parameters, the cutoff distance was specified to 12 Å. Electrostatic interactions were calculated using the Particle Mesh Ewald (PME) summation scheme. Turn on switching for the van der Waals interactions, which were calculated with a switching function from 10 Å to 12 Å. SHAKE method [45] was used for constraining the bonds with hydrogen. The number of time steps between each output was 2 fs. Set up the temperature and pressure controllers. The Langevin was turned on and the value of the Langevin damping coefficient was set to 5/ps. The value of the Langevin temperature was set equal to the target temperature for the simulation of temperature control. The pressure control of the system was set to 1 atm.(6)The system was subjected to energy minimization for 1000 steps by steepest descents and subsequently equilibrated for 500 ps, and then the equilibrated system was subjected to molecular dynamics simulations for 2 ns each at five different temperatures, namely, 300 K, 400 K, 500 K, 550 K, and 600 K. The coordinates were saved at every 500 time steps.

3.1. Global Structural Stability

MD simulations generate an ensemble of conformations and thus include valuable information of the protein dynamics. In the following we present a detailed analysis of the four molecular dynamics trajectories, in water, generated for the protein of Barnase. The RMSD of the backbone atoms of the protein from the starting structure over the course of simulation may be used as a measure of the conformational stability of a protein during the simulation. The plots of RMSD of the protein versus time at different temperatures are shown in Figure 1(a). The plots show that the MD simulation of enzyme at 300 K is very stable throughout the simulation time. In the trajectory run at 400 K the backbone RMSD increases slightly from the starting conformation, which fluctuates, between 1.1 Å and 1.7 Å during the simulation. The average value of RMSD is about 1.20 Å in 400 K simulation, slightly above the value of 0.96 Å for the 300 K structure simulation. The curve corresponding to the 500 K simulation fluctuates more and displays a sharp rise (1.73 Å) about 662 ps, after that RMSD increases further and oscillates between 2.36 Å and 3.71 Å for the remainder of the simulation. At 600 K, the RMSD rises to about 3 Å over the first 200 ps of the simulation and reaches a value of 3.78 Å at around 900 ps. Therefore, it records a descent to about 2.56 Å at 1184 ps and another sharp rise to 4.16 Å after 1295 ps. The rise in RMSD indicates large changes to the protein structure and some disruption of the tertiary structure of the protein.

(a)

(b)

3.2. Structural Flexibility

Following molecular dynamics simulations at multiple temperatures, it was of interest to determine which general regions of the Barnase polypeptide chain exhibit hypersensitivity to thermal stress. The results are shown in Figure 1(b). There is a relatively small bump-like peak for the structure at 300 K. At 400 K, the loops between β-sheet1 and β-sheet2; β-sheet2 and β-sheet3 show significant increase together with the N- and C-terminals. The interesting finding is a dramatic increase in RMSF values between residues 5 and 7 at 500 K, which implies the beginning of the unfolding process for the N-terminal loop region and the N-terminal of α-helix1 region. As the temperature increases, these above-mentioned peaks generally become more pronounced. This pattern is especially noted in the 600 K simulation. As can be seen, most of the changes occurred in the loop region, N-terminal, and C-terminal. The regular secondary structure regions such as α-helix and β-sheet showed much less mobility during the simulations. The curves observed for four temperature simulations exhibited more or less similarly distributed fluctuations. Only at high temperature simulations most of the residues become highly mobile; therefore the curve shows a lot more fluctuation. This is due to the loss of secondary structure at these high temperatures. These observations clearly reveal the different behavior of the residues of Barnase molecule in response to increasing thermal stress and give an indication of the regions of the Barnase polypeptide chain which are most sensitive or responsive to heating.

3.3. Secondary Structure Analysis

At high temperature simulations, dramatically increased RMSD values were observed for the loop regions and both terminal regions in Barnase, which indicate these regions extensive local conformational changes upon thermal unfolding. A close analysis of the time evolution of the secondary structure (Figure 2) can present further information about its structural flexibilities.

(a)

(b)

(c)

(d)

Figure 2(a) reveals that α-helices, β-sheets, and loops observed within Barnase structure are maintained stably throughout the whole simulation period at 300 K. The overall conformation, hydrophobic core compactness, and secondary structural elements are all stable, and there is no water penetration into the protein. The simulation at 300 K also appears to agree well with the stability of the RMSD and RMSF curves (Figure 1). In case of 400 K simulation, the structure exhibits slight deviation from starting structure. There is a high degree of similarity between the graphs corresponding to the simulations at 300 K and 400 K. Only marginal structural fluctuations were observed and no significant structural changes.

At 500 K, protein structure fluctuation is significantly more pronounced. In this time period the dominant structural change was the expansion of the protein in response to the temperature increase, and the packing density in the three hydrophobic cores decreased. During the simulation, the edge residues, particularly those at or near the N-terminal part of α-helix1, are less stable and unfold first. At 2 ns simulation, α-helix1 maintains a regular shape; α-helix2 was perturbed, but mostly at the termini; α-helix 3 was gradually unfolded, after 1687 ps; α3 was lost completely. Among the five β-sheets, the first obvious observation is that β-sheet1 of the structure is partially unfolded after 959 ps, and a shortening of β-sheet1 is observed for residues 50, 51, and 56. A similar shortening also occurs at β2, β3, and β4 for residues 71, 75, 91, and 96. As shown in Figure 2(c), it is found that the β-sheet began to unfold at the edges and associated turns, and the center of the sheet is mostly stronger than the edges. In addition, the loops and turns unfold to various degrees. However, in spite of these important fluctuations in the protein, it appears that the main chain still shows essentially the same overall fold as in the native structure, and most native secondary structure elements remain present until the end of simulation (Figure 3).

When the temperature is increased to 550 K, the structure of the protein shows a continuous and progressive unfolding. The N terminus begins to unfold during the first 250 ps. This is followed by partial denaturation of α-helix1; the α-helix1 lost one turn at the N-terminal about 1000 ps and moved away from the rest of the protein during most of the simulation. α-helix3 was unfolded in the beginning of the simulation. In the β-sheet, its disruption starts at the edges of the β-sheet and near the irregular element of the β-sheet1 (β-bulge at residues 53 and 54); it is promoted by an increase in the twist and an influx of water molecules, and the β-sheet1 lost mostly about 850 ps. Actually, from the time dependence of secondary structure as well as the overlap view of tertiary structure (Figure 3), we can find that central β-sheet3 seems to be the most stable in five β-sheets. Meanwhile, rearrangements of secondary structural elements were observed along with the simulation and additional α-helix and β-sheet developed in the structure. However, these nonnative interactions were not stable enough as they would be disrupted over time. When the temperature was raised to 600 K, unfolding began almost immediately. Destruction of the native protein structure occurred very fast, as also indicated by the RMSF and RMSD values (Figure 1). The protein was highly coiled in the early stage simulation, and only a few rearrangement secondary structures remained in the protein. The MD trajectories can provide a detailed view of the conformational transitions in the early stage of thermal unfolding. Figure 3 shows snapshots of protein Barnase from the different temperature trajectories. The corresponding snapshots at 300 K are also shown as references.

3.4. Intramolecular Contacts

3.4.1. Hydrogen Bonding Pattern

Hydrogen bond is one of the factors influencing the thermal stability of protein. In the hydrogen bond calculations a distance cutoff of 3.0 Å and an angle cutoff of 20° were applied. The average numbers of hydrogen bonds are 29, 24, 17, and 18 for the 300 K, 400 K, 500 K, and 600 K simulations, respectively. Thus, as the simulation temperature is increased, there is a concomitant decrease in number of intact hydrogen bonds. This is reasonable as the structures become more distorted as the simulation temperature is raised. It is also evident from the plot (Figure 4) that although the number of hydrogen bonds varies in different temperature it is steadily maintained throughout the simulations except for 600 K simulation. The interesting finding is the rapid increase of hydrogen bonding number along with the simulation at 600 K after 94 ps. Due to a large distortion of regular secondary structural elements and unpacking of the hydrophobic cores, some of water molecules are inserted in hydrophobic cores and participate in hydrogen bonds as both donors and acceptors with the main-chain polar groups. This observation suggests the gradual destabilization of the protein in concert with increasing thermal stress.

(a)

(b)

(c)

(d)

3.4.2. Salt Bridge Analysis

To further probe the stability behavior of Barnase under thermal stress, we analyzed another important intramolecular contact, namely, salt bridge. Charged residues in globular proteins frequently form salt bridges. The electrostatic contribution of salt bridge has been suggested to be important for protein stability. Furthermore, the statistical analysis of salt bridges from mesophilic and thermophilic organisms has shown a higher frequency of complex salt bridges in thermophilic proteins, suggesting that they have a special role in thermostabilization.

In the structure of Barnase, nine salt bridges, Asp8-Arg110, Asp12-Arg110, Asp54-Lys27, Glu73-Lys27, Asp75-Arg83, Asp75-Arg87, Asp93-Arg69, Asp93-Lys66, and Glu60-Lys62, can be identified with the help of VMD. Interestingly, among these salt bridges, there are four salt bridge networks. Networks of ionic interactions occur when more than two ionic residues interact, and an increased occurrence has been suggested to be essential in explaining the enhanced thermal stability of protein [46]. In order to estimate the behavior of unfolding under thermal stress, the lifetime/occupancy of these pairs were analyzed in detail. In dynamic simulation at 300 K, these bridges were found to be stable during the period of 2.0 ns. Figure 5 shows the distance as a function of time in protein unfolding trajectory among these salt bridges. From the profiles it was found how changes in these salt bridge network interactions are correlated with the initial events associated with tertiary structure unfolding in response to thermal stress.

(a)

(b)

(c)

(d)

(e)

Figure 5

Oxygen-nitrogen distance of salt bridges plots. Time evolution of distance between selected atoms during simulation. (a) Distance between Asp8-Arg110 and Arg110-Asp12, which form salt bridge network in the core1. (b) Distance between Asp54-Lys27 and Lys27-Glu73; these residues form salt bridges close to the core2. (c) Distance between Arg83-Asp75 and Asp75-Arg87, (d) Lys66-Asp93 and Asp93-Arg69, and (e) Glu60-Lys62; these above-mentioned atom pairs form salt bridge close to the core3.

The salt bridge network of Asp8-Arg110-Asp12 is located in the main hydrophobic core1, which is formed by the packing of α-helix1 against the β-sheet and it is thought to be the major stabilizing element of Barnase. In hydrophobic core1, the two residues Asp8 and Asp12 located in α-helix1 could form salt bridge with Arg110 located in C-terminal. The C-terminal is docked to α-helix1 in the simulation, and this stable docking is dominated by strong electrostatic interactions between Arg110 and two acidic residues on the helix, Asp8 and Asp12. An interesting finding is that the Asp8-Arg110-Asp12 double salt bridge is not stable very much in the native structure of Barnase. As shown in Figure 5(a), the salt bridge network of Asp8-Arg110-Asp12 may be kept within a short distance in solution after 600 ps at 400 K. When the temperature is increased to 500 K, the salt bridge Asp8-Arg110 was maintained within a relatively short distance during the 1750 ps simulation, while the salt bridge Asp12-Arg110 was maintained within a relatively short distance only about 400 ps during the whole simulation, and the two residues fell apart eventually after 600 ps.

The rupture of the double salt bridge initiates the separation of the α-helix1 and β-sheet. The side chain of lle109 has moved away from the aromatic ring of Phe7 and the Asp8, Asp12. Then, some water moves into the center of the core; the inward motion of Lys98 is coupled to the Arg110 outward motion of the side chain. The exterior strands of the β-sheet were solvated by water molecules that replace some of the hydrogen bonds between β-sheets 4 and 5. The break of two salt bridges lead to the increase in accessible surface area and partial penetration of the water molecules and thus core1 undergoes a partial opening.

In the 550 K simulation, ruptures and restorations of the salt bridge Asp8-Arg110 were observed along the first half of unfolding simulations, and the Asp8 and Arg110 side chains begin to recover during the 1110 to 1890 ps period. The separation of the Asp12 and Arg110 side chains begins at initial stage of simulation, and the side chains of two residues came within a relatively short distance during the 381 to 1005 ps and 1099 to 1500 ps period of the simulation, but they fell apart eventually. The denaturation of the N-terminal part of α-helix1 (Phe7), the unfolding of the edges of the β-sheet, the denaturation of the C-terminal (Ile109), and the separating motion of loop1 (which contributes the Leu20 and Tyr24 sidechains to core1) contribute significantly to the unfolding of the main hydrophobic core1. Meanwhile, accompanying the solvation of hydrophobic core upon thermal unfolding, α-helix1 and strands of the β-sheet also undergo dramatic structural distortion changes again.

The salt bridge network of Asp54-Lys27-Glu73 was found to be stable at 300 K and 400 K simulations (Figure 5(b)), and they can hold the stable conformation of the spatial neighborhood. When the temperature rose to 500 K, the salt bridge of Asp54-Lys27 was maintained within 4 Å except for transient separation at about 1344 ps, while another salt bridge of Lys27-Glu73 only maintained the short distance from 400 ps to 600 ps and they were clearly separated from each other after about 650 ps. At 550 K, Glu73 could sometimes interact with Arg83 located in loop4 and Arg87 located in β-sheet3 during the simulation. The disruption of salt bridge Lys27-Glu73 means the loss of interaction between β-sheet2 and α-helix2; then, the sliding motion of loop4 is coupled to the movement of β-sheet2, which results in a slight increase in the β-sheet twist. Ala30, Leu42, and Ile51 side chains move away from the center of core2. The above-mentioned movement and the loss of native hydrogen bonds between β-sheet1 and β-sheet2 result in the separation of core2 from the rest of the protein. Concomitant with these changes, there was an obviously unfolding in the hydrophobic core2 (Figure 3).

The double salt bridges between Arg83, Arg87, and Asp75 were found to be very stable, up to 500 K simulation (Figure 5(c)). This salt bridge network is close to the hydrophobic core3 and provides a moderate stable spike holding the conformation of loop4, β-sheet2, and β-sheet3. Moreover, during the course of thermal simulation at 550 K, the salt bridge of Asp75-Arg87 located in the middle of β-sheet2 and β-sheet3, as expected, was found to be more stable than the salt bridge of Asp75-Arg83, which is located in between β-sheet2 and loop4. The two residues Asp75 and Arg83 came within a short distance during the first 1400 ps simulation; later they were completely separated. Since the presence of salt bridge network Asp75-Arg83-Arg87 and interstrand hydrogen bonds, simulation suggests that the most stable interaction appears in β-sheet2 and β-sheet3.

Another salt bridge network Lys66-Asp93-Arg69 is also close to the core3. In the network, the two residues Lys66 and Arg69 are all located in loop3, and the residue Asp93 is located in the β-turn connecting strands3 and 4. It is evident from the plot that salt bridge Asp93-Arg69 is maintained throughout the simulation (within 4.0 Å) without any transient separation from 300 K to 500 K (Figure 5(d)). Increasing the temperature to 550 K results in Asp93-Arg69 side chain contact breaks after 1750 ps and does not reform. In contrast, the stability of Lys66-Asp93 is not very stable at high temperature. The Lys66 side chain transiently separates from the Asp93 side chain during about the 1300 to 1400 ps period at 400 K simulation. However, the occupied time of simulation for this salt bridge descended to about 50% with an average distance of 5 Å at 500 K. The rupture of Lys66-Asp93 and several hydrogen bonds between β-sheet3 and β-sheet4 are coupled to the separating motion of β-sheet4 and β-sheet5 from the rest of the β-sheet, and the packing was looser than in the native state (Figure 3). In the 550 K simulations Lys66-Asp93 is stable only for 500 ps or so during the initial simulation. The presence of salt bridge networks Arg83-Asp75-Arg87 and Lys66-Asp93-Arg69 has stabilized the native state of hydrophobic core3 at elevated temperatures. Compared with other hydrophobic cores, core3 is the last one to unfold. For the stability of core3 against thermal unfolding, from the time dependence of secondary structure as well as the overlap view of tertiary structure (Figures 2 and 3), we can also find that the central three β-sheets, β2, β3, and β4, seem to be the most stable in protein, which provide the key stabilization interactions within hydrophobic core.

Finally, there is salt bridge Glu60-Lys62 located in the outer domain of hydrophobic core3, and the two residues are all located in loop3. Loop residues showed increased mobility relative to sheet or helical residues (Figure 1(b)). The salt bridge of Glu60-Lys62 could not be well-formed as it was disrupted in the beginning and recovered for a while in the middle of simulation along the four temperature simulations (Figure 5(e)). On the other hand, Glu60 preferred to form salt bridge with the neighboring residue Arg59. There occurred a positioning switch of the side chain of Glu60 from Lys62 to Arg59 in the simulation, leading to the departure of Glu60 from Lys62. Meanwhile, due to the fluctuation of the loop region, the backbone conformation and the tertiary packing of the loop were considerably influenced.

On the basis of these results, we concluded that surface salt bridge does stabilize the native state of the protein at elevated temperatures, and the complex salt bridge contribution to the overall protein stability is more than the individual pairs. Several experiments also confirm that hyperthermophilic proteins generally possess not only an increased number of surface salt bridges, but also an increased number of salt bridge networks [15, 47–49].

4. Conclusion

We have here investigated the electrostatic stability of noncovalent interactions in the context of temperature adaptation of Barnase by using MD simulations. The results show that hydrogen bond is very sensitive to heat, while salt bridge is comparatively stable. Nine salt bridges have been identified (Asp8-Arg110, Asp12-Arg110, Asp54-Lys27, Glu73-Lys27, Asp75-Arg83, Asp75-Arg87, Asp93-Arg69, Asp93-Lys66, and Glu60-Lys62) as critical salt bridges. These salt bridges are of fundamental importance in maintaining the structural integrity of the protein structure. Among these nine pairs, two salt bridge networks (Arg83-Asp75-Arg87 and Lys66-Asp93-Arg69) have been found to be extremely stable throughout the simulation up to 500 K. Two salt bridge networks located in core3 outer domain add more stability towards thermostable core3 region. The strength and the number of salt bridges present in a protein and whether they are involved in networks or not are important for the overall structural stability. Their presence can have a large impact on the structural integrity modulating molecular plasticity. The present study attempts to gain a deeper understanding of the noncovalent intramolecular interaction factors conferring thermostability of Barnase. It can offer a general picture of the first steps of unfolding and may help to design biotechnologically improved thermostable proteins.

Acknowledgments

This work was supported in part by the Fundamental Research Funds for the Central Universities under Grant no. JUSRP111A46 and the National Natural Science Foundation of China under Grant no. 61170119. Ming Li thanks the supports in part by the National Natural Science Foundation of China under the Project Grant nos. 61272402, 61070214, and 60873264.

References

S. Robic, M. Guzman-Casado, J. M. Sanchez-Ruiz, and S. Marqusee, “Role of residual structure in the unfolded state of a thermophilic protein,” Proceedings of the National Academy of Sciences of the United States of America, vol. 100, no. 20, pp. 11345–11349, 2003.
View at: Publisher Site | Google Scholar
A. Razvi and J. M. Scholtz, “Lessons in stability from thermophilic proteins,” Protein Science, vol. 15, no. 7, pp. 1569–1578, 2006.
View at: Publisher Site | Google Scholar
S. Kumar and R. Nussinov, “How do thermophilic proteins deal with heat?” Cellular and Molecular Life Sciences, vol. 58, no. 9, pp. 1216–1233, 2001.
View at: Google Scholar
K. A. Luke, C. L. Higgins, and P. Wittung-Stafshede, “Thermodynamic stability and folding of proteins from hyperthermophilic organisms,” FEBS Journal, vol. 274, no. 16, pp. 4023–4033, 2007.
View at: Publisher Site | Google Scholar
M. M. Gromiha, M. C. Pathak, K. Saraboji, E. A. Ortlund, and E. A. Gaucher, “Hydrophobic environment is a key factor for the stability of thermophilic proteins,” Proteins, vol. 81, no. 4, pp. 715–721, 2013.
View at: Publisher Site | Google Scholar
A. Szilágyi and P. Závodszky, “Structural differences between mesophilic, moderately thermophilic and extremely thermophilic protein subunits: results of a comprehensive survey,” Structure, vol. 8, no. 5, pp. 493–504, 2000.
View at: Publisher Site | Google Scholar
S. Basu and S. Sen, “Do homologous thermophilic-mesophilic proteins exhibit similar structures and dynamics at optimal growth temperatures? A molecular dynamics simulation study,” Journal of Chemical Information and Modeling, vol. 53, no. 2, pp. 423–434, 2013.
View at: Publisher Site | Google Scholar
A. D. Meruelo, S. K. Han, S. Kim, and J. U. Bowie, “Structural differences between thermophilic and mesophilic membrane proteins,” Protein Science, vol. 21, no. 11, pp. 1746–1753, 2012.
View at: Google Scholar
E. Papaleo, M. Olufsen, L. de Gioia, and B. O. Brandsdal, “Optimization of electrostatics as a strategy for cold-adaptation: a case study of cold- and warm-active elastases,” Journal of Molecular Graphics and Modelling, vol. 26, no. 1, pp. 93–103, 2007.
View at: Publisher Site | Google Scholar
A. V. Gribenko, M. M. Patel, J. Liu, S. A. McCallum, C. Wang, and G. I. Makhatadze, “Rational stabilization of enzymes by computational redesign of surface charge-charge interactions,” Proceedings of the National Academy of Sciences of the United States of America, vol. 106, no. 8, pp. 2601–2606, 2009.
View at: Publisher Site | Google Scholar
J. M. Vinther, S. M. Kristensen, and J. J. Led, “Enhanced stability of a protein with increasing temperature,” Journal of the American Chemical Society, vol. 133, no. 2, pp. 271–278, 2011.
View at: Publisher Site | Google Scholar
P. Dürrschmidt, J. Mansfeld, and R. Ulbrich-Hofmann, “Differentiation between conformational and autoproteolytic stability of the neutral protease from Bacillus stearothermophilus containing an engineered disulfide bond,” European Journal of Biochemistry, vol. 268, no. 12, pp. 3612–3618, 2001.
View at: Publisher Site | Google Scholar
C. N. Pace, H. Fu, K. L. Fryar et al., “Contribution of hydrophobic interactions to protein stability,” Journal of Molecular Biology, vol. 408, no. 3, pp. 514–528, 2011.
View at: Publisher Site | Google Scholar
G. I. Makhatadze, V. V. Loladze, D. N. Ermolenko, X. Chen, and S. T. Thomas, “Contribution of surface salt bridges to protein stability: guidelines for protein engineering,” Journal of Molecular Biology, vol. 327, no. 5, pp. 1135–1148, 2003.
View at: Publisher Site | Google Scholar
J. H. Missimer, M. O. Steinmetz, R. Baron et al., “Configurational entropy elucidates the role of salt-bridge networks in protein thermostability,” Protein Science, vol. 16, no. 7, pp. 1349–1359, 2007.
View at: Publisher Site | Google Scholar
A. S. Thomas and A. H. Elcock, “Molecular simulations suggest protein salt bridges are uniquely suited to life at high temperatures,” Journal of the American Chemical Society, vol. 126, no. 7, pp. 2208–2214, 2004.
View at: Publisher Site | Google Scholar
J. Chen, H. M. Yu, C. C. Liu, J. Liu, and Z. Y. Shen, “Improving stability of nitrile hydratase by bridging the salt-bridges in specific thermal-sensitive regions.,” Journal of Biotechnology, vol. 164, pp. 354–362, 2012.
View at: Publisher Site | Google Scholar
J. Li, Y. Matsumura, M. Shinjo, M. Kojima, and H. Kihara, “A stable α-helix-rich intermediate is formed by a single mutation of the β-sheet protein, src SH3, at pH 3,” Journal of Molecular Biology, vol. 372, no. 3, pp. 747–755, 2007.
View at: Publisher Site | Google Scholar
S. P. He, Z. T. Ding, and F. Liu, “Output regulation of a class of continuous-time Markovian jumping systems,” Signal Processing, vol. 93, no. 2, pp. 411–419, 2013.
View at: Publisher Site | Google Scholar
S. He and F. Liu, “Robust stabilization of stochastic Markovian jumping systems via proportional-integral control,” Signal Processing, vol. 91, no. 11, pp. 2478–2486, 2011.
View at: Publisher Site | Google Scholar | Zentralblatt MATH
S. P. He and F. Liu, “Finite-time H∞ control of nonlinear jump systems with time-delays via dynamic observer-based state feedback,” IEEE Transactions on Fuzzy Systems, vol. 20, no. 4, pp. 605–614, 2012.
View at: Publisher Site | Google Scholar
C. Cattani, A. Ciancio, and B. Lods, “On a mathematical model of immune competition,” Applied Mathematics Letters, vol. 19, no. 7, pp. 678–683, 2006.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
C. Cattani and A. Ciancio, “Separable transition density in the hybrid model for tumor-immune system competition,” Computational and Mathematical Methods in Medicine, vol. 2012, Article ID 610124, 6 pages, 2012.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
G. Strittmatter, J. Janssens, C. Opsomer, and J. Botterman, “Inhibition of fungal disease development in plants by engineering controlled cell death,” BioTechnology, vol. 13, no. 10, pp. 1085–1089, 1995.
View at: Google Scholar
S. Agarwal, B. Nikolai, T. Yamaguchi, P. Lech, and N. V. Somia, “Construction and use of retroviral vectors encoding the toxic gene barnase,” Molecular Therapy, vol. 14, no. 4, pp. 555–563, 2006.
View at: Publisher Site | Google Scholar
S. Leuchtenberger, A. Perz, C. Gatz, and J. W. Bartsch, “Conditional cell ablation by stringent tetracycline-dependent regulation of barnase in mammalian cells,” Nucleic Acids Research, vol. 29, no. 16, article e76, 2001.
View at: Publisher Site | Google Scholar
J. L. Ivanova, E. F. Edelweiss, O. G. Leonova, T. G. Balandin, V. I. Popenko, and S. M. Deyev, “Application of fusion protein 4D5 scFv-dibarnase:barstar-gold complex for studying P185HER2 receptor distribution in human cancer cells,” Biochimie, vol. 94, no. 8, pp. 1833–1836, 2012.
View at: Publisher Site | Google Scholar
T. Konuma, T. Kimura, S. Matsumoto et al., “Time-resolved small-angle X-ray scattering study of the folding dynamics of barnase,” Journal of Molecular Biology, vol. 405, no. 5, pp. 1284–1294, 2011.
View at: Publisher Site | Google Scholar
J. Yin, D. Bowen, and W. M. Southerland, “Barnase thermal titration via molecular dynamics simulations: detection of early denaturation sites,” Journal of Molecular Graphics and Modelling, vol. 24, no. 4, pp. 233–243, 2006.
View at: Publisher Site | Google Scholar
A. Caflisch and M. Karplus, “Molecular dynamics simulation of protein denaturation: solvation of the hydrophobic cores and secondary structure of barnase,” Proceedings of the National Academy of Sciences of the United States of America, vol. 91, no. 5, pp. 1746–1750, 1994.
View at: Google Scholar
A. Li and V. Daggett, “Molecular dynamics simulation of the unfolding of barnase: characterization of the major intermediate,” Journal of Molecular Biology, vol. 275, no. 4, pp. 677–694, 1998.
View at: Publisher Site | Google Scholar
A. Langham and Y. N. Kaznessis, “Molecular simulations of antimicrobial peptides,” Methods in Molecular Biology, vol. 618, pp. 267–285, 2010.
View at: Publisher Site | Google Scholar
T. Darden, D. York, and L. Pedersen, “Particle mesh Ewald: an $N \cdot \log (N)$ method for Ewald sums in large systems,” The Journal of Chemical Physics, vol. 98, no. 12, pp. 10089–10092, 1993.
View at: Publisher Site | Google Scholar
T. Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Interdisciplinary Applied Mathematics, Springer, New York, NY, USA, 2nd edition, 2010.
View at: Publisher Site | MathSciNet
R. Kubo, “The fluctuation-dissipation theorem,” Reports on Progress in Physics, vol. 29, no. 1, article 306, pp. 255–284, 1966.
View at: Publisher Site | Google Scholar | Zentralblatt MATH
W. G. Hoover, “Canonical dynamics: equilibrium phase-space distributions,” Physical Review A, vol. 31, no. 3, pp. 1695–1697, 1985.
View at: Publisher Site | Google Scholar
S. Nose, “A molecular dynamics method for simulations in the canonical ensemble,” Molecular Physics, vol. 52, no. 2, pp. 255–268, 1984.
View at: Publisher Site | Google Scholar
G. J. Martyna, D. J. Tobias, and M. L. Klein, “Constant pressure molecular dynamics algorithms,” The Journal of Chemical Physics, vol. 101, no. 5, pp. 4177–4189, 1994.
View at: Google Scholar
D. Frenkel and B. Smit, Understanding Molecular Simulations: From Algorithms to Applications, Academic Press, San Diego, Calif, USA, 2nd edition, 2002.
View at: Zentralblatt MATH
L. Kalé, R. Skeel, M. Bhandarkar et al., “NAMD2: greater scalability for parallel molecular dynamics,” Journal of Computational Physics, vol. 151, no. 1, pp. 283–312, 1999.
View at: Publisher Site | Google Scholar
A. D. MacKerell Jr., D. Bashford, M. Bellott et al., “All-atom empirical potential for molecular modeling and dynamics studies of proteins,” Journal of Physical Chemistry B, vol. 102, no. 18, pp. 3586–3616, 1998.
View at: Google Scholar
S. Baudet and J. Janin, “Crystal structure of a barnase-d(GpC) complex at 1.9 Å resolution,” Journal of Molecular Biology, vol. 219, no. 1, pp. 123–132, 1991.
View at: Publisher Site | Google Scholar
W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, “Comparison of simple potential functions for simulating liquid water,” The Journal of Chemical Physics, vol. 79, no. 2, pp. 926–935, 1983.
View at: Google Scholar
W. Humphrey, A. Dalke, and K. Schulten, “VMD: visual molecular dynamics,” Journal of Molecular Graphics, vol. 14, no. 1, pp. 33–38, 1996.
View at: Publisher Site | Google Scholar
J.-P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, “Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes,” Journal of Computational Physics, vol. 23, no. 3, pp. 327–341, 1977.
View at: Google Scholar
M. Olufsen, E. Papaleo, A. O. Smalås, and B. O. Brandsdal, “Ion pairs and their role in modulating stability of cold- and warm-active uracil DNA glycosylase,” Proteins, vol. 71, no. 3, pp. 1219–1230, 2008.
View at: Publisher Site | Google Scholar
G. Pappenberger, H. Schurig, and R. Jaenicke, “Disruption of an ionic network leads to accelerated thermal denaturation of D-glyceraldehyde-3-phosphate dehydrogenase from the hyperthermophilic bacterium Thermotoga maritima,” Journal of Molecular Biology, vol. 274, no. 4, pp. 676–683, 1997.
View at: Publisher Site | Google Scholar
K. S. P. Yip, T. J. Stillman, K. L. Britton et al., “The structure of Pyrococcus furiosus glutamate dehydrogenase reveals a key role for ion-pair networks in maintaining enzyme stability at extreme temperatures,” Structure, vol. 3, no. 11, pp. 1147–1158, 1995.
View at: Google Scholar
N. Dhaunta, K. Arora, S. K. Chandrayan, and P. Guptasarma, “Introduction of a thermophile-sourced ion pair network in the fourth beta/alpha unit of a psychophile-derived triosephosphate isomerase from Methanococcoides burtonii significantly increases its kinetic thermal stability,” Biochimica et Biophysica Acta, vol. 1834, no. 6, pp. 1023–1033, 2013.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2013 Zhiguo Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

2439

Downloads

1619

Citations