**2. General techniques: Horseradish modeling and system building**

The X-ray crystal structure from horseradish peroxidase (HRP) was collected from Protein Data Bank (www.rcsb.org/, PDB code 1H5A) [45], and reconstructed using i-tasser metaserver [46]. The best model was selected based on lowest RMSD and higher C-value, to achieve the best initial coordinates. To improve the model system, all missing atoms and sidechain positions were modeled against available plant peroxidases on the PDB. A penta-coordinated resting state heme group was crystallized as Fe3+, and axial coordination with a histidine groups was retained.

The Protein Preparation Wizard routine (Schrödinger Maestro v2019-4, New York, 2019), was applied as a preparation step to correctly assign missing hydrogen atoms and protonation states of ionizable residues. Finally, PROPKA2.0 sub-routine was applied to whole model system [47] and pH 7 was fixed accordingly to our previous studies.

Enzyme structural parameters were described using the ff14SB force field [48]. The full systems comprised ~60,000 atoms in a cubic box of 15 Å length including TIP3P explicit water and Na + ions to ensure overall charge neutrality. Non-bonded interactions were calculated within a 12 Å cutoff, and long-range electrostatics were treated using the Particle-Mesh Ewald method [49]. SHAKE algorithm was enabled to constrain all bonds involving hydrogen during simulations. Conventional MD protocol for each system replica comprised: 1) 5000 steepest descent minimization steps followed by 10,000 conjugate gradient minimization steps; 2) 500 ps of progressive NVT heating from 0 to 300 K; 3) 5 ns of NVT equilibration, and 4) 10 ns of NPT dynamics at 300 K and 1 bar.

#### **2.1 Gaussian accelerated molecular dynamics (GaMD)**

Conventional MD simulations were carried out to initially relax full system coordinates, using AMBER18 software package [50]. Then, the GaMD module implemented in AMBER v. 18 was applied to perform extra 50 ns of short classic MD simulation to collect the statistics for calculating GaMD acceleration parameters. Extra 50 ns of short equilibration was applied after adding the boost potential, and finally three independent 500 ns GaMD production simulations with randomized initial atomic velocities. This method applies Gaussian functions (boost potentials) on the total potential energy term, which disturb the potential energy surface allowing the exploration of different energy states. In addition, functions that directly disturb the dihedral angles of the amide bonds were applied, promoting conformational changes. All GaMD simulations were performed with a dual-boost level by setting the reference energy to the lower bound. One boost potential is applied to the dihedral energetic term and the other to the potential energy term. The average and SD of the system potential energies were calculated every 500,000 steps (1 ns) for all simulation systems. The upper limit of the boost potential SD, σ0, was set to 6.0 kcal·mol−1 for both the dihedral and the total potential energetic terms [51–53]. The system temperature was ~298 K and 1 atm pressure, with integration steps of 2 fs. Three different system coordinates were extracted from GaMD trajectories, each one representing a different energy and conformational state.

#### **2.2 Building of nanomaterial structures**

The selected nanomaterials were Au nanoclusters (AuNCs) modified with two types of linkers: 4-aminothiophenol (ATP) and mercaptobenzoic acid (MBA); each cluster had 2.1 nm diameter and a nucleus composed of 96 Au atoms (Au314 (SR) 96). The nanostructures were modeled as reported elsewhere [54]. Six gold

**107**

**Figure 4.**

*Design of Bioelectrochemical Interfaces Assisted by Molecular Dynamics Simulations*

ers, i.e., ATP amino groups and central Fe(III) atom of HRP heme (NH3

produced model can be used to screen different applications.

further molecular dynamics analysis (**Figure 4**).

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

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

*Horseradish peroxidase homology model reconstruction and MD system. (A) Modeling pipeline applied to* 

*reconstruction of HRP, (B) final solvated HRP model for MD simulations.*

nanoclusters (6 AuNCs) were placed as a thin layer, emulating our previous experi-

Adaptive Steered Molecular Dynamics (ASMD) simulations [55, 56] started from the earlier selected conformational states, derived from the GaMD analysis. Due to the chemical complexity of the enzyme surface, the simulations were driven by the electrostatic surface (PME) [49] for each conformational state and the 6 AuNCs layer. Hence, the scans were performed at speeds of 10 Å / ns, which comprised a 10-step profile of 1 Å between the Fe(III) of the active site and the ATP amino group or MBA group carboxylic of the 6 AuNCs. Three independent replicas were performed through gradual scans, involving measurements of the free energy of coupling between HRP and 6 AuNCs. The free energy was measured by shortening the distance between the central Fe(III) atom and the amino groups of the link-

+

➔ Fe3+).

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

**2.3 Molecular coupling**

**3. Results and discussion**

**3.1 HRP model system**

mental system (~ 20 nm gold nanoparticles) [24].

nanoclusters (6 AuNCs) were placed as a thin layer, emulating our previous experimental system (~ 20 nm gold nanoparticles) [24].
