**3.1 HRP model system**

The direct applications for protein structures produced by molecular modeling techniques such as homology modeling, include identification of structural and functional regions within a protein, and can be exploted as a lead for further experimental studies such as mutation analysis, catalytic and electrochemical characterization [27]. Furthermore, if homology modeling is combined with other computational methods such as molecular dynamics or quantum mechanics, the produced model can be used to screen different applications.

In our case, the available information comes from the genetic sequence of HRP enzyme and incomplete structures deposited on PDB, hence, it is necessary to reconstruct and evaluate our final model system based on homology models.

Our initial model was reconstructed with i-tasser metaserver, using as template the resting state horseradish peroxidase (1H5A), and building a waterbox for further molecular dynamics analysis (**Figure 4**).

#### **Figure 4.**

*Horseradish peroxidase homology model reconstruction and MD system. (A) Modeling pipeline applied to reconstruction of HRP, (B) final solvated HRP model for MD simulations.*

### **3.2 Conformational sampling**

The homology model system of HRP was subjected to an initial conformational search with GaMD method for 500 ns of simulation which is comparable to conventional MD simulations on the order of microseconds. This data allow us to reconstruct an energy hypersurface from three structural and energetic variables: total energy of protein, RMSD of backbone and superficial hydrogen bonds count. The total energy of the proteins indicates how close is our model system to a minimum energy state; RMSD of the alpha carbons is related to changes in relative positions of the atoms, regarding to the X-ray crystal model (**Figure 5**); and the hydrogen bonds count between aspartic acid (D), glutamic acid (E), lysine (K) and arginine (R) residues, indicate how prone the enzyme is to bond to AuNCs linkers functional groups (NH3 + or COOH). With this approach, energetically more favorable conformations were selected under two criteria: (1) Low RMSD values, *i.e.*, preserved HRP structure necessary for catalysis; (2) High solvent exposure of groups necessary for the esterification reaction, i.e., selecting those conformations where the probability of having the previous residues exposed to the solvent is higher. At the end of this sampling, three HRP conformations with the necessary intermolecular interactions were selected; with this approach we try to theoretically predict the formation of the amide bond in the solvent exposed residues, while the over-all enzyme structure remains correctly folded (**Figure 5**).

Using the GaMD approach, we were able to explore the conformational diversity along different reaction coordinates of the HRP enzyme, and to predict more efficient couplings to AuNCs. This approximation over the boost potentials allows a reconstruction with less noise on the potential energy, which results in a more robust method to deal with artifacts during variations in the potential energy [57, 58]. Hence, our data suggest that during simulation, there are different enzyme conformations which preferentially bind to AuNCs and increase the interactions with the linker functional groups.

#### **Figure 5.**

*Conformational sampling of HRP. The hypersurface was reconstructed using the second order McLauren cumulative expansion method. The free energy values were extracted from each conformation, depending on its structural variation and the solvent-accessible hydrogen bonds of amino acids D, E, K, and R.*

**109**

**Figure 6.**

*Design of Bioelectrochemical Interfaces Assisted by Molecular Dynamics Simulations*

tion and simultaneously preserve structure and reactivity of the enzyme.

particles modified with 4-aminothiophenol (**Figure 6A**) [24].

The main objective of this work is to give insights using molecular modeling methodologies, for the design of more efficient BEI, based on HRP enzyme and AuNCs. The structural description of HRP with molecular dynamics methodologies can help to design experimental strategies to increase the success of the immobiliza-

In order to estimate the immobilization efficiency between GaMD selected structures of HRP and AuCNs functionalized with linkers, free energy calculations were performed with adaptive steered molecular dynamics. The prediction model system is intended to recreate the interface previously described in Section 1.2 (**Figure 1**). GC was used as current collector and HRP was coupled to functionalized-gold nano

For all the calculations, the average electrostatic surface potential was directly extracted from GaMD, since electrostatic term is computed for each integration step. HRP showed two well characterized electrostatic isosurfaces, with a prominent

*Gold nanoparticles experimental and theoretical model systems. (A) SEM image of gold nanoparticle electrodeposit (average particle diameter ≈ 18 nm) in glassy carbon electrode (GC), (B) AuNCs modified with* 

