The design of bioelectrochemical interfaces (BEI) is an interesting topic that recently demands attention. The synergy between biomolecules and chemical components is necessary to achieve high molecular selectivity and sensitivity for the development of biosensors, synthesis of different compounds, or catalytic processes. For most BEI, the charge transfer process occurs in environments with particular chemical conditions; modeling these environments is a challenging task and requires multidisciplinary efforts. These interfaces can be composed of biomolecules, such as proteins, DNA, or more complex systems like microorganisms. Oxidoreductases enzymes are good candidates, among others, due to their catalytic activities and structural characteristics. In BEI, enzymes are immobilized on conductive surfaces to improve charge transfer processes. Covalent immobilization is the most common method to prolong lifetime or modulate the detection process. However, it is necessary to implement new methodologies that allow the selection of the best candidates for a more efficient design. Homology modeling of oxidoreductases combined with Molecular Dynamics (MD) simulation methods are alternative and already routinely used tools to investigate the structure, dynamics, and thermodynamics of biological molecules. Our motivation is to show different techniques of molecular modeling (Homology Modeling, Gaussian accelerated molecular dynamics, directed adaptive molecular dynamics and electrostatic surface calculations), and using horseradish peroxidase as a model to understand the interactions between biomolecules and gold nanoclusters (as current collector). Additionally, we present our previous studies considering molecular simulations and we discuss recent advances in biomolecular simulations aimed at biosensor design.
- bioelectrochemical interfaces
- homology modeling
- covalent immobilization
- gold nanoclusters
- molecular dynamics
For some years now, the design and construction of bioelectrochemical interfaces (BEI), ranging from electrochemical biosensors (EC) for analytical applications, biofuel cells (BFC) to the development of biocomputing systems for information processing, have been topics where the scientific community has considerable participation. The implementation and integration of biomolecules with electronic elements or conductive surfaces to produce bioelectronic devices is a current topic and gathers great importance for developing biosensors that have an essential role in clinical applications, food quality control analysis and forensic medicine [1, 2, 3]. These applications always involve high sensitivity, low-sample volume, and low-cost production. In most BEI, charge transfer processes are carried out in environments with high ionic concentration, necessary for their function. These interfaces can be composed of macromolecules such as proteins, peptides, or more complex systems such as bacteria. Since the activity and lifetime depend on the correct interaction between the conductive surfaces (current collector) and the biomaterials, to achieve good bioelectrochemical responses, the enzymes need to improve the orientation towards the surface of the current collector [4, 5].
1.1 Covalent immobilization through alkanethiol linkers
The immobilization of various biomaterials is an important issue for BEI design, the covalent immobilization of enzymes stands out among the different enzymatic coupling strategies since the bonds formed through the linkers anchored to the current collector promotes direct charge transfer responses [6, 7]. Alcakanothiols are organic molecules widely used to establish self-assembled monolayers (SAMs) on the surface of gold electrodes [8, 9, 10, 11]. The thiol groups chemisorb on gold electrodes forming strong thiolate-gold bonds (Au-S). The resulting monolayer may be used to design molecular scaffolds to couple enzymes [8, 12]. The main drawback of covalent immobilization systems is the risk of having a low surface concentration of enzymes that are either active or correctly oriented for direct charge transfer; which may lead to a low-efficiency bioelectrochemical response . Therefore, having active enzymes and a good bioelectrochemical response becomes an important task. Aromatic molecules like 4-aminothiophenol, structurally mimics enzymatic substrates and can promote stable and direct contacts, necessary for efficient bioelectrochemical reactions.
Researches have proposed nanomaterials as support matrices for the design of BEI to increase the surface concentration of active enzymes on the current collector. During the last decade, nanomaterials, coupled with enzymes, have had significant relevance in the design of biosensors [5, 14, 15, 16]. Current advances in synthesis methodologies of these materials allow having a wide variety of nanomaterials with different sizes, shapes, surface charges, and physicochemical characteristics [4, 14, 17].
Some of these nanomaterials can be modified in different ways to improve biocompatibility. Biocompatibility refers to the biological nature events that do not interfere with those of electronic signal transduction and vice versa. Metallic nanoparticles are promising materials that increase the electroactive area and improve the sensitivity and stability of the attached enzymes on the electrodes, bringing the enzymatic active site close to the electrode (the redox cofactor should not exceed a 20 Å distance) to achieve direct charge transfer reactions [14, 18, 19, 20, 21, 22, 23].
In our group, we have established the conditions to improve the electrochemical response of horseradish peroxidase (HRP) coupled to gold nanoparticles modified with the aromatic alkanothiol 4-aminotiophenol, using glassy carbon (GC) as current collector . We proposed that a critical factor in HRP bioelectrochemical response is to promote conformational changes and proper enzyme orientation when coupled to the electrode. The protein conformation around the prosthetic group will determine the redox potential of the Fe3+/Fe2+ couple. The activity towards H2O2 also depends on the protein structural arrangements [25, 26]. We found that a proper environment for the enzyme activity was achieved by increasing the distance between gold nanoparticles. The voltammetric studies of the prosthetic heme group showed significant differences between the enzyme immobilized on randomly deposited gold nanoparticles and the enzyme immobilized on a well-dispersed gold nanoparticles deposit (Figure 1). The dispersed electrode improved the electrochemical response of the enzyme; this fact showed that the distance between enzymes is important, probably because longer distances decreased the steric impediments between the enzymes, and it appears to be a better immobilization strategy. With these results, we proposed that an essential part of BEI design is related to the structure of the biological systems. The electrochemical sensing is dependent on the recognition or transformation of the substrate, as previously observed .
The advent of high-resolution and robust techniques as X-ray crystallography [28, 29] and cryo-EM [30, 31, 32] has enormously contributed to the baggage of structural information available for bioelectrochemical applications. These techniques describe the atomic positions of enzymes, as well as the conformational dynamics resulted from recognition and binding processes, allowing the prediction of the possible electrochemical behavior of enzymes immobilized on electrodes.
1.2 Molecular modeling methodologies
The use of computational tools such as molecular mechanics and quantum mechanics for the study of chemical and biological reactions, involves a mathematical treatment of a large number of particles (hundreds of thousands of atoms that build up molecules like proteins, nucleic acids, lipids, and sugars) [33, 34]. In general terms, these methods are ideal for obtaining chemical and physical properties from three-dimensional molecular models [35, 36, 37]. It is necessary to correlate the structural models with their corresponding chemical, catalytic or biological properties; towards the rationalization and design of molecules for specific applications. Mainly, molecular Dynamics (MD) methods are based on the solution of Newton’s second law for all the atoms in the system, and help to predict the behavior of a biomolecule during a specific time. The integration of motion equations allows to analyze the trajectories corresponding to position, velocity, and acceleration of each particle of the system in any fraction of time . Unfortunately, the simplification of considering atoms as spheres, and bonds as springs or tensors (with different parameters depending on the type of bond), makes difficult to study bond-breaking and formation, or reaction mechanisms. Quantum mechanics (QM) mathematically describes the fundamental behavior of matter on an electronic scale by solving the Schrödinger equation for each atom in the system. QM also describes the behavior of atoms and molecules, in terms of their chemical reactivity, geometry, and their optical, electrical, magnetic, and mechanical properties [39, 40]. The coupled implementation of these techniques could describe the biomolecules involved in charge transfer processes, and design strategies for efficient BEI.
1.3 Molecular dynamics as a tool to predict the electrochemical activity in proteins
A challenging task in BEI design is the selection of robust biomolecules capable of tolerate non-physiological conditions, without loosing efficiency and conformation . The protein folding is governed by a series of molecular interactions between the amino acids and the surrounding chemical environment . Therefore, predicting their behavior with molecular models contributes to the optimization of resources before the development of BEI.
In our previous work the VP6 rotavirus capsid protein was encapsulated with the ionic polymer Nafion on GC. We demonstrated that this electrode could transfer electric charge applying an external potential, when using the redox probe potassium ferricyanide (K3[Fe(CN)6]). In parallel, we modeled by MD the electrostatic and conformational states of VP6 (Figure 2). This model suggested a route for ionic conductivity, where the electrostatic protein surface displayed a negative charge which interacted with the ferricyanide redox probe, promoting the charge transfer reaction. In order to show if this electrochemical activity was particular to proteins in general, under the same conditions as the experiments with VP6, bovine serum albumin (BSA), whose primary biological function is to transport different molecules through the bloodstream, was tested. However, we did not observe the charge transfer reaction of potassium ferricyanide when this protein was encapsulated (Figure 2). These results demonstrated that VP6 protein could be used as conductive scaffold for the development of different BEI applications .
Calculations of polarized electrostatic surface potentials from homology models of viral proteins, surprisingly showed that several capsids could display field effects as the recorded on VP6. In contrast, this effect cannot be displayed on non-charged proteins as BSA (Figure 3). Therefore, our previous data on viral capsids can serve as a workflow to identify candidate proteins or enzymes for construction of BEI devices. Altogether, theoretical and experimental reports on BEI [41, 42, 43] are examples of how structural information could help to elucidate the behavior of biological systems during electrochemical reactions. Nonetheless, the aforementioned studies depend on structural information from biomolecules, and certain biomolecules are very complex and highly labile, which difficult elucidation of atomic positions. However, since new protein sequences are continuously available, predictive methodologies as homology modeling becomes a valuable tool to asset structural information.
The aim of this multidisciplinary research is to perform molecular simulations on horseradish peroxidase as a model system, to generate data that can be used as selection criteria for design and further experimental validation in BEI development.
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) , and reconstructed using i-tasser metaserver . 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  and pH 7 was fixed accordingly to our previous studies.
Enzyme structural parameters were described using the ff14SB force field . 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 . 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 . 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, 52, 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 . Six gold nanoclusters (6 AuNCs) were placed as a thin layer, emulating our previous experimental system (~ 20 nm gold nanoparticles) .
2.3 Molecular coupling
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)  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 linkers, i.e., ATP amino groups and central Fe(III) atom of HRP heme (NH3+➔ Fe3+).
3. Results and discussion
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 . 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).
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,
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.
3.3 Interaction between modified 6 AuNCs and HRP
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) .
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 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.
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.
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, A.C. (National Supercomputing Center-CNS) grants for facilities at Thubat-Kaal 2 machine; very special thanks to Engineer Jesus Alaniz (CNS-IPICyT A.C) and Engineer Rogelio Elvira Morán (IER-UNAM) for their technical assistance (supercomputing and SEM analysis, respectively). W.I. García-García and G.A. Huerta-Miranda thank to CONCAYT for the PhD scholarship. A. Vidal-Limon thanks to CONACYT A1-S-15877 for the financial support of his academic fellowship.