In Silico Analysis of Glucose Oxidase H516r and H516d Mutations for an Enzymatic Fuel Cell

Glucose oxidase (GOx) is an oxido-reductase enzyme that catalyzes glucose into hydrogen peroxide and glucono delta-lactone (GDL). GOx has the potential to be used in the medical field. Numerous research concerning the usage of GOx to create enzymatic biofuel cells have been done, nevertheless the results obtained have not been optimal. This research aims to increase the Km values of GOx in order to increase its potential as a material for an enzymatic fuel cell. The amino acid histidine in position 516 is a residue in the active site that plays an important part in the process of glucose oxidation. In this research we mutated H516 by in silico twice resulting in the mutants R516 and D516. The mutations resulted in a change of the docking area for both mutants and in the docking affinity for H516D resulting in higher Km values. This shows that the H516 residue plays an important part in the functions of glucose oxidase and mutation into aspartate could improve glucose oxidase based enzymatic fuel cells.


INTRODUCTION
Advancements in the fields of nanoscience and nanotechnology have given a new direction on the design and the development on micro and nano sized tools and materials for the fields of industry, pharmacy, clinical, environmental and security. To maintain the usage of miniature electronic devices for a long period of time in a natural condition, a miniature energy source is needed. To overcome this problem a miniature biofuel cell can be used as an alternative energy source for microelectronic devices. Biofuel cells generally are divided into two types which are microbe-based and enzymebased. The primary advantage of enzyme based biofuel cells is its capacity to create a smaller sized microbial based biofuel cells but with the same amount of energy created (Ivnitski, Branch et al., 2006).
Enzymatic biofuel cells are a source of energy which have a high potential to be used on small electronic devices and biosensors such as the ones that can be used in bodies. The first fuel cell, which involved electrolysis of water, was discovered by Grove in 1839 (Song, Penmasta et al., 2011). In the last half century many innovations concerning enzymatic biofuel cells were achieved but it is still lacking in certain aspects notably its short lifetime, low voltage, and low efficiency because of only utilizing one enzyme for oxydation (Ivanov, Vidaković-Koch et al., 2010).
Glucose oxidase (GOx) is an oxidoreductase enzyme that catalyzes glucose into hydrogen peroxide and glucono delta-lactone (GDL). This enzyme is created by some types of fungus and insects. It also displays antibacterial activity in the presence of oxygen and glucose (Wong, Wong et al., 2008). GOx has the potential to be used in the medical field. It is commonly used with peroxidase reaction to visualize by colorimetry the formation of H 2 O 2 , this enables the determination of the amount of glucose in serum or blood plasma (Ambade, Sharma et al., 1998). A similar assay can also be done to observe the concentration of glucose in fermentation (Marco, Longo et al., 2016) and bioreactors (Goldrick, Lee et al., 2017). GOx based biosensors use electrodes replacing O¬2 to obtain electrons needed to oxidize glucose thus resulting in electrical currents proportional to the concentration of glucose (Blanford 2013).
Numerous research concerning the usage of GOx to create enzymatic biofuel cells have been done, this is due to its catalitical activity which can be done under physiological conditions. Nevertheless, the results obtained have not been optimal. The research of Chen and Yu (2009) concerning the usage of GOx as a material for a microfuel cell resulted in a Km value 1.28 mM, whereas the research of , Rohmayanti (2017), concerning the potential of GOx as a material for biosensors resulted in Km values of respectively 0.77 mM, 2.8 mM (Yu and Chen 2009, Ambarsari, Setyawati et al., 2016, Rohmayanti, Ambarsari et al., 2017. This research aims to increase the Km values of GOx in order to increase its potential as a material for an enzymatic fuel cell. The active site of the enzyme or the region of GOx where the reaction with the substrate takes place is buried in a pocket. It is defined by the residues Glu412, His516, and His559 (Petrović, Frank et al., 2017). The residue intended to be mutated is the amino acid Histidine at codon position 516 as no research has been conducted before to analyse the effects of the mutation of this residue on the function of GOx.
The residue intended to be mutated is the amino acid Histidine at codon position 516 located at the active site of the protein. Protein mutation is a way to find out the contribution of an amino acid in affecting the structure and function of a protein (Yuen and Liu 2007). This mutation was done in silico and the effects of the mutation were analysed through molecular docking simulation of GOx as the receptor and the glucose β-D-glucose using the wild-type GOx IPBCC 08.610 isolated according to previous research obtained from the NCBI database with the accession code MH953586.1 (Rohmayanti, Ambarsari et al., 2017). The advantages we hope will be able to be gained from this research is that the information obtained from the molecular docking simulation can be used as a basis for an vitro experiment to prove the validity of the results in a real life environment.

Ligand Preparation
The ligand used is β-D-glucose. The 2D ligand file is downloaded through the website pubchem in *.sdf format and converted to *.pdb format using the software Open Babel. Polar hydrogen atoms are added using Discovery Studio Visualizer 3.5 Client and saved as Ligand.pdb.

Receptor Molecule Preparation
The receptor used in this research is a glucose oxidase IPBCC 08.610 structure that has been isolated (Kurniatin, Ambarsari et al., 2020). The sequence is obtained from the NCBI online database with the accession code MH593586.1. The receptor data obtained in *.fasta format was changed to a 3D structure in *.pdb format using the web-based software SWISS-MODEL (Waterhouse, Bertoni et al., 2018) and saved as Receptor.pdb. The ligand and water still attached to the protein were removed, then atom hydrogens are added to the protein using the software Discovery Studio Visualizer 3.5 Client (Pilot 2016). The receptor is then saved in .pdb format.
Receptor Mutation using Chimera 1.11.2 (Yang, Lasker et al., 2012) The receptor file in .pdb format is opened using Chimera 1.11.2 (Pettersen, Goddard et al., 2004). Mutation is done by highlighting the residues which are to be mutated, this can be done by opening the sequence window from the menu and finding the desired residue in that sequence. Before mutating the residue we need to show all the atoms in its proximity by using the zone settings to show all the atoms within 5 angstroms. We then proceed to mutate the residue through the structure editing menu and changing the rotamer into the desired one where we would be given the choice to choose the position and structure of the resulting mutation. The mutation with the highest probability which does not collide with other atoms were chosen. The resulting mutated receptor is then saved in .pdb format. The residue His-516 is mutated into arginine (H516R) and aspartate (H516D). The receptor is then minimized using YASARA (Krieger, Joo et al., 2009).

Analysis of Mutant Structure Changes (Reddy, Vijayasarathy et al., 2006)
Mutant receptor structure is analyzed using the Ramachandran plot in the website by submitting our structure to http://mordred.bioc.cam.ac.uk/~rapper/rampag e.php and comparing it to the original receptor analysed with the same method. Codon level structure is analyzed by changing the data format from *.pdb to *.fasta using OpenBabel (O'Boyle, Banck et al., 2011) followed by reverse translation using BioEdit (Hall 1999).

Reseptor.pdbqt and Ligan.pdbqt Preparation and Docking Parameters (Syahputra 2014)
H516R.pdb is opened using Chimera 1.11.2, then non-polar hydrogen atoms are removed. The resulting file is then opened using AutoDockTools 1.5.6. Gasteiger charges are then calculated. The receptor is chosen as a Macromolecule and saved as H516R.pdbqt. Ligand.pdb is opened using AutoDockTools 1.5.6. The number of torsions of ligands is determined and saved as Ligand.pdbqt. Docking parameters are done with Set Map Types by choosing Ligand.pdbqt and Mutatedreceptor.pdbqt. The coordinates of the test ligand and the natural ligand is done in positions 45. 341, 9.731, and 53.837 (x, y, z) with a docking size of 20,18, and 28 (x, y, z) (Fikrika, Ambarsari et al., 2016).

Simulation of Ligand docking with Enzyme Receptor (Syahputra 2014)
The docking simulation is done with the aid of the Command Prompt (CMD). The receptor, ligands and Autodock Vina software are placed in the same folder. The docking command is written in CMD. Each docking simulation is done with the help of CMD. The test ligand is docked with each receptor (wildtype or mutant) in the active sites of the reseptor (targeted docking). The active sites of the reseptor are Glu412, His516, and His559 . The docking procedure is repeated using the command num_modes 20 times. The results obtained are in log.txt form. The best affinity energy (ΔG) is obtained in the first mode.

(Laskowski and Swindells 2011, Syahputra 2014)
The file of the docking result is opened using Discovery Studio Visualizer 2016. The first mode with the best affinity energy is chosen and copied into the receptor file. The receptor file chosen contains the ligand structure in .pdb format. The file is then analysed using Ligplot+ 1.4.5. The result of the analysis which consists of a 2D drawing of the interaction between the ligand and the protein shows the hydrogen bonds, the hydrophobic interaction with the length/ distance of the bonds. Determination of the Michaelis constant (Km) from the affinity energy (ΔG) of the ligand is done by using the formula ΔG = -RT ln Km.

RESULTS AND DISCUSSION
Glucose Oxidase IPBCC 08.610 3D

Structure Modelisation
The receptor used in this research is a sequence isolated previously, which is a published Aspergillus niger strain glucose oxidase on the website NCBI with the accession code MH593586.1 (Kurniatin, Ambarsari et al., 2020). The 3D structure of the glucose oxidase IPBCC 08.610 sequence was obtained through the process of homology modelling using the software SWISS MODEL (Waterhouse, Bertoni et al., 2018) with the known protein structure of glucose oxidase 1GAL as a template. The resulting 3D structure consists of a 583 residue glucose oxidase enzyme with a flavin-adenine dinucleotide (FAD) ligand which consists of 39 residue ( Figure 1). The sequence of the glucose oxidase in the file starts at residue 25 to residue 605 thus residue numbering appearing in subsequent figures in this document are 22 numbers higher than what it is in reality. Homology modeling is the construction of an atomic model of a target protein based solely on the target's amino acid sequence and the experimentally determined structures of homologous proteins, referred to as templates (Wedemeyer MJ 2019). The existing glucose oxidase A. niger file from the online database. The Protein Data Bank 1GAL was used as the template for the homology modelling of the glucose oxidase IPBCC 08.610 sequence (Burley, Berman et al., 2019). The resulting 3D model has a structure highly similar to that of the known 1GAL structure so we can conclude that the 3D structure modelisation was a success (Figure 1).

Figure 2. Ramachandran Plot of wild-type receptor
In order to determine the stability of the 3D model of the receptor created, a Ramachandran plot analysis was conducted. A Ramachandran plot is divided into four quadrants which are quadrant I or the most favoured regions, quadrant II or the additionally allowed regions, quadrant III or the generously allowed regions, and the quadrant IV or the disallowed regions (Hollingsworth and Karplus 2010). The results of the analysis of the wild type receptor shows that 95.7% of the residues are situated in the favoured region, 4.3% are in the allowed region and 0% are in the outlier region. This result shows us that the 3D model of the wildtype receptor created is of high quality and stable ( Figure 2).

Residue Mutation and Structural Changes
Mutations are changes in the genetic sequence and are the primary cause of diversity among organisms, these changes can occur at many different levels with a wide range of consequences. The mutations in this research are done in silico by using the program Chimera 1.11.2 through the means of residue substitution. Chimera is an extensible program for interactive visualization and analysis of molecular structures and related data, including density maps, supramolecular assemblies, sequence alignments, docking results, trajectories, and conformational ensembles (Yang, Lasker et al., 2012). The success of the technique used can be seen from the structure of the protein with the Ramachandran plot and also from the codon sequence. There are several types of mutations, the smallest of which are called point mutation, in which a single base pair is changed into another base pair. There are two types of substitutions according to the properties of the amino acid, conservative and non-conservative (Studer, Dessailly et al., 2013). The properties of an amino acid can be seen from its structure and polarity.
Mutation of the enzyme glucose oxidase IPBCC 08.610 was done twice, both on the same residue which was H516. Mutation is done by substituting a residue into another residue using the software Chimera 1.11.2, the details of the mutation can be seen in Table 1. An amino acid residue is coded by 3 nucleotides. Residue H516 is positioned 1548-1550 in the nucleotide sequence. The mutant enzyme H516R is mutated through single point mutation from the histidine (H) coding codon CAT into the Arginine (R) coding codon CGT. The mutant enzyme H516D is also mutated using single point mutation from the histidine (H) coding codon CAT into the Aspartate (D) coding codon GAT. The minimized state of the 3D structures of H516R and H516D can be seen in Figure 3.  The stability of the structures are analysed using ramachandran plots. The results of the ramachandran plots show that for the natural enzyme glucose oxidase IPBCC 08.610 97.8% of amino acid residues are in the favourable region, 2.1% are in the allowed region and 0.2% are in the outlier region ( Figure 2). Results for the mutant enzyme H516R show 97.4% amino acid residues in the favourable region and 2.6% in the allowed region whereas for the mutant enzyme H516D 95.7% of amino acid residues are in the allowed region and 4.3% are in the favourable region (data not shown). A residue percentage in the favoured region close to 100% shows a high quality sterochemical model (Reddy, Vijayasarathy et al., 2006). This shows that the receptors resulting from the mutations are high quality and can be used for further analysis.
The changes in nucleotide sequences can be analysed using reverse-translation with the software BioEdit. BioEdit is one of the most commonly used programs in molecular biology studies, its main usage is as sequence alignment editor. The GOx IPBC.08.610 contains 583 amino acids thus 1749 nucleotides in its sequence. The changes done resulting from the mutations are all substitutions of base pairs, no insertions or deletions of nucleotide bases were observed. According to Clark et al., (2005) conservative substitutions usually do not cause a big impact towards the protein whereas non conservative substitutions result in a bigger change towards the protein (Clark, Tavaré et al., 2005).

Molecular Docking and Active Site Determination
The molecular docking done in this research is between the ligand β-d-glucose with the wild type enzyme glucose oxidase IPBCC 08.610 as a receptor and the mutant receptors H516R and H516D. Before the simulation is conducted, the receptor and ligand must first undergo preparation, then be ready for molecular docking procedures.
In order to determine the site where the molecular docking will be done, a grid-box must be set up. This can be done through different ways such as redocking the natural ligand or by creating one around the residues that define the binding site. Redocking for this receptor is not possible because it is a simulated 3D model we created through SWISS-MODEL which does not contain cocrystallized ligands. Therefore the grid-box was created by using the active site amino acids. The active site was the region of an enzyme where substrate molecules bind and undergo a chemical reaction. The active site consists of residues that form temporary bonds with the substrate (binding site) and residues that catalyse a reaction of that substrate (catalytic site) (Albert 2010). The validation of a grid-box is carried out by the repeated docking between reseptor and β-d-glucose using the command num_modes 20 times. The best affinity energy (ΔG) is obtained as the best model for a grid-box.
Glucose oxidase from the mold Aspergillus niger oxidizes β-D-glucose with a wide variety of oxidizing substrates. It is an oxido-reductase that catalyses the oxidation of glucose to hydrogen peroxide and D-gluconoδ-lactone. The active site of the enzyme contains, in addition to FAD, three amino acid side chains that are intimately involved in catalysis: His516 with a pK(a)=6.9, and Glu412 with pK(a)=3.4 which is hydrogen bonded to His559, with pK(a)>8 (Leskovac, Trivic et al., 2005). Thus the grid box is created around the amino acids His516, His559, and Glu412.
After the molecular docking procedure is conducted, the binding site of the wild type receptor is analyzed along with the mutant receptors and compared with the binding site from the glucose oxidase 1GAL with the ligand β-D-glucose obtained from a previous research data (Meyer, Wohlfahrt et al., 1998). This comparison is done to ensure the validity of the 3D model created in this research by analysing the similarity in the binding site through the Binding Site Similarity (%BSS) in Table 3 and verifying the presence of the amino acids known to be actively involved in the oxidization process of β-D-glucose.
The binding site from the wild-type GOx enzyme and its mutants with the β-Dglucose ligand compared with the binding site 1GAL with the ligand β-D-glucose shows slight differences through their %BSS values. Compared to 1GAL, the wild type receptor shows 6 matching residues out of 10 resulting in a BSS value of 60%. The mutations show an increase in %BSS values with this comparison too with the H516R mutation having a %BSS value of 80% and the H516D mutation with a %BSS value of 70%. This change in BSS is due to the difference in residues with hydrogen bonds and hydrophobic bonds present in the wild type receptor (IPBCC 08.610 GOx) and the mutant receptors (H516 R and H516D) as mentioned before. The difference in %BSS and the hydrogen, hydrophobic bonds show that the mutated residue impacted the binding site. The result of this molecular docking with the autodock vina software is in the form of a -log‖ file which shows the affinity energy (ΔG) and the root mean square deviation (RSMD). According to the results here, for the binding site analysis are the ones with the lowest binding affinity (ΔG) with an RSMD value of under 2 Å, the RSMD is used to show the change in position of the ligand binding with the receptor. The position of the ligand is fine if the RSMD value is lower than 2 Å, an RSMD value over 2 Å shows that the position of the molecular docking does not conform with that of the X-Ray conformation. The RSMD values obtained and used as the basis of active site analysis in this studies for the wildtype receptor, the H516R receptor and the H516D receptor are 1.491 Å, 1.259 Å, 1.538 Å, respectively. These values all fall under the parameter of being under 2 Å, thus the position of the ligand with the active sites are relativily close (Table 3). The binding site of molecular docking simulation were analysed using the software Ligplot+ 1.45 that provides us the data concerning the type of bonds, the length of the hydrogen bonds, the atoms on the ligand that are binded with the receptor, and the amino acid residues on the enzyme receptors that interact with the ligand. The results show that the receptor IPBCC 08.610 has hydrogen bonds with the amino acid residues Asn 514, Arg 512, His 516, His 559, Asn 107, Tyr 68 and hydrophobic bonds with the residue Tyr 515 (Figure 4).

Figure 5.
Visualisation of molecular docking simulation between glucose oxidase (1GAL) and β-D-glucose (Meyer, Wohlfahrt et al., 1998) To verify the validity of the 3D model created, the active site of the model and its residues are compared to ones of a published model. The model used as comparison for this research is the 1GAL glucose oxidase receptor docked with a β-d-glucose ligand in the research done by Meyer et. al (1998) (Meyer, Wohlfahrt et al., 1998). The receptor 1GAL had hydrogen bonds with the residues His 516, His 559, Tyr 68, Thr 110, Asn 514, Phe 414, Arg 512 and hydrophobic bonds with the residues Asp 424, Trp 426, Tyr 515 ( Figure 5). The mutant receptor H516R has hydrogen bonds with the residues Asn 514, Arg 516, Tyr 68, His 559 and hydrophobic bonds with the residues Gly 108, Phe 414, Thr110, Trp 426 ( Figure 6). The mutant receptor H516D has hydrogen bonds with the residues Arg 512, Tyr 68, Asn 514, Asp 516, Thr 110, Gly 108, His 559, and hydrophobic bonds with the residues Trp 448 ( Figure 6).
From the comparison done of the test receptors with the 1GAL data from previous research, we can see that the amino acids His 516 (or the mutated versions) and His 559 are all present bound to the receptors (Figures 4  and 6). This shows that our 3D model created contains the same active site as the one from a known X-ray crystallized model thus validating the model simulated. The other amino acid known to be involved in the catalysis reaction Glu 412 is not present in our analysis of the active site. This is due to it not interacting directly with the ligand, but hydrogen bonded with His 559.

Affinity Energy, Michaelis Constante and Impact on Enzymatic Fuel Cells
The stability of ligand interaction with the receptor can be analyzed using the affinity energy obtained. Chemical reactions are separated into two according to energy needs, endergonic and exergonic. The release of energy due to the reaction is mathematically written with a negative value. The lower the affinity energy obtained, the reaction will be more spontaneous and the conformation more stable (Aquino Neto and De Andrade 2013). The stability of interaction between ligand and receptor can be interpreted as the strength of the inhibitor to compete with other inhibitors or the natural substrate of that receptor.
Binding affinity is the strength of the binding interaction between a single biomolecule (e.g. protein or DNA) to its ligand/binding partner (e.g. drug or inhibitor). Binding affinity is typically measured and reported by the equilibrium dissociation constant (Kd), which is used to evaluate and rank order strengths of bimolecular interactions. The smaller the Kd value, the greater the binding affinity of the ligand for its target. The larger the Kd value, the more weakly the target molecule and ligand are attracted to and bind to one another (Gupta, Chaudhary et al., 2015). Binding affinity is influenced by non-covalent intermolecular interactions such as hydrogen bonding, electrostatic interactions, hydrophobic and Van der Waals forces between the two molecules. In addition, binding affinity between a ligand and its target molecule may be affected by the presence of other molecules. In a chemical reaction the Michaelis constant Km is numerically equal to the substrate concentration at which the reaction rate is half of Vmax (the maximum reaction rate achieved by the system). Since the Km value is equal to the Kd value in a reaction and the values of the reaction is not known in a molecular docking experiment, we assume that Km=Kd for this study An enzymatic fuel cell is a type of fuel cell based on enzymes as catalysts instead of expensive metals. Currently, enzymatic biofuel cells are inefficient, have a short lifespan and do not produce much power (Elouarzaki, Cheng et al., 2018), in order to overcome this problem a numerous amount of studies have been done. One of the types of enzymatic biofuel cells available are the ones using glucose oxidase as the anode. According to the research by Chen and Yu (2009) (Yu and Chen 2009) the maximum power density (Pmax) of a biofuel cell reaches a saturated value when [Glucose] ≥ 20 mM in his research, this explains why a biofuel cell with 100 mM glucose has no response to additional glucose in his research. Thus an increase in Km value can decrease this saturation value and enhance the fuel utilization.
The affinity energy between the ligand β-D-glucose and the wild-type receptor IPBCC 08.610 is -6.5 kcal/mol. The affinity energi (ΔG) between the mutant receptor H516R and the ligand β-D-glucose is the same as the wildtype receptor -6.5 kcal/mol. The mutant receptor H516D has an affinity energi (ΔG) with the ligand β-D-glucose of -6.2 kcal/mol, this means that the histidine to aspartate mutation has increased the affinity energy of the enzyme (Table 4). The affinity energy from the molecular docking can be used to calculate the Michaelis constant (Km). The Km value can be calculated by using the formula: ΔG = RT ln(Km/Co).
The value of R used is 1.986 kal/mol and the value of T used is 298.15L. Co is the concentration in thermodynamics standard which is 1 mol/L (Kastritis and Bonvin 2013). According to calculations, the Km value resulting from the three enzymes are 17.187 µM, 17.187 µM, and 28.517 µM, respectively. These results show us that the resulting Km value of the IPBCC 08.610 enzyme is equal to the Km value of the H516R enzyme thus this mutation does not result in a significant difference towards the reaction of glucose binding with the receptor. On the other hand, the mutation H516D results in an increase of Km value compared to the Km value of the original enzyme therefore showing that this mutation has impacted the reaction of glucose binding with the receptor and shows us that a H516D mutation histidine to aspartate can improve the quality of glucose oxidase based enzymatic fuel cells.

CONCLUSION
The 3D model created based on the glucose oxidase enzyme IPBCC.08.610 with the NCBI accession code MH593586.1 has the same active sites with the enzyme with the known configuration 1GAL. His 516 mutation into Arg 516 shows no significant differences whereas the mutation into Asp 516 shows a slight increase in Km values. The molecular docking simulation shows that the H516D mutation has the potential to increase the effectiveness of glucose oxidase as an enzymatic fuel cell material.