*4-aminothiophenol linker used for free energy calculations.*

*DOI: http://dx.doi.org/10.5772/intechopen.93884*

**3.3 Interaction between modified 6 AuNCs and HRP**

*Design of Bioelectrochemical Interfaces Assisted by Molecular Dynamics Simulations DOI: http://dx.doi.org/10.5772/intechopen.93884*

## **3.3 Interaction between modified 6 AuNCs and HRP**

*Homology Molecular Modeling - Perspectives and Applications*

The homology model system of HRP was subjected to an initial conformational

or COOH). With this approach, energetically

search with GaMD method for 500 ns of simulation which is comparable to conventional MD simulations on the order of microseconds. This data allow us to reconstruct an energy hypersurface from three structural and energetic variables: total energy of protein, RMSD of backbone and superficial hydrogen bonds count. The total energy of the proteins indicates how close is our model system to a minimum energy state; RMSD of the alpha carbons is related to changes in relative positions of the atoms, regarding to the X-ray crystal model (**Figure 5**); and the hydrogen bonds count between aspartic acid (D), glutamic acid (E), lysine (K) and arginine (R) residues, indicate how prone the enzyme is to bond to AuNCs

more favorable conformations were selected under two criteria: (1) Low RMSD values, *i.e.*, preserved HRP structure necessary for catalysis; (2) High solvent exposure of groups necessary for the esterification reaction, i.e., selecting those conformations where the probability of having the previous residues exposed to the solvent is higher. At the end of this sampling, three HRP conformations with the necessary intermolecular interactions were selected; with this approach we try to theoretically predict the formation of the amide bond in the solvent exposed residues, while the over-all enzyme structure remains correctly folded (**Figure 5**). Using the GaMD approach, we were able to explore the conformational diversity along different reaction coordinates of the HRP enzyme, and to predict more efficient couplings to AuNCs. This approximation over the boost potentials allows a reconstruction with less noise on the potential energy, which results in a more robust method to deal with artifacts during variations in the potential energy [57, 58]. Hence, our data suggest that during simulation, there are different enzyme conformations which preferentially bind to AuNCs and increase the

*Conformational sampling of HRP. The hypersurface was reconstructed using the second order McLauren cumulative expansion method. The free energy values were extracted from each conformation, depending on its* 

*structural variation and the solvent-accessible hydrogen bonds of amino acids D, E, K, and R.*

+

**3.2 Conformational sampling**

