proteins STRUCTURE O FUNCTION O BIOINFORMATICS
Glucose oxidase from Penicillium amagasakiense: Characterization of the transition state of its denaturation from molecular dynamics simulations Guido Todde,1 Sven Hovm€ oller,1 Aatto Laaksonen,1,2,3 and Francesca Mocci1,2* 1 Department of Materials and Environmental Chemistry, Stockholm University, Stockholm, Sweden 2 Dipartimento di Scienze Chimiche e Geologiche, Cagliari University, Cagliari, Italy 3 Stellenbosch Institute of Advanced Studies (STIAS), Wallenberg Research Centre at Stellenbosch University, South Africa
ABSTRACT Glucose oxidase (GOx) is a flavoenzyme having applications in food and medical industries. However, GOx, as many other enzymes when extracted from the cells, has relatively short operational lifetimes. Several recent studies (both experimental and theoretical), carried out on small proteins (or small fractions of large proteins), show that a detailed knowledge of how the breakdown process starts and proceeds on molecular level could be of significant help to artificially improve the stability of fragile proteins. We have performed extended molecular dynamics (MD) simulations to study the denaturation of GOx (a protein dimer containing nearly 1200 amino acids) to identify weak points in its structure and in this way gather information to later make it more stable, for example, by mutations. A denaturation of a protein can be simulated by increasing the temperature far above physiological temperature. We have performed a series of MD simulations at different temperatures (300, 400, 500, and 600 K). The exit from the protein’s native state has been successfully identified with the clustering method and supported by other methods used to analyze the simulation data. A common set of amino acids is regularly found to initiate the denaturation, suggesting a moiety where the enzyme could be strengthened by a suitable amino acid based modification. Proteins 2014; 82:2353–2363. C 2014 Wiley Periodicals, Inc. V
Key words: protein denaturation; mutation; unfolding; MD simulation; cluster analysis.
INTRODUCTION Glucose oxidase (GOx) is a flavoenzyme catalyzing the oxidation of b-D-glucose to d-gluconolactone and the reduction of molecular oxygen to hydrogen peroxide. It is not found active on a-D-glucose. GOx has important applications in food and medical industries. Prototypes of miniature enzymatic fuel cells using GOx as catalyst have already been produced.1 GOx is found in many organisms such as insects and fungi2,3 and the extraction and characterization of GOx from different sources is an active research area. For biotechnological applications it is extracted mainly from the Aspergillus and Penicillium fungi. To improve its catalytic properties and structural stability is an important area of bioengineering.4–9 Although Aspergillus niger is currently the most common source as well as the best characterized, GOx from Penicillium amagasakiense (GOx-pa) has been shown to catalyze the glucose oxidation more efficiently. This makes it a better candidate to be used in bioelectrochemical nanodevices.10
C 2014 WILEY PERIODICALS, INC. V
The GOx-pa is a homodimer formed by two identical chains (see Fig. 1). It belongs to the a/b-fold class according to the SCOP11 classification. The monomers of GOx-pa contain 587 amino acids each, and one prosthetic group, flavin adenine dinucleotide (FAD). In the crystallized protein (PDB)12 ID: 1GPE13 each monomer has a molecular weight of 66.7 kDa with four N-glycosylation sites at positions ASN93, ASN165, ASN357, and ASN392, respectively. Three different sugars are bound at these positions: N-acetyl-D-glucosamine, b-D-mannose, and a-D-mannose. Their binding order is shown in Supporting Information (Fig. S1).
Additional Supporting Information may be found in the online version of this article. *Correspondence to: Francesca Mocci, Dipartimento di Scienze Chimiche e Geologiche, Cagliari University, Cagliari, Italy. E-mail: [email protected] Received 12 December 2013; Revised 22 March 2014; Accepted 29 April 2014 Published online 9 May 2014 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/prot.24596
PROTEINS
2353
G. Todde et al.
Figure 1 Cartoon representation of GOx. Colors legend: monomers (ice-blue and blue), FAD (green), sugars (pink), and dimer interface (red). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Unfortunately both GOx-an and GOx-pa suffer from relatively short operational lifetimes when used in nanodevices due to their relatively fragile structure. To understand how the inactivation of this protein occurs, several experimental studies have focused on the thermally and/ or chemically induced denaturation of GOx in aqueous solution.14–17 However, results obtained in these experimental studies do not allow to identify at the residue level where the denaturation starts to take place. This would be of great help to propose protein modifications capable of strengthening the protein. Atomistic molecular dynamics (MD) simulations have been used in the past to study unfolding and to identify the possible pathways in denaturation processes.18,19 In computer simulations the denaturation is conveniently induced and speeded up by increasing the temperature, or using solvents including agents known to destabilize the protein structure like urea or guanadinium.18,20–24 The pressure influence on protein denaturation process has also been studied.25–27 Moreover, artificial forces acting on selected groups of atoms have been used to trigger a denaturation.28,29 Concerning the relationship between high-temperature unfolding and unfolding under native conditions, one very important observation from the simulations of thermal denaturation is that “increasing the temperature accelerates protein unfolding without changing the pathway”.21 This means that it can be assumed that the thermal unfolding follows an Arrhenius law behavior, and the increase in temperature can efficiently compensate the limited timescale accessible to MD simulations.
2354
PROTEINS
MD studies of protein denaturation have been comprehensively reviewed recently by Toofanny and Daggett.30 Furthermore previous denaturation studies31 have shown that in order to adequately sample the denaturation process, an order of 5–10 simulations is required. Understanding of how the breakdown process starts and proceeds on molecular level can help to artificially improve the stability of proteins.18,31–34 Information obtained from detailed denaturation studies have already led to a design of a more robust enzyme, by strengthening its weak regions and transition points. In fact, mutating key residues close to structural transition points along the denaturation process have been proven to be capable to strengthen the whole structure of the enzyme and/or to make the moiety around the catalytic center more stable. Yoda et al.32 studied goat lactoalbumin by MD simulations and Circular Dichroism (CD) experiments and found that a single amino acid substitution was capable of increasing the stability of the native state by about 3.5 kcal/mol. The structural disruptions observed in the simulations were found to be consistent with the experimental results. Also the results obtained by Pikkemaat et al.33 on haloalkane dehalogenase and of Badieyan et al.34 on a glycoside hydrolases show that MD simulations are capable of identifying the protein domains that can be strengthened. In this study, we simulate the thermally induced denaturation process of GOx-pa at different temperatures from room temperature up to 600 K and follow the denaturation pathways in a systematic way. To the best of our knowledge, our work is the first of its kind performed for such large protein as GOx-pa. The results from our simulations would allow us to rationally propose mutations to obtain thermally a more stable protein. Here we focus on determination and characterization of the transition state (TS) of the denaturation process. The full denaturation process of GOx will be discussed in our upcoming study. Although already shortly discussed here, we will propose amino acid based modifications to strengthen the enzyme to be used in bionanodevices in a later communication. MATERIALS AND METHODS All the MD simulations were performed using Gromacs 4.5.335–38 and the amber99sb39 force field. Parameters for the cofactor FAD that were not present in the Amber force field were taken, except for the charges, from the work of Schneider and Suhnel.40 The partial atomic charges were calculated following the method used in amber99sb, using the restrained electrostatic potential procedure41 at the Hartree–Fock theory level with 6–31G* basis set. During the geometry optimization, the cofactor was kept in an elongated conformation, close to that observed in the crystal structure, by freezing the relevant dihedral angles. The starting
MD Study of GOx Denaturation
structure used for all simulations was the crystal structure solved by X-ray diffraction at 1.80 A˚ resolution obtained from the Protein Data Bank12 (PDB ID: 1GPE13). The simulations of the dimer were performed at four different temperatures: 300, 400, 500, and 600 K. Although the commonly used force field parameters are expected to be valid at ambient temperatures, numerous studies show that the same parameters can be applied also when working at very much higher temperatures.18–21,30,34 We have performed three independent simulations at both 500 and 600 K starting from different structures, randomly chosen from the first 10 ns produced at 300 K. Since GOx is a homodimer, every simulation of the dimer, in fact consists of two simulations of the monomer. Hence, totally 14 simulations of the monomer were produced at high temperature. All simulations followed the same protocol. The equilibration procedure, described in the following, is similar to that used by Daggett and Beck.31 The potential energy of the starting structure was briefly minimized (500 steps) with a mixed conjugate-gradient and steepestdescent algorithm. The resulting structure was placed in an empty cubic box subsequently filled with 68,000 TIP3P42 water molecules. To make the system neutral 32 water molecules were replaced by 32 sodium ions. The system was then energy-minimized (water only) before a short MD simulation (10 ps) of the solvent was performed keeping the solute fixed. Then the protein was minimized keeping the water constrained and a final minimization of the entire system was performed. At this point a 20-ns MD simulation at 300 K in NPT ensemble was carried out. The system was heated linearly (20 K per ps) in all simulations in the NVT ensemble. The lengths of the production runs were 20 ns at 300 K, 500 ns at 400 K, 150 ns at 500 K, and 15 ns at 600 K. All production runs were carried out in the NVT ensemble (except at 300 K where NPT was used) using 1 fs as timestep and the trajectories were saved every 1 ps (10 ps at 400 K). The Nose-Hoover thermostat43,44 was applied in all simulations and the Parrinello–Rahman barostat45,46 in the 300 K simulation. In all simulations the van der Waals interactions were computed within a cutoff distance of 10 A˚ and the Coulombic interactions were treated with the PME47,48 method. Bonds involving hydrogen atoms were treated by LINCS algorithm.49 All simulated systems are composed of 225,000 particles. All the produced trajectories sum up for a total simulated time of 1015 ls. Converted to a single AMD Opteron 6172 processor, the total computational time used would have been about 580,000 CPU hours. Analysis of the structural modifications involved the calculation of contacts and hydrogen bonds (HBs) between residues. A contact is defined in the following manner: each heavy atom (C, N, O, and S) is considered on average equal to a volume of 19.5 A˚3,50 so the vol-
Figure 2 Ca root mean square fluctuation. (A) and (B) show RMSF at 300 K (violet), respectively, for the first and second monomer. The crystallographic B-factors13 are converted to RMS fluctuation by the equation RMSF 5 (3B/8p2)1/2 and are shown by the orange line. The segments on the top of each plot represent the position of b-strands (black) and a-helices (red), while the vertical blue lines represent the position of the glycosylation sites. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
ume of each amino acid is computed. An average volume of 162.8 A˚3 is obtained and this volume is assimilated to a sphere with a radius of 3.4 A˚. Hence, two nonconsecutive residues are considered in contact when the distance between their centers of mass is shorter than 7 A˚ (twice the mean sphere radius in order to avoid false counts from the second shell of neighbors). A HB is formed when the distance between donor (NH, OH) and acceptor (N,O) is shorter than 3.4 A˚ and the angle Acceptor-Donor-Hydrogen is smaller than 60 . The secondary structure (SS) is assigned according to the Dictionary of Protein Secondary Structure program.51,52
RESULTS AND DISCUSSION Reference simulation at 300 K
The simulation performed at 300 K was used as a reference to test the stability of the force field and of the adopted simulation protocol. During the 20 ns simulated at 300 K the protein was found to be stable. The Ca RMSD (see Fig. S2 in Supporting Information) of the single monomers was
Glucose oxidase from Penicillium amagasakiense: Characterization of the transition state of its denaturation from molecular dynamics simulations Guido Todde,1 Sven Hovm€ oller,1 Aatto Laaksonen,1,2,3 and Francesca Mocci1,2* 1 Department of Materials and Environmental Chemistry, Stockholm University, Stockholm, Sweden 2 Dipartimento di Scienze Chimiche e Geologiche, Cagliari University, Cagliari, Italy 3 Stellenbosch Institute of Advanced Studies (STIAS), Wallenberg Research Centre at Stellenbosch University, South Africa
ABSTRACT Glucose oxidase (GOx) is a flavoenzyme having applications in food and medical industries. However, GOx, as many other enzymes when extracted from the cells, has relatively short operational lifetimes. Several recent studies (both experimental and theoretical), carried out on small proteins (or small fractions of large proteins), show that a detailed knowledge of how the breakdown process starts and proceeds on molecular level could be of significant help to artificially improve the stability of fragile proteins. We have performed extended molecular dynamics (MD) simulations to study the denaturation of GOx (a protein dimer containing nearly 1200 amino acids) to identify weak points in its structure and in this way gather information to later make it more stable, for example, by mutations. A denaturation of a protein can be simulated by increasing the temperature far above physiological temperature. We have performed a series of MD simulations at different temperatures (300, 400, 500, and 600 K). The exit from the protein’s native state has been successfully identified with the clustering method and supported by other methods used to analyze the simulation data. A common set of amino acids is regularly found to initiate the denaturation, suggesting a moiety where the enzyme could be strengthened by a suitable amino acid based modification. Proteins 2014; 82:2353–2363. C 2014 Wiley Periodicals, Inc. V
Key words: protein denaturation; mutation; unfolding; MD simulation; cluster analysis.
INTRODUCTION Glucose oxidase (GOx) is a flavoenzyme catalyzing the oxidation of b-D-glucose to d-gluconolactone and the reduction of molecular oxygen to hydrogen peroxide. It is not found active on a-D-glucose. GOx has important applications in food and medical industries. Prototypes of miniature enzymatic fuel cells using GOx as catalyst have already been produced.1 GOx is found in many organisms such as insects and fungi2,3 and the extraction and characterization of GOx from different sources is an active research area. For biotechnological applications it is extracted mainly from the Aspergillus and Penicillium fungi. To improve its catalytic properties and structural stability is an important area of bioengineering.4–9 Although Aspergillus niger is currently the most common source as well as the best characterized, GOx from Penicillium amagasakiense (GOx-pa) has been shown to catalyze the glucose oxidation more efficiently. This makes it a better candidate to be used in bioelectrochemical nanodevices.10
C 2014 WILEY PERIODICALS, INC. V
The GOx-pa is a homodimer formed by two identical chains (see Fig. 1). It belongs to the a/b-fold class according to the SCOP11 classification. The monomers of GOx-pa contain 587 amino acids each, and one prosthetic group, flavin adenine dinucleotide (FAD). In the crystallized protein (PDB)12 ID: 1GPE13 each monomer has a molecular weight of 66.7 kDa with four N-glycosylation sites at positions ASN93, ASN165, ASN357, and ASN392, respectively. Three different sugars are bound at these positions: N-acetyl-D-glucosamine, b-D-mannose, and a-D-mannose. Their binding order is shown in Supporting Information (Fig. S1).
Additional Supporting Information may be found in the online version of this article. *Correspondence to: Francesca Mocci, Dipartimento di Scienze Chimiche e Geologiche, Cagliari University, Cagliari, Italy. E-mail: [email protected] Received 12 December 2013; Revised 22 March 2014; Accepted 29 April 2014 Published online 9 May 2014 in Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/prot.24596
PROTEINS
2353
G. Todde et al.
Figure 1 Cartoon representation of GOx. Colors legend: monomers (ice-blue and blue), FAD (green), sugars (pink), and dimer interface (red). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Unfortunately both GOx-an and GOx-pa suffer from relatively short operational lifetimes when used in nanodevices due to their relatively fragile structure. To understand how the inactivation of this protein occurs, several experimental studies have focused on the thermally and/ or chemically induced denaturation of GOx in aqueous solution.14–17 However, results obtained in these experimental studies do not allow to identify at the residue level where the denaturation starts to take place. This would be of great help to propose protein modifications capable of strengthening the protein. Atomistic molecular dynamics (MD) simulations have been used in the past to study unfolding and to identify the possible pathways in denaturation processes.18,19 In computer simulations the denaturation is conveniently induced and speeded up by increasing the temperature, or using solvents including agents known to destabilize the protein structure like urea or guanadinium.18,20–24 The pressure influence on protein denaturation process has also been studied.25–27 Moreover, artificial forces acting on selected groups of atoms have been used to trigger a denaturation.28,29 Concerning the relationship between high-temperature unfolding and unfolding under native conditions, one very important observation from the simulations of thermal denaturation is that “increasing the temperature accelerates protein unfolding without changing the pathway”.21 This means that it can be assumed that the thermal unfolding follows an Arrhenius law behavior, and the increase in temperature can efficiently compensate the limited timescale accessible to MD simulations.
2354
PROTEINS
MD studies of protein denaturation have been comprehensively reviewed recently by Toofanny and Daggett.30 Furthermore previous denaturation studies31 have shown that in order to adequately sample the denaturation process, an order of 5–10 simulations is required. Understanding of how the breakdown process starts and proceeds on molecular level can help to artificially improve the stability of proteins.18,31–34 Information obtained from detailed denaturation studies have already led to a design of a more robust enzyme, by strengthening its weak regions and transition points. In fact, mutating key residues close to structural transition points along the denaturation process have been proven to be capable to strengthen the whole structure of the enzyme and/or to make the moiety around the catalytic center more stable. Yoda et al.32 studied goat lactoalbumin by MD simulations and Circular Dichroism (CD) experiments and found that a single amino acid substitution was capable of increasing the stability of the native state by about 3.5 kcal/mol. The structural disruptions observed in the simulations were found to be consistent with the experimental results. Also the results obtained by Pikkemaat et al.33 on haloalkane dehalogenase and of Badieyan et al.34 on a glycoside hydrolases show that MD simulations are capable of identifying the protein domains that can be strengthened. In this study, we simulate the thermally induced denaturation process of GOx-pa at different temperatures from room temperature up to 600 K and follow the denaturation pathways in a systematic way. To the best of our knowledge, our work is the first of its kind performed for such large protein as GOx-pa. The results from our simulations would allow us to rationally propose mutations to obtain thermally a more stable protein. Here we focus on determination and characterization of the transition state (TS) of the denaturation process. The full denaturation process of GOx will be discussed in our upcoming study. Although already shortly discussed here, we will propose amino acid based modifications to strengthen the enzyme to be used in bionanodevices in a later communication. MATERIALS AND METHODS All the MD simulations were performed using Gromacs 4.5.335–38 and the amber99sb39 force field. Parameters for the cofactor FAD that were not present in the Amber force field were taken, except for the charges, from the work of Schneider and Suhnel.40 The partial atomic charges were calculated following the method used in amber99sb, using the restrained electrostatic potential procedure41 at the Hartree–Fock theory level with 6–31G* basis set. During the geometry optimization, the cofactor was kept in an elongated conformation, close to that observed in the crystal structure, by freezing the relevant dihedral angles. The starting
MD Study of GOx Denaturation
structure used for all simulations was the crystal structure solved by X-ray diffraction at 1.80 A˚ resolution obtained from the Protein Data Bank12 (PDB ID: 1GPE13). The simulations of the dimer were performed at four different temperatures: 300, 400, 500, and 600 K. Although the commonly used force field parameters are expected to be valid at ambient temperatures, numerous studies show that the same parameters can be applied also when working at very much higher temperatures.18–21,30,34 We have performed three independent simulations at both 500 and 600 K starting from different structures, randomly chosen from the first 10 ns produced at 300 K. Since GOx is a homodimer, every simulation of the dimer, in fact consists of two simulations of the monomer. Hence, totally 14 simulations of the monomer were produced at high temperature. All simulations followed the same protocol. The equilibration procedure, described in the following, is similar to that used by Daggett and Beck.31 The potential energy of the starting structure was briefly minimized (500 steps) with a mixed conjugate-gradient and steepestdescent algorithm. The resulting structure was placed in an empty cubic box subsequently filled with 68,000 TIP3P42 water molecules. To make the system neutral 32 water molecules were replaced by 32 sodium ions. The system was then energy-minimized (water only) before a short MD simulation (10 ps) of the solvent was performed keeping the solute fixed. Then the protein was minimized keeping the water constrained and a final minimization of the entire system was performed. At this point a 20-ns MD simulation at 300 K in NPT ensemble was carried out. The system was heated linearly (20 K per ps) in all simulations in the NVT ensemble. The lengths of the production runs were 20 ns at 300 K, 500 ns at 400 K, 150 ns at 500 K, and 15 ns at 600 K. All production runs were carried out in the NVT ensemble (except at 300 K where NPT was used) using 1 fs as timestep and the trajectories were saved every 1 ps (10 ps at 400 K). The Nose-Hoover thermostat43,44 was applied in all simulations and the Parrinello–Rahman barostat45,46 in the 300 K simulation. In all simulations the van der Waals interactions were computed within a cutoff distance of 10 A˚ and the Coulombic interactions were treated with the PME47,48 method. Bonds involving hydrogen atoms were treated by LINCS algorithm.49 All simulated systems are composed of 225,000 particles. All the produced trajectories sum up for a total simulated time of 1015 ls. Converted to a single AMD Opteron 6172 processor, the total computational time used would have been about 580,000 CPU hours. Analysis of the structural modifications involved the calculation of contacts and hydrogen bonds (HBs) between residues. A contact is defined in the following manner: each heavy atom (C, N, O, and S) is considered on average equal to a volume of 19.5 A˚3,50 so the vol-
Figure 2 Ca root mean square fluctuation. (A) and (B) show RMSF at 300 K (violet), respectively, for the first and second monomer. The crystallographic B-factors13 are converted to RMS fluctuation by the equation RMSF 5 (3B/8p2)1/2 and are shown by the orange line. The segments on the top of each plot represent the position of b-strands (black) and a-helices (red), while the vertical blue lines represent the position of the glycosylation sites. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
ume of each amino acid is computed. An average volume of 162.8 A˚3 is obtained and this volume is assimilated to a sphere with a radius of 3.4 A˚. Hence, two nonconsecutive residues are considered in contact when the distance between their centers of mass is shorter than 7 A˚ (twice the mean sphere radius in order to avoid false counts from the second shell of neighbors). A HB is formed when the distance between donor (NH, OH) and acceptor (N,O) is shorter than 3.4 A˚ and the angle Acceptor-Donor-Hydrogen is smaller than 60 . The secondary structure (SS) is assigned according to the Dictionary of Protein Secondary Structure program.51,52
RESULTS AND DISCUSSION Reference simulation at 300 K
The simulation performed at 300 K was used as a reference to test the stability of the force field and of the adopted simulation protocol. During the 20 ns simulated at 300 K the protein was found to be stable. The Ca RMSD (see Fig. S2 in Supporting Information) of the single monomers was
