J. Phys. Chem. B 2008, 112, 15127–15133
15127
Dynamic Redox Environment-Intensified Disulfide Bond Shuffling for Protein Refolding in Vitro: Molecular Simulation and Experimental Validation Diannan Lu† and Zheng Liu* Department of Chemical Engineering, Tsinghua UniVersity, Beijing 10084, China, ReceiVed: May 26, 2008; ReVised Manuscript ReceiVed: September 13, 2008
One challenge in protein refolding is to dissociate the non-native disulfide bonds and promote the formation of native ones. In this study, we present a coarse-grained off-lattice model protein containing disulfide bonds and simulate disulfide bond shuffling during the folding of this model protein. Introduction of disulfide bonds in the model protein led to enhanced conformational stability but reduced foldability in comparison to counterpart protein without disulfide bonds. The folding trajectory suggested that the model protein retained the two-step folding mechanism in terms of hydrophobic collapse and structural rearrangement. The disulfide bonds located in the hydrophobic core were formed before the collapsing step, while the bonds located on the protein surface were formed during the rearrangement step. While a reductive environment at the initial stage of folding favored the formation of native disulfide bonds in the hydrophobic core, an oxidative environment at a later stage of folding was required for the formation of disulfide bonds at protein surface. Appling a dynamic redox environment, that is, one that changes from reductive to oxidative, intensified disulfide bond shuffling and thus resulted in improved recovery of the native conformation. The above-mentioned simulation was experimentally validated by refolding hen-egg lysozyme at different urea concentrations and oxidized glutathione/reduced glutathione (GSSG/GSH) ratios, and an optimal redox environment, in terms of the GSSG to GSH ratio, was identified. The implementation of a dynamic redox environment by tuning the GSSG/GSH ratio further improved the refolding yield of lysozyme, as predicted by molecular simulation. Introduction Disulfide bond is important for maintaining the native conformation and biological functions of many proteins.1,2 De novo introduction of disulfide bonds by genetic manipulation is an effective means of enhancing protein stability. However, promotion of disulfide shuffling, an essential procedure in protein folding in vivo and in vitro that results in the formation of native disulfide bonds and dissociation of non-native ones, remains inadequately understood, particularly at a molecular level. As pointed out by Anfinsen in the 1970s,3 when a suitable environment is created, a fully reduced and unfolded protein can fold into its native conformation, which is characterized by a global state of minimum free energy and the formation of all native disulfide bonds. On the other hand, it has been experimentally demonstrated that the special complexities in the folding of disulfide bond-containing proteins in terms of the presence of diversified folding pathways, various oxidative intermediate states, and consequently final products.4,5 The formation and breakage of disulfide bonds is known to be ubiquitous; however, the influence of this process on oxidative folding is not well understood. To date, two models have been proposed to explain the mechanism of protein folding.6,7 One is the folded-precursor model, which attributes the formation of native disulfide bonds to a consequential event of the folding process. The other is the quasi-stochastic mechanism, which assumes that the formation of disulfide bonds, both native and non-native, begins from the initially unfolded conformations. * To whom correspondence should be addressed. E-mail: liuzheng@ mail.tsinghua.edu.cn. Tel: 86-10-6277 9876. Fax: 86-10-6277 0304. † E-mail:
[email protected].
There remain experimental observations that cannot be interpreted by either of these models. Decreased entropy in the unfolded state has been identified as the main reason for the increased stability of disulfide-bonded proteins.8-10 According to the proximity rule proposed by Camacho and Thirumulai,9 the effect of entropy on protein stability becomes stronger with increased sequence separation of the pair of linked residue. The individual protein can deviate significantly from this ideal behavior, as proven by the difficulties encountered in engineering disulfide bond-stabilized proteins. In practice, in comparison to the absence of a disulfide bond, the presence of a disulfide bond in the folded state complicates the folding process. Most efforts from molecular simulation of the folding of protein with disulfide bonds focus on how the oxidative or reductive cysteine affects protein stability. In these studies, the cysteine residues are set at the beginning and remain fixed throughout the folding process.11-13 Recently, Camacho and Thirumalai14 have attempted to explore dynamic disulfide bond formation using a 23-bead 2D lattice model with Monte Carlo methods and have shown significant effects of the native and non-native folding intermediates on the folding pathways. More recently, Chinchio et al.15 have developed a “UNRES” force field to describe the dynamic formation and breakage of disulfide bonds, in which the potential function is obtained by introducing a transition barrier that separates the bonded and nonbonded states of half-cystine residues. An improved prediction has been achieved for helical proteins in which on one hand, the formation of a disulfide bond stabilizes the folding intermediate but on the other hand, it also reduces the folding rate. In vivo folding experiments have established that in the early stage of folding, disulfide bond formation is error-prone, that
10.1021/jp804649g CCC: $40.75 2008 American Chemical Society Published on Web 10/29/2008
15128 J. Phys. Chem. B, Vol. 112, No. 47, 2008
Lu and Liu long-range van der Waals attractions. The bond energy is related to the fluctuation of bond lengths, bond angles, dihedral angles, and disulfide bonds. Specifically, bond fluctuation is described by the harmonic potential as
Vb )
∑ kb(r - σ)2
(1)
bonds
Figure 1. Native structure of the HPNC protein. Grey, hydrophobic; red, polar; blue, neutral; yellow, cysteine.
is, the wrong cysteines are connected16,17 or the correct cysteines are paired but in a temporal order that hinders folding.18-20 The slow spontaneous folding of disulfide bond-containing proteins in vitro is incompatible with the time scale of secretion in vivo.21,22 The eukaryotic cell overcomes this problem by using a specialized redox environment, that is, the environment in the endoplasmic reticulum (ER) that is equipped with a biocatalyst such as protein disulfide isomerase (PDI) for disulfide bond formation and isomerization.23-26 However, in vitro oxidative protein refolding remains empirical and optimal conditions for proteins with different disulfide bonds vary significantly from one to the another one.27-33 This study is an extension of the authors’ efforts in establishing a dynamic environment facilitating protein folding, which mimics the molecular chaperone that recognizes and captures the misfolded intermediates while promotes the formation of the native ones.34 A model protein containing two disulfide bonds was proposed, in which one disulfide bond is located in the core while the other is at the surface of the protein. Morse potential was applied to describe the dynamic formation and breakage of disulfide bonds during Langevin dynamics simulation of disulfide bond shuffling. The effects of redox conditions on protein stability, folding kinetics, and yield of native conformation were examined as function of denaturant concentrations. Dynamic control of the redox condition to promote disulfide bond shuffling, which mimicked PDI in vivo, resulted in improved oxidative folding in comparison to that obtained under steady-state operation. The above-mentioned simulations were experimentally validated using hen-egg lysozyme, a protein containing four pairs of disulfide bonds. A dynamic redox environment (from reductive to oxidative) was created by tuning the GSSG/GSH ratio and improved folding was observed, as predicted by molecular simulation. Models and Methods Simulations. HPNC Protein Model. In this study, a modified HPN model, termed the HPNC model, is proposed (shown in Figure 1) in which all residues are treated as spherical beads of equal size that are tangentially connected in the sequence H6C1H2N3(PH)4N3C1H5C1H2N3(PH)4L1C1L1. In comparison to the HPN model protein that has been extensively applied in the simulation of protein folding and aggregation,35-40 the HPNC model contains 46 beads that are divided into four groups and are denoted as hydrophobic (H), polar (P), neutral (N), and cysteine (C). The side chains of the amino-acid residues are not explicitly considered. The native structure of the model protein (shown in Figure 1) was obtained by a simulated annealing method. According to previous reports, the Hamiltonian of the model protein includes the bond energy, excluded volume effects, and
where r is the center-to-center distance between two nearest neighboring amino-acid residues, σ is the equilibrium bond length, and kb ) 100εh is the bond spring constant. For a tangentially connected chain, the equilibrium bond length is the same as the bead diameter. In this study, εh and σ are 5.0 kJ/ mol and 0.38 nm, respectively. The bending potential specifies the energy associated with the fluctuation of the bond angles. It is also expressed in a harmonic form as
Vθ )
∑
kθ(θ - θ0)2
(2)
bondangles
where θ denotes the bond angle, θ0 ) 105°, is the equilibrium bond angle, and kθ ) 10εh/(rad)2. The dihedral potential is described by
Vφ )
∑ [kφ(1)(1 + cos φ) + kφ(2)(1 + cos 3φ)]
dihedralangles
(3) where kφ(1) ) 0 and kφ(2) ) 0.2εh if more than one of the four beads defining the dihedral angle φ is neutral (N), and kφ(1) ) kφ(2) ) 1.2εh otherwise. The dihedral potential exhibits a minimum in the trans or gauche configurations. These energy minima are responsible for the formation of β-turns in the native structure of the model protein. The Morse potential is introduced here to describe the formation and breakage of disulfide bonds, which only occurs between cysteines (C)
VD(r) ) Dij[1 - exp(-βij(r - σ))]2
(4)
where Dij is the depth of the well, which reflects the oxidative capacity of the folding environment. A lower Dij indicates a more reductive environment for protein folding and vice versa. βij defines the steepness of the well, which is fixed as 2.63 nm-1, and σ is the equilibrium distance of the disulfide bond, which is fixed as 0.38 nm. The “nonbonded” pair interaction includes the excluded volume effect and the long-range van der Waals attraction, which takes the form of Lennard-Jones (LJ)-like potential when applied to amino-acid residues separated by at least two peptide bonds
VLJ )
∑ 4εh[Aij( σr )
12
i,j
( σr ) ]
- ε*pBij
6
(5)
where εh represents unit energy, and the dimensionless parameters Aij and Bij depend on the identities of the interacting beads. To facilitate direct comparison with results from previous singlechain simulations, we adopted Aij ) 1 and Bij ) 1 for
Redox Environment-Intensified Disulfide Bond Shuffling in Vitro interactions between HH pairs; Aij ) 1/3 and Bij ) -1 for PH, PC, HC, CC, and PP; and Aij ) 1 and Bij ) 0 for the HN, PN, CN, and NN pairs. The dimensionless parameter ε*p reflects the hydrophobic interaction strength between hydrophobic beads of the model protein. High εp* implies a strong hydrophobic interaction between protein beads. In practice, ε*p can be tuned by the denaturant concentration, and a high ε*p value indicates a low denaturant concentration and vice versa. Simulation Method. Langevin dynamics and the VelocityVerlet algorithm were applied to examine structural transitions in different redox environments. For direct comparison with previous simulations, the friction coefficient was set as γ ) 0.05. The protein configuration is updated in time step of 0.0002 ps. As in the previous case of a single confined protein,41 velocity Langevin dynamics coupled with the Velocity-Verlet algorithm was performed using Gromacs 3.3 as the platform.42 The force function is
mi
d2ri 2
dt
) -miξi
dri + Fi(r) + ri◦ dt
(6)
where mi ) 2.2 × 10-22 g/mol is the average mass of an amino acid, ri is the position of bead i, ξi ) 1/(5.0 ps) is the friction constant, Fi(r) is the configurational force, and r°i is the random force that is related to the temperature through the fluctuationdissipation theorem. The weighted histogram analysis method (WHAM) was used to characterize the thermodynamics of HPNC model protein.36 Trajectories for WHAM analysis were carried out at different ε*p from 0.0 to 1.4 in increments of 0.02 and 275 K. At each value of ε*p, isolated HPNC protein was equilibrated for 8 ns, followed by sampling for 2 ns. One hundred such cycles were carried out at each ε*p to calculate the ensemble averages. Three key properties are calculated during the course of simulations: the total potential energy (V), protein’s radius of gyration (Rg), and structure overlap function (χ). The total potential energy includes the contributions from the bond energy, long-range van der Waals interactions, and disulfide bond potential
V ) Vb + Vθ + Vφ + VD + VLJ
(7)
The radius of gyration is defined as
R2g
2 ) N(N - 1)
∑∑ i