linkers functional groups (NH3

interactions with the linker functional groups.

**108**

**Figure 5.**

The main objective of this work is to give insights using molecular modeling methodologies, for the design of more efficient BEI, based on HRP enzyme and AuNCs. The structural description of HRP with molecular dynamics methodologies can help to design experimental strategies to increase the success of the immobilization and simultaneously preserve structure and reactivity of the enzyme.

In order to estimate the immobilization efficiency between GaMD selected structures of HRP and AuCNs functionalized with linkers, free energy calculations were performed with adaptive steered molecular dynamics. The prediction model system is intended to recreate the interface previously described in Section 1.2 (**Figure 1**). GC was used as current collector and HRP was coupled to functionalized-gold nano particles modified with 4-aminothiophenol (**Figure 6A**) [24].

For all the calculations, the average electrostatic surface potential was directly extracted from GaMD, since electrostatic term is computed for each integration step. HRP showed two well characterized electrostatic isosurfaces, with a prominent

**Figure 6.**

*Gold nanoparticles experimental and theoretical model systems. (A) SEM image of gold nanoparticle electrodeposit (average particle diameter ≈ 18 nm) in glassy carbon electrode (GC), (B) AuNCs modified with 4-aminothiophenol linker used for free energy calculations.*

negatively charged surface derived from the presence of 12 exposed acidic residues. Each electrostatic potential was plotted as isosurface over each HRP conformation.

The free energy calculations between negatively charged isosurface of HRP and AuNCs modified with ATP linker (positively charged NH3 + ) showed coupling values of 60 Kcal · mol−1, with a minimum distance of 22 Å from the Fe(III) of the heme group (~ 2 Å to the AuNCs NH3 + functional group) (**Figure 7A**). On the other hand, the positively charged isosurface of HRP showed higher coupling energies of ~120 Kcal · mol−1 at a minimum distance of ~21.5 Å to the same AuNCs functional group. For the AuNCs modified with MBA, the coupling energy of the negatively charged isosurface of HRP was 130 Kcal · mol−1 at 33 Å distance, (~3 Å from surface functional groups), while the positive HRP surface showed free coupling energies of ~75 Kcal · mol−1 and minimum distances of 32.5 Å from Fe(III) heme group (3.5 Å from surface functional group) (**Figure 7B**).

The difference between both HRP isosurfaces with the positively charged surface imposed by ATP was ~70 Kcal · mol−1 which means that the enzyme electrostatic potential imposed a clear effect for the coupling to AuNCs. However, an overall difference of 15 Kcal · mol−1 between both linkers showed that the coupling assays were energetically more stable for the negative isosurfaces of the HRP enzyme than the positively charged.

**Figure 7.**

*Coupling free energy profiles of HRP isosurfaces and AuNCs. (A) Free energy profile of HRP coupled to 6 AuNCs-ATP, (B) molecular interaction model of negative isosurface of HRP and 6 AuNCs-ATP, (C) free energy profile of HRP coupled to 6 AuNCs-MBA, (D) positive isosurface of HRP and 6 AuNCs-MBA.*

**111**

**Acknowledgements**

**4. Conclusions**

**Figure 8.**

*Design of Bioelectrochemical Interfaces Assisted by Molecular Dynamics Simulations*

The free energy calculations suggested that HRP showed an energy minimum around the bond distance between the carboxyl groups and the amino groups [59, 60]. Our data suggest that the intermolecular interactions guided by negative electrostatic surfaces of HRP, are of lower free energy for the positive charged AuNCs-ATP, than the negative charged AuNCs-MBA (**Figure 8**). The high electrostatic attraction between HRP and AuNCs-ATP would promote an efficient electrochemical response. On the other hand, the aforementioned interactions between HRP and negative charged linkers, like MBA, are less stable and result in higher free energy

The findings of this study on HRP coupled to gold nanoclusters, indicate that the polarized electrostatic isosurface potentials are key factors to select the most efficient linker for coupling. The resulted interaction energy and distance between HRP and AuNCs-ATP are adequate to promote the formation of covalent bonds between acidic residues and amino functional groups. The evidence from this study points towards the idea that molecular simulation methods, such as homology modeling and molecular dynamics are valuable tools to take into account for design of BEI. Our previous multidisciplinary work of VP6 capsids has led us to conclude that molecular dynamics simulations elucidate structural determinants to understand the behavior of biomolecules on BEI. An implication for using homology models coupled to molecular dynamics, is the possibility of widely sample the conformational space of biomolecules probes before electrochemical experimentation. Further theoretical and experimental studies are necessary to describe the interaction with other functional group linkers and validate by electrochemical techniques the real effect of the charge difference between AuNCs-ATP and AuNCs-MBA on the redox response of this enzyme, respectively.

The authors acknowledge the financial support of CONACYT A1-S-15877, PAPIIT-DGAPA-UNAM-IN204218 and Centro Nacional de Supercomputo-IPICyT,

coupling profiles and lower electrostatic attraction for this interfaces.

*Different interactions displayed between HRP and AuNCs-linkers scaffolds.*

*DOI: http://dx.doi.org/10.5772/intechopen.93884*

*Design of Bioelectrochemical Interfaces Assisted by Molecular Dynamics Simulations DOI: http://dx.doi.org/10.5772/intechopen.93884*

**Figure 8.** *Different interactions displayed between HRP and AuNCs-linkers scaffolds.*

The free energy calculations suggested that HRP showed an energy minimum around the bond distance between the carboxyl groups and the amino groups [59, 60]. Our data suggest that the intermolecular interactions guided by negative electrostatic surfaces of HRP, are of lower free energy for the positive charged AuNCs-ATP, than the negative charged AuNCs-MBA (**Figure 8**). The high electrostatic attraction between HRP and AuNCs-ATP would promote an efficient electrochemical response. On the other hand, the aforementioned interactions between HRP and negative charged linkers, like MBA, are less stable and result in higher free energy coupling profiles and lower electrostatic attraction for this interfaces.
