Computational Modeling Methods for Understanding the Interaction of Lignin and Its Derivatives with Oxidoreductases as Biocatalysts

Abstract

This chapter will be presented as follow. First, a brief introduction to structure and characterization of lignin and its derivatives is presented, as well as their importance as chemical scaffolds for obtaining value-added products in chemical, food, pharmaceutical and agriculture industry. Second, an extensive review of different reports using computational modeling methods—like molecular dynamics simulations, quantum mechanics and hybrid calculation methods, among others—in the understanding of enzyme-substrate interaction and biocatalysis will be presented. Third, and as last part of chapter, some hand picked examples from literature will be chosen as successful cases where the interplay between experiment and computation has given as a result protein engineered oxidoreductases with improved catalytic capabilities.

1. Introduction

Lignin, one of the more abundant and recalcitrant substrates available in biosphere, is a complex aromatic heteropolymer formed by radical polymerization of guaiacyl, syringyl, and p-hydroxyphenyl units linked by β-aryl ether linkages, biphenyl bonds and heterocyclic linkages, among others [1, 2]. This natural substrate and its derivatives are of upmost interest in biotechnological and chemical industry because the high-value products that can be obtained from their chemical or biological transformation. The principal strategies for depolymerizing lignin to produce valuable chemicals are catalytic reduction, cracking or hydrolysis reactions, catalytic oxidation and enzymatic reactions [3, 4, 5, 6]. The latter approach has proven to be more efficient compared to existing physical and chemical pretreatments of biomass that involve the use of costly equipment and reagents, while biocatalytic processes may represent an alternative environmentally friendly approach with lower energy requirements and cost [7]. In this regard, the isolation of lignin-degrading microbial strains may lead to the discovery of novel biocatalysts, like peroxidases and laccases among other enzymes that could be potentially useful for lignin valorization [8, 9]. Obtaining compounds derived from biotransformation processes, with interesting and new applications in chemical, cosmetic, pharmaceutical and food industries, has been the primary interest of several research groups across the chemistry and biology research community. Microorganisms, or their enzymatic systems, catalyze biotransformation processes and these are usually carried out by growing cultures, purified enzymes or immobilized cells, among others [10]. Those biological chemical reactions can be used in order to introduce chiral centers, to resolve racemates, to convert a particular functional group among several groups with similar reactivities, to functionalize a nonactivated carbon regioselectively and to convert labile molecules due to the mild reaction conditions under which they take place [10]. A very good example are the ligninolytic enzymes from white-rot fungus Phanerochaete chrysosporium, that have been used in several biotransformation and bioremediation processes like the oxidation, by the extracellular enzyme lignin peroxidase (LiP) of diphenyl compounds biphenyl, biphenylene, dibenzofuran, dibenzo-p-dioxin and diphenyl ether, which constitute an important class of environmentally persistent pollutants [11]. More recently, the transformation of three anti-inflammatory drugs such as diclofenac, ibuprofen and naproxen were carried out by pellets of P. chrysosporium in fed-batch bioreactors operating under continuous air supply or periodic pulsation of oxygen [12]. The authors found that, under operating conditions of the reactor, high elimination percentages for the pharmaceutical compounds (80–99%) were achieved. On the other hand, high-value products derived from the natural substrate of ligninolytic enzymes, lignin, have also been reported. First-generation fine chemicals can be obtained from lignosulfonate by oxidation, vanillin being the most important fine chemical directly produced from the hydrolysis-oxidation [13]. Vanillin is the main component of the natural vanilla extract. It is used as a flavoring agent in foods, beverages, pharmaceuticals and in the fragrance industry [14]. Besides its value as flavoring agent, vanillin is the starting reagent in the synthesis of several second-generation and important fine chemicals. It can be used in synthetic processes leading to several pharmaceutical chemicals like cyclovalone (digestant or choleretic) [15], etamivan (analeptic and a central nervous system and respiratory system stimulant) [16] and levodopa (an antiparkinson agent) [17]. According to the authors, even though only a very few of the first-generation fine chemicals from oxidation of lignosulfonates have been isolated, a huge variety of organic processes to pharmaceutical chemicals are based on them. Future discovery and optimization of more selective depolymerization processes of lignosulfonates may extend the employment of such fine chemicals as building blocks [13]. However, the chemical processes for obtaining vanillin (and other important and related chemical building blocks) from lignosulfonates have been very questioned due to the huge amounts of waste effluents produced that, combined with the growing public awareness on environmental issues, are leading to unsustainable effluent-treatment costs [18, 19, 20]. A clean and friendly environmental alternative, for the synthesis of such kind of valuable chemical building blocks, is the biotransformation catalyzed by ligninolytic enzymes from P. chrysosporium or other microorganisms. Of special interest are oxidoreductases that take advantage from the incorporation of different cofactors—such as heme, flavin and metal ions—to catalyze redox reactions. As recently reported, classical fungal oxidoreductases comprise basidiomycete ligninolytic peroxidases, and ascomycete and basidiomycete multicopper oxidases (MCO, mainly laccases) with different redox potentials and abilities to act on lignin-derived products [21]. Among those, ligninolytic peroxidases have been known for some 40 years, and representatives of the three main types—lignin peroxidase (LiP, EC 1.11.1.14), manganese peroxidase (MnP, EC 1.11.1.13) and versatile peroxidase (VP, EC 1.11.1.16)—have been extensively studied due to their biotechnological potential for the chemical modification and degradation of lignin and other recalcitrant compounds [22]. The molecular biology and structure-function of these lignin-degrading heme peroxidases have been already carefully described [23]. Other interesting oxidoreductases are the so-called dye-decolorizing peroxidases (DyPs, EC 1.11.1.19) that have been recently described and structurally characterized in basidiomycetes [24]. Its structural convergence with ligninolytic peroxidases, besides a long-range electron transfer (LRET) mechanism for oxidation of bulky lignin-derived and dye substrates, make them worth for study as another example of evolutionary convergence between phylogenetically unrelated enzymes oxidizing recalcitrant structures including lignin [21]. Finally, different copper-containing oxidoreductases have been related to lignocellulose degradation, being one of the most studied the laccases (EC 1.10.3.2) that, together with peroxidases, are the most thoroughly studied oxidoreductases in wood-rot fungi and, by far, the largest number of biotechnological applications have been reported for [25, 26]. The present chapter summarizes the most relevant experimental/computational evidence that has been published recently about these key oxidoreductases as biocatalysts, that are helping in the biotransformation and valorization of lignin, and its derivatives, as substrates into high-value products. Moreover, it will be emphasized the use of computational modeling methods in the understanding of substrate-enzyme interactions and their helping role in protein engineering as a key tool for tune the efficiency and specificity of this important kind of enzymes.

2. Computational modeling methods in enzyme-substrate studies

Different computational methodologies could be used in order to study and to understand the underlying mechanisms that make biocatalysts so efficiently degrading pollutants or biotransforming substrates into high-value products. Therefore different computational methods currently available could be of key importance in understanding those sophisticated processes taking place at the molecular level. All of them differ both in computational cost as well as in accuracy. For instance, methods like molecular docking attempt to predict the structure (or structures) of the intermolecular complex formed between two or more molecules. Docking is widely used to suggest the binding modes of protein inhibitors. Most of the docking algorithms are able to generate a large number of possible structures, so they also require a means to score each structure to identify those of most interest. The ‘docking problem’ is thus concerned with the generation and evaluation of plausible structures of intermolecular complexes [27, 28]. The main disadvantage of docking methods is their scoring function, because it may become insufficient to adequately represent all binding forces participating in the protein-substrate binding. However those methods offer us an important structural starting point for the application of more accurate methods that can take into account the dynamical behavior, energetics profile and electronic characteristics of the system under study. Once one has the structural starting point for protein-substrate obtained from molecular docking, the next computational method that can be applied is molecular dynamics (MD) simulations. This approach is very powerful in the description of very large systems like DNA-drug complexes [29, 30], polymer structures [31, 32], membrane proteins [33, 34], protein-ligand complexes [35, 36], etc. The simplicity of the potential energy function of this method makes to run very long simulations properly, in the nano and micro second time scale, and to extract valuable structural information and energy data from trajectories obtained. Other more refined molecular dynamics approaches like MD-FEP allow to calculate absolute enzyme-substrate binding free energies [37]. As a matter of fact, in all these methods it is mandatory to perform a parameterization of the ligand prior to run the molecular dynamics simulations. This is a severe step back in the simulation process and the successful and accurate representation of enzyme-substrate binding process, because the quality of the simulation is strongly dependent on the quality of the ligand parameterization. This problem is totally alleviated using quantum mechanics (QM) or quantum mechanics/molecular mechanics (QM/MM) methods, which lie at the other end of the computational methods spectrum as the more computational expensive approaches like ab initio Car-Parrinello molecular dynamics, QM and derived electronic structure methods in general. These methods are very accurate in the prediction of the reactivity, structural and thermodynamic properties of small molecular systems, but their application to the understanding of enzyme-substrate interaction and to the enzyme catalysis is very scarce yet. This can be explained in part by the high computational cost they demand and the limitations to include explicit solvation effects and conformational searching in a fast and cost affective way. On the other hand, in the last two decades, the hybrid QM/MM calculations have became the more popular methods within computational chemistry community. The broad range of applications goes from calculation of spectroscopic and excited-states properties, studies of enzymatic reactions, surface and chemical catalysis and, in a minor extend, protein-ligand interactions (drug design). The successfully application of the QM/MM method in those areas relies on the combination of accurate QM methods and fast MM computational algorithms. For a complete review of the method and a survey of biological applications, see Refs. [38, 39, 40] and references therein. In the next sections, the role and support of/from those computational methodologies in the experimental study and understanding of main oxidoreductases as biocatalysts degrading lignin and its derivatives will be reviewed.

2.1. Laccase

As pointed out before, laccases are the most studied member of oxidoreductases because their broad ranges of chemical and biotechnological applications [41, 42]. From the computational point of view, there is also a growing amount of research that have attempted to understand their catalytic mechanism and also to tuned their substrate affinity with biotechnological, environmental or chemical application purposes. Regarding laccase catalytic mechanism, it has been studied for decades but is still not completely elucidated, especially in terms of the reduction of dioxygen to water. Its key structural features have been also under investigation by several groups using techniques such as X-ray diffraction, electron paramagnetic resonance (EPR) spectroscopy and site-directed mutagenesis. For a more complete view of this topic see Ref. [43]. In 2013, the available literature about theoretical studies on the active-site structure, spectroscopic and thermodynamic properties, and reaction mechanism of multicopper oxidases (MCO) were collected by Rulíšek et al. [44]. They emphasized that complicated electronic structure of the trinuclear copper cluster (TNC) make it hard to obtain quantitatively correct spectroscopic data from the density functional theory (DFT) methods. Therefore, multi-reference wave function methods (e.g. CASPT2 and MRCI-S) would be necessary to obtain quantitative splitting of the lowest doublet and quartet states, reliable excitation energies of ligand-field states, and qualitatively correct predictions of the g tensors for these spin-frustrated systems. More recently, this same group reported a work that combined quantum mechanical and molecular mechanical free-energy perturbation (QM/MM-FEP) methods in combination with explicit solvent simulations to study the reaction mechanism of the multicopper oxidases. They studied redox potentials (RP), acidity constants, isomerization reactions, as well as water and O₂ binding reactions. A full reaction mechanism of the multicopper oxidases with atomic detail was proposed. They also showed that the two copper sites in the protein communicate so that redox potentials and acidity constants of one site are affected by up to 0.2 V or 3 pK_a units by a change in the oxidation state of the other site [45]. On the other hand, Hong et al. [46] studied the redox potentials and reorganization energies of the type 1 (T1) Cu site in four multicopper oxidases namely laccase from Trametes versicolor, a laccase-like enzyme isolated from Bacillus subtilis, CueO required for copper homeostasis in Escherichia coli and the small laccase (SLAC) from Streptomyces coelicolor. They used first principles DFT and QM/MM simulations to calculate redox potentials and reorganization energies that were in good agreement with experiments. Moreover, they proposed a mutated MCO with improved performance taking into account that an estimate of the tunability of the T1 Cu potential, and as a relatively high cathodic redox potential can be achieved with Cu in a well-defined protein environment. In a similar way, geometrical distortions (distances and angles) in the copper coordination sphere were introduced in a model of the three-coordinated T1 Cu site of laccase from T. versicolor in order to evaluate the effect on redox potential. All calculations were done within a DFT framework, where the best approach to increase redox potential was by distortion of the dihedral angle ω (defined as C_{methylthiolate} –S–Cu–N_ImA), rationalized as a decrease in the overlap of imidazole orbitals in the redox-active molecular orbital (β-LUMO) [47]. More recently, the pH dependence and mutants on the laccase redox potentials at T1 site were studied with QM/MM approaches. First, the protein optimized geometries were obtained with RI-BP86/def2-SVP/def2-TZVP(Cu):CHARMM QM/MM model and then the RPs were obtained at the M06/6-311++G**/SDD(Cu) level of theory. The authors found that the oxidation state of the TNC affected the T1 site RP by about 0.2–0.3 V, depending on the protein protonation state. Moreover, they also calculated the RP of a F463M laccase mutant protein to probe if their adopted methodology could describe the change upon mutation of this residue. Notably, they were able to reproduce the experimental change in the RP of a related enzyme (0.09–0.1 V) very closely (0.08 V) [48]. Other relevant research works have been published with the aim to understand the differential affinity of several substrates against laccase T1 site and to explain the electron transfer pathways and the enzyme residues involved in that process. For instance, Bello et al. [49] performed MD simulations of three substrates (2,6-dimethoxyphenol (DMP), syringaldazine (SGZ) and 2,2-azino-bis(3-ethylbenzothiazoline-6-sulfonic acid) (ABTS)) bounded to T1 site of MCOs from Thermus thermophilus, B. subtilis and E. coli. The authors found that binding modes of the electron-donor molecules to the electron transfer binding site were primarily attributed to hydrophobic contacts that are supposed to play an important role in substrate specificity. Moreover, some MCO-substrate complexes showed an electron-donor molecule conformation in which an electron could be directly transferred to the histidines coordinating T1 Cu, although for others additional electron transference pathways could be also feasible through the participation of charged residues during electron transfer. More recently, five lignin-related compounds, namely 2,6-dimethoxyphenol, ferulic acid, guaiacol, sinapic acid and vanillyl alcohol, were selected to demonstrate the key binding mechanisms between T. versicolor laccase and lignin. The authors performed molecular docking and molecular dynamics simulations to study and compare the interaction between various lignin model compounds and laccase with the aim to proposed the molecular basis for their binding, which might be helpful in the understanding of lignin-degrading reaction [50]. Based on molecular docking binding energies, of enzyme-substrate complexes, and their stability through 10 ns MD simulations, the authors concluded that hydrophobic interactions seemed to be necessary to the interaction of lignin/lignin model compounds with laccase, while hydrogen bonds (H-bonds) were alternative. In a similar research work, Awasthi et al. [51] also performed molecular docking, MD simulations and, additionally, molecular mechanics-Poisson Boltzmann surface area (MM-PBSA) analyses of several lignin model compounds namely sinapyl alcohol (monomer), guaiacyl 4-O-5 guaiacyl (dimer), syringyl β-O-4 syringyl β-O-4 sinapyl alcohol (trimer) and guaiacyl β-O-4 syringyl β-β syringyl β-O-4 guaiacyl (tetramer) bounded to plant and fungal laccases from Populus trichocarpa and T. versicolor, respectively. The authors concluded, from structural changes at the C-terminal region, differences in the interacting residues of the binding site, different binding modes of tetramer at the active site, differences in the root mean square deviation (RMSD), interaction energies, hydrogen-bond interactions and binding free-energy analyses during MD simulations; that their computational results supported a situation favorable for degradation of lignin polymer by fungal laccase meanwhile their synthesis is by plant laccase. The inhibition of laccases binding site has also been studied. For instance, the inhibition of laccase enzymatic catalytic activity by formetanate hydrochloride (FMT), a carbamate pesticide, was investigated by cyclic voltammetry and quantum-chemical calculations based on DFT with a protein fragmentation approach. The authors aimed to develop enzymatic electrochemical biosensors for inexpensive analytical detection of pesticides [52]. In other similar work, the authors performed molecular docking calculations to understand the inhibition mechanism of medicarpin on the lignin degradation enzyme laccase from T. versicolor. Those preliminary results would be useful in developing wood-preserving formulations based on medicarpin, according to authors [53]. Finally, Christensen et al. applied multiple-seed MD simulations to a T. versicolor laccase in response to variable ionic strengths, temperatures and glycosylation status. The persistence of backbone hydrogen bonds was identified as a key descriptor of structural response to environment, whereas solvent-accessibility, radius of gyration and fluctuations were only locally relevant [54]. When moving to study of laccases from rational protein design through directed evolution and computational approaches, the work done by the groups of Guallar and Camarero is a reference point. For instance, in 2015 Monza et al. [55] reported a novel computational approach that combined efficient conformational sampling and quick reactivity scoring in order to understand how substrate oxidation was improved during a directed evolution experiment of a fungal laccase from Pycnoporus cinnabarinus. Two substrates were used to screen the laccase activity, namely ABTS and DMP. First, the enzyme-substrate conformational space nearby the T1 pocket was sampled with Protein Energy Landscape Exploration (PELE) [56], a Monte Carlo algorithm that combines protein structure prediction algorithms with ligand and protein structural perturbations [57]. Afterwards, the lowest energy conformations were randomly selected and their reactivity was scored by estimating the amount of spin density localized on the substrate with QM/MM calculations. Authors drew three main conclusions from computational calculations: (1) the mutations accumulated during directed evolution increased laccases enzymatic activity by affecting substrate binding rather than the metal redox potential; (2) large conformational sampling revealed significant changes in the protein-ligand energy landscape upon mutation and (3) the oxidation rate of a target substrate can be improved by fine-tuning the binding event. Moreover, the oxidation (substrate’s spin density) increase as a result of an enhanced electrostatic stabilization of the radical species could be confirmed by quantum-chemical calculations. In a very recent work, researchers from same group evaluated the robustness of abovementioned computational approach to estimate activity, emphasizing the importance of the binding event in laccase reactivity [58]. To do so, they used SGZ and four phenols derivatives, with 4OH, 4MeO, 4Me and 4Cl group substituents, as substrates bounded into the T1 Cu site of Myceliophthora thermophila laccase (MtL) and P. cinnabarinus laccase (PcL). Authors concluded that redesigning substrate binding at the T1 pocket, guided by in silico methodologies, is a more consistent option that concentrated efforts on increasing the redox potential of the enzyme. The strengths and weaknesses of the protocol were also discussed as well as the significance of the methodology in protein engineering. In Section 3, two additional successful examples of protein engineering aided by computer simulations, and reported by same group, will be analyzed and discussed with more detail.

2.2. Versatile peroxidase

Computational work in VPs has been mostly focused on identifying the type of amino acid radicals that these enzymes use to degrade substrates. QM calculations at a DFT level (mainly B3LYP functional) have been used together with experimental data from multifrequency electronic paramagnetic resonance (EPR) and pulse electron nuclear double resonance (ENDOR) experiments to study the nature of the tryptophan radical (Trp164) of VP from Pleurotus eryngii [59]. Here, Pogni et al. through the calculations of g-tensors, hyperfine (hf) couplings and spin densities identified that the radical form of Trp164 most likely corresponds to a neutral radical and not a cationic radical as it was thought to be for Trp171 from LiP [60]. However, it is very likely that the tryptophan cationic radical is indeed an intermediate which then becomes deprotonated at the indole nitrogen by Glu243, which is forming a hydrogen bond at this position, to form the more stable neutral radical [59]. This study was followed by a more robust QM/MM investigation on the same system and also considering the mutant W164Y [61], where EPR parameters, structural information and spin densities were compared with experimental data to confirm the nature of the neutral tryptophan radical. The calculations using the B3LYP/CHARMM level of theory showed that deprotonation of the cation radical (though not detected experimentally) occurs through an exothermic reaction with a barrier height of approximately 2–4 kcal/mol. QM/MM MD simulations were also performed to assess the stability of the tyrosyl radical. It was observed a spontaneous proton transfer of the phenoxyl hydrogen atom to Glu243, however, the hydrogen bond formed between the phenoxyl oxygen in the W164Y variant and the newly protonated Glu243 was found to be not stable featured by the formation of new hydrogen bonds with different water molecules. The more rigid hydrogen-bond network found for the tryptophan neutral radical and the more flexible vicinity around the tyrosyl radical can be correlated with the fact that this last mutant exhibits a slower self-reduction process and oxidation of substrates, which could be attributed to a slower re-protonation of the phenoxyl oxygen to recover the reduced Tyr residue. The analysis of the electrostatic potential at this site revealed the importance of specific residues that influence the EPR parameters, which in the case of Trp164^• were found to be Glu243 and Lys253, while for Tyr^•, Glu243, Arg257 and Lys253 were identified. Following this study, in a more recent investigation [62], the flexibility and the formation of new hydrogen-bond networks in three different double variants of P. eryngii VP was studied in more detail by QM/MM MD simulations. The calculations showed that the incorporation of specific amino acids at the vicinity of the tyrosyl radical such as Leu and Glu (W164Y/R257L and W164Y/R257E mutants) helped to rigidize the hydrogen bond between the Tyr^• species with Glu243, but the variant W164Y/R257A produced the same effect that the single mutant W164Y where the hydrogen bond with Glu243 is easily broken. This can be attributed to the fact that in the first two double mutants neither Leu257 nor Glu257 compete to form hydrogen bonds with the phenoxyl oxygen of Tyr^•, as it is in the case of the single mutant through the amino acid Arg257. Besides, the inclusion of the less bulky Ala residue instead of Arg would also promote the entrance of more water molecules that facilitate the breaking of the hydrogen bond with Glu243. Thus, and in the light of these results, QM/MM methods were crucial to model these types of amino acid radicals where effects of the protein environment such as electrostatic interactions and steric effects are highly important for the stabilization of the radical species. Other recent computational studies [63] have explored the binding of veratryl alcohol (endogenous substrate of ligninolytic enzymes) at the surface of P. eryngii VP, aiming at modeling the recognition event of this molecule and then the LRET pathway to the heme group. For this, exploration of the ligand binding process through the PELE [56] method mentioned above, allowed the observation of protein-ligand complexes to initiate the study of the LRET. This latter aspect was investigated at a QM/MM level using the e-pathway approach [64], in which different residues are ranked according to their electron affinity, and in this way, the most favorable amino acids for hosting the radical can be identified. The results from this study revealed that Trp164, Trp244 and Phe198 are directly involved in the LRET, fact that was confirmed by site-directed mutagenesis experiments. More recently, Sáez-Jiménez et al. [65] reported an experimental and computational study of an evolved VP from P. eryngii that was obtained by directed evolution. The authors could provide a structural-functional basis for the improved alkaline stability of the 2-1B variant (containing E37K/H39R/V160A/T184 M/Q202L/D213A/G330R mutations). Using 20 ns MD simulations of native VP and 2-1B in solution and pH = 8.0, and starting from solved crystal structures, they could explain the structural enzyme integrity of the variant (with penta-coordinated heme iron) under alkaline condition. As alkaline inactivation has been related to the formation of a hexa-coordinated iron complex with the distal histidine residue, they monitored the evolution of that distance looking at possible correlations with the stability of the enzyme. The authors concluded that introduction of three basic residues in VP (Lys37, Arg39 and Arg330) led to new connections between heme and helix B (where the distal histidine residue is located), and formation of new electrostatic interactions, that avoided the hexa-coordination of the heme iron and thereby the alkaline inactivation of the variant.

2.3. Lignin peroxidase (LiP)

Similar to versatile peroxidase, QM/MM calculations at the level B3LYP/CHARMM have been carried out to elucidate the nature of the tryptophan radical (Trp171) present in LiP from P. chrysosporium [60]. These calculations were based on the experimental work performed by Smith et al. [66], where catalytic activity toward veratryl alcohol (VA) was engineered in Coprinus cinereus peroxidase (CiP). This was achieved by introducing the catalytic tryptophan residue (Trp178 in CiP) and some specific acid residues that form part of the microenvironment around it, mimicking the surroundings of Trp171 in P. chrysosporium LiP. Smith et al. also studied different LiP variants (E168Q, E250Q and E250Q/E168Q) to compare the results obtained with the new-engineered CiP. Thus, in the theoretical work of Bernini et al. [60] three systems were studied namely pristine LiP, the LiP variant E250Q/E168Q and the CiP variant D178W/R257E/R271D to allow direct comparison with the experimental data available for the last two systems. The values for g-tensors and hf coupling constants were calculated for the corresponding neutral and cationic radicals. This allowed confirming that neutral radicals (Trp171^• and Trp168^•) were present in both LiP and CiP variants, respectively. In the former case, the mutations provided a rather neutral electrostatic potential around Trp171 stabilizing the neutral radical, however, in the case of the CiP variant, since the electrostatic potential is more similar to pristine LiP, the stabilization of the radical comes from the presence of K260, residue that is absent in pristine LiP. This later residue would facilitate the entrance of more water molecules around the radical for its proper stabilization through hydrogen-bonding interactions. No experimental data are available in the case of pristine LiP to establish the nature of the Trp171 radical, though QM/MM calculations showed that deprotonation of Trp171 would be energetically unfavorable and that Trp171^•+ might exist in pristine LiP due to the highly negative electrostatic potential around it [60]. QM and QM/MM calculations have also been carried out to understand the binding and the electron transfer mechanisms of VA at Tyr181 (redox-active amino acid) in Trametopsis cervina LiP [67]. It was observed that VA is stabilized by sandwich π stacking interactions with Phe89 and Tyr181 showing a spin density distribution largely maintained on the VA ligand. QM/MM calculations including in the QM region the heme group and its axial ligands, VA, Tyr181 and Phe89, showed that two unpaired electrons are present in the iron-oxo moiety while another unpaired electron is shared among the porphyrin ring, Tyr181 and VA, observing, therefore, electron transfer from the Tyr181-VA complex to the porphyrin radical, giving rise to compound II. An interesting characteristic of T. cervina LiP is the formation of adduct between Tyr181 and VA through covalent bonding. Besides experimental findings, its existence was confirmed by EPR-based quantum calculations and it was further characterized by QM calculations to study the most favorable binding position of VA at the ring of Tyr181. On the other hand, in a similar research work [68], the binding of VA at Trp171 in P. chrysosporium LiP was studied with different computational methods and QM/MM calculations were used to compute spin densities for different VA-Trp171^•+ complexes at different separation distances. The goal was to obtain a cut-off distance for the occurrence of the electron transfer between these two fragments. This information was then used to analyze MD simulations of close-distance complexes that allowed identifying crucial interactions for the binding of VA at Trp171. Other recent QM/MM studies [69] have been focused on calculating the redox potential of compound I in LiP, looking for identifying specific amino acids that influence its oxidative power, and in this way propose new mutants with higher oxidative abilities, or with the capacity of functioning at softer pH conditions. Castro et al. [69] identified that acidic (mainly Asp and Glu) residues around the oxoferryl complex tend to lower the redox potential by electrostatic stabilization while positively charged residues destabilize the complex and therefore increase the redox potential. Two amino acids were in silico mutated (E40Q and D183N) which were chosen based on electrostatic and structural considerations and gave estimates of redox potential increases of 140 and 190 mV, respectively, and therefore represent novel mutations that can be tested experimentally. In an alternative approach, experimental work together with QM calculations have been used to engineer new mutants of P. chrysosporium LiP (isoenzyme H8) that increase the lignin-degrading efficiency by stabilizing radicals that relay the electron transfer to the heme [70]. Pham et al. found that Trp251 is a site where radical products (guaiacol radicals) from a non-phenolic lignin dimer are covalently bound, causing therefore suicidal inhibition of the degradation process. DFT-based QM calculations were used to calculate ΔG° barriers for the critical redox centers H176/Heme, W171 and W251, which confirmed W251 as a possible electron relay amino acid in the LRET process. Thus, mutants such as A242D and T208D that would tend to stabilize the W251 radical by electrostatic and hydrogen-bonding interactions were proposed and characterized, with the former one exhibiting until 21.1- and 4.9-fold higher increments in k_cat and k_cat/k_M values, respectively. In order to get more insights into the electron transfer mechanism operating in LiP, QM/MM calculations to analyze the most favorable electron transfer pathway have been performed using the e-pathway approach, as was done for versatile peroxidase [63]. It was confirmed the involvement of Trp251 in the electron transfer to the heme, proposing the final route Trp171-Phe205-Trp251 to be the one involved in the LRET. The computational results were supported by the experimental mutation W251A, which provided threefold lower catalytic efficiency toward the oxidation of VA. This finding allowed concluding that both VP and LiP, which come from different basidiomycetes, use an analog pathway in the LRET. Some research works have been done with regard to application of other computational methods to understand the stability of wild type LiP and mutants or binding of lignin compounds to its heme group or reactive surface Trp residue. For example, Gerini et al. investigated the dynamical and structural properties of LiP and its W171A mutant in aqueous solution using MD simulations. They found that in both cases, the enzyme retained its overall backbone structure and all its non-covalent interactions in the course of the MD simulations. Moreover, these authors performed steered molecular dynamics docking simulations which have shown that LiP natural substrate (VA) could easily approach the heme edge through the access cannel [71]. With a similar goal, some molecular dynamics calculations were performed on LiP and horseradish peroxidase (HRP) structures in order to evaluate the effect of calcium ions on protein structure [72]. Chen et al.determined and verified, by MD simulations, the robustness and stability of the docking predicted 3D structures of laccase-lignin, LiP-lignin and MnP-lignin complexes [73]. From the point of view of electronic structure methods, Elder et al. used ab initio molecular orbital calculations for studying the LiP active site and VA at the STO-3G level using the unrestricted Hartree-Fock approximation for open-shell species [74]. In that book chapter, the author carefully studied the electronic and energetic properties of reactant species in LiP active site (Compounds I and II and VA^•+ as well as their complexed forms). Due to the general biological significance of the metalloporphyrins and the enzyme systems in which they occur, including LiP, the author made a revision of all previous reports where these compounds have been studied using the methods of theoretical chemistry. Earlier the same author reported the oxidation of a lignin model compound via electron transfer. They found that oxidation of VA appears to be controlled by relative endothermicity associated with a proton transfer step within each of the possible oxidative pathways. The transfer of an electron from the lignin model compound to the VA^•+was slightly endothermic, occurring by way of a non-bonded intermediate, with slightly lower energy than either products or reactants [75]. Finally, ten Have et al. used calculated ionization potentials (IP) to determine the oxidation of vanillin precursors by LiP. The conversion of a series of phenolic derivatives, O-methyl ethers, to the Cα-Cβ cleavage product 3,4-dimethoxybenzaldehyde by LiP allowed them to hypothesize that similar compounds with an IP < 9.0 eV would undergo the same Cα-Cβ cleavage reaction upon incubation with LiP [76].

2.4. Manganese peroxidase

Computational research in manganese peroxidase has being centered on the study of the binding interactions of this enzyme with lignin model compounds as well as the role played by the C-terminal tail in the catalytic differences seen among MnPs subfamilies. For the former analysis, computational studies using molecular docking have been carried out to describe the binding modes of lignin (using a lignin model substrate) to LiP, laccase and MnP from P. chrysosporium, and MD simulations were performed to verify the stability of these systems. Chen et al. [73] determined for the lignin-MnP system that in the best docking pose the lignin model compound was positioned at the center of the binding pocket found by the cavity detection algorithm of Molegro Virtual Docker (MVD). The most important residues for the binding of lignin were Arg42, His173 and Arg177, all located next to the heme group. Hydrophobic contacts were the leading interactions. It is worth noting that according to the crystal structure of MnP [77], the residues Glu35, Glu39 and Asp179 are also implied in the binding of Mn²⁺ and the residue Arg42 could be involved in the stabilization of this ion through its interaction with Glu39. MD simulations of 3000 ps for this system were performed showing stability with respect to the starting conformation after 1000 ps, with a mean RMSD of 2.4 and 1.4 Å for the backbone of MnP and lignin, respectively. Also, the total energy during the molecular dynamics remained without major variations, confirming the stability of the system. The role of the C-terminal tail in the catalysis and stability of MnP from Ceriporiopsis subvermispora was studied through dynamic ligand diffusion and QM/MM techniques [78]. In a study of over 31 fungal genomes [79], 3 MnP subfamilies were defined depending on the length of the C-terminal tail: short, long and extralong MnPs. Experimental results showed that extralong tails of MnPs contribute to the Mn²⁺ oxidation and to a major acidic stability, but also showed that only short MnPs have the capacity of oxidizing ABTS independently (in the absence of Mn²⁺) [78]. Fernández-Fueyo et al. using PELE [56] and electronic couplings calculations by the fragment charge difference (FCD) approach [80], provided computational information for the ABTS oxidation at the heme-propionate channel and the relation between the C-tail length and the MnPs catalytic properties. Since only small differences exist between long and extralong MnPs, there was no representative for the long MnP subfamily. Calculations were performed using the crystal structure of the extralong isoenzyme MnP6 [78] and its in silico shortened-tail form. Local exploration of the ABTS diffusion at the entrance of the heme-propionate channel using PELE showed that for the shortened-tail form of MnP6 the most favorable structures had ABTS closer to the active site. In contrast, for the extralong isoenzyme, the tool failed to provide a structure with good interaction energy, and ABTS was located at a longer distance from the active site. It was observed that the C-tail interferes with the entrance of the ligand into the heme-propionate channel and that longer Gluh35-Asph179 distances enable the approach of ABTS in the shortened-tail form. The C-tail in the extralong MnP seems to limit the mobility of the residues at the entrance of the channel, and also Asph179 and Gluh35 block the access of ABTS even in the most energetically favourable position. Moreover, QM/MM calculations (B3LYP/OPLS2005) using the FCD approach were performed in the best positions predicted for both systems, showing electronic coupling values of one order of magnitude higher for the shortened-tail MnP (4.315E−5 eV) than for the extralong MnP (4.65E-6 eV), demonstrating the higher capacity of the shortened-tail MnP to oxidize ABTS. Rational enzyme engineering of MnP6 from C. subvermispora has also been carried out with the help of computational methods [81]. Here, this highly stable peroxidase at low pH was rationally mutated and activated for the oxidation of ABTS based on the analysis of MnP4, a less stable but ABTS-oxidizing MnP from Pleurotus ostreatus [82]. Acebes et al. started the engineering research using PELE for the active site inspection on both systems, ABTS-MnP6 and ABTS-MnP4. Results of these explorations showed that in the energetically minimum structure of MnP4 located at the main heme channel, two histidines, H220 and H142, interacted forming hydrogen bonds with ABTS’s negatively charged sulfonates. These two histidines were not found in MnP6, but instead G139 and N218 were found in their respective positions. Thus, in silico G139H and N218H mutations were built in MnP6. Further PELE exploration of the ligand diffusion of ABTS on MnP4 and MnP6 double mutant showed similar energy values, distance and binding modes between these two systems. These results were also supported by electronic coupling calculations using the FCD method, where this double mutant presented an increase of the electronic coupling value of almost three orders of magnitude compared to wild type MnP6 and about ~four times lower than in MnP4. Experimental results confirmed the activation of MnP6 toward the oxidation of ABTS.

2.5. Dye-decoloring peroxidase

DyP-type peroxidases were described 10 years ago like a new family of the heme-peroxidases [83]. The structural difference of them with the other class of fungal, plant and animal peroxidases is that they do not have well conserved amino acids around the heme moiety. Initial studies in DyPs exhaustively characterized this region and after being reported the cloning and expression of Auricularia auricula-judae DyP (AauDyP) in E. coli [84], studies of the substrate oxidation sites in AauDyP were reported. In this context, computational studies have helped to obtain a better characterization of AauDyP and to give molecular details in the different applications that have been explored in this family of peroxidases. In a first approximation, simulations using PELE [56] were performed to study the binding of the substrate RB19 to the structure of AauDyP [85]. The main interactions between RB19 and DyP were found to be with residues Trp105, Tyr147/337, Trp207, Tyr285, Trp377 and the heme entrance channel. The investigation aimed at elucidating which of these residues were involved in the oxidation of the substrate. For this, QM/MM pairwise comparison between Trp377 and the others residues was carried out where spin densities were calculated after subtracting one electron from the systems. Other QM/MM calculations including the heme group in its compound I state and RB19 in the QM region showed that Trp377 is the preferential oxidation site due to a higher concentration of the spin density on this residue. However, analysis of the spin density also showed that other surface amino acids such as Tyr337 could also function as potential oxidizing sites but with lower activities [85]. This was confirmed by site-directed mutagenesis experiments and enzymatic kinetic characterization of the mutants W377S, Y147S, Y337S, Y147F/Y337F and G169L, which revealed the existence of two oxidation sites, a high turnover site (residue Trp377) and a low turnover site (probably the heme-access channel). The data from EPR experiments proved the existence of the radical species Tyr337^• and Trp377^• and encouraged to propose that the oxidation in DyP-type peroxidases of bulky substrates, such as RB19, is ruled by a long-range electron transfer (LRET) mechanism. Finally, by means of the QM/MM e-pathway approach the residues Trp377, Pro310, Arg309, Arg306, Ile305 and His304, were identified as the residues involved in the LRET pathway in A. aricula-judae DyP [85]. Other studies of the mutants W377S, Y147S, Y337S and Y147F/Y337F in A. aricula-judae DyP evaluated by Baratto et al. [86] through EPR spectroscopy and QM/MM calculations have been performed. The calculation of g-tensor values via single point calculations on the B3LYP/AMBER-optimized geometries of the Tyr337^• and Tyr147^• radicals provided new insights into the molecular interactions with their surroundings. They found that Tyr337 presents a strong hydrogen-bonding interaction between the phenoxyl oxygen and the phenoxyl proton transferred to the carboxylic acid oxygen of the nearby Glu354. On the other hand, Tyr147 is involved in weaker hydrogen-bonding interactions with solvent molecules and the Trp377 radical is present in its neutral form, similar to Trp164 in P. eryngii VP. Additional evidence from EPR spectroscopy experiments helped to conclude that Tyr337 is the secondary radical species formed in A. aricula-judae DyP at pH 3, which is simultaneously created with the Trp377 radical, which is the main radical site. Thus, if Trp377 is replaced, Tyr337 constitutes the main redox site, and once Tyr337 is removed, a potential LRET pathway originated from Tyr147 is activated. This fact evidenced the activation and deactivation of different LRET pathways in A. aricula-judae DyP [86]. Other more recent studies [87], have also attempted to characterize the oxidative sites present in AauDyP by means of molecular dynamics simulations, where energy differences in the charge transfer processes among the oxidative centers have been identified. Trp377 is recognized as the most kinetically favorable amino acid for hosting the electron hole, confirming the previously mentioned findings, otherwise the residue Tyr229, or the cluster formed by Trp105-Tyr147-Tyr337, are also thermodynamically favored. Based on the structural information provided by the crystal structure of A. aricula-judae DyP, it was proposed as a protein engineering strategy, the enlargement of the heme pocket by the single mutations F359G and L357G, and the evaluation of the sulfoxidizing reactivity of mutants toward methyl-phenyl sulfide (MPS) and methyl-p-tolyl sulfide (MTS). It is worth mentioning that the tuning of the stereoselectivity in sulfoxidation reactions has many applications in organic synthesis and in enzymatic reactions. The mutant F359G was observed to be highly stereoselective and reactive producing the S sulfoxide of MPS and MTS, while native DyP had no sulfoxidation activity, and the L357G variant produced racemic mixtures of sulfoxides R and S enantiomers but with lower catalytic efficiencies. PELE simulations showed that the sulfoxidation ability present in the mutants was due to a better access of the substrates to the reactive heme cofactor. In addition to this, QM/MM calculations performed at the DFT M06-L(lacvp*)/OPLS level with QSite [88] elucidated that the high stereoselectivity and reactivity observed in the mutant F359G could be explained by the substantial difference between the spin density populations on the substrate molecules for the two enantiomers, (0.8% pro R-enantiomer, 21.8% pro S-enantiomer) in MPS and (12.5% pro R-enantiomer, 25.8% pro S-enantiomer) in MTS [89]. Other study including computational methodologies was reported by Strittmatter et al. [24]. Docking experiments were performed in AauDyPI using N-heterocyclic derivatives like substrates to explore its binding interaction in the heme pocket. The observation of competitive inhibition by imidazole with regard to H₂O₂ conversion and changes in the Soret band of AauDyPI in the presence of imidazole derivatives was correlated with the computational data and permitted to predict which compound was a good candidate for coordination or not with the heme iron [24].

3. Rational protein design successful cases

In this section, we aim to briefly present two relevant and successful cases of rational protein design, where computational tools played a key role in improvement of catalytic activity of laccase as a key representative of oxidoreductase family. The first case was reported in 2016 by Santiago et al. [90]. That research work described the computer-aided engineering of a double mutant laccase that showed improved kinetic parameters (k_cat and k_M), when compared to parent enzyme, for oxidation of several arylamines. The authors studied the biological oxidation of four substrates (aniline (ANL), p-phenylenediamine (PPD), ABTS and N,N-dimethyl-p-phenylenediamine (DMPD)) that are of industrial interest for synthesis of conductive polyaniline (PANI), in case of ANL, or precursors of high-value reagents like methylene blue as potential antimalarial agent, in case of DMPD. The computational protocol used to design the double mutant of laccase was as follows: (1) determining of the binding event of substrates nearby T1 site by means of Monte Carlo-based molecular simulations using software PELE [56]; (2) estimating the spin density fraction localized on the substrate by means of QM/MM calculations on randomly selected structures; and (3) designing the laccase-binding site in order to improve substrate oxidation rate. This computational protocol was validated experimentally by estimating the kinetic parameters of substrates and comparing them with parent laccase. Notably, in case of ANL a 2-fold k_cat increase was achieved by introducing two mutations at T1 site namely N207S/N263D that stabilized the oxidized form of ANL thereby boosting electron transfer to the Cu T1 site. Moreover, the new mutant was tested at working conditions for the enzymatic synthesis of PANI from ANL, which showed an increment of absorbance of 36% compared to the parent when production of PANI was followed at 800 nm. Finally, an in silico cross validation was performed for aniline derivative DMPD. After applying the same computational protocol and measuring the kinetic parameters for biological oxidation of this substrate, by mutant and parent laccases, the authors found that substrate spin density increased from 39 to 68% by the presence of N207S/N263D mutations in the engineered laccase which also was consistent with the observed increase in k_cat (from 459 ± 18 to 741 ± 48 s⁻¹). The authors concluded that applied computational approach provided a reliable in silico screening, reducing experimental validation to precisely design laccase mutants.

In a second case of rational protein engineering, the substrate binding pocket of a chimeric laccase (3A4) was redesigned in order to improve the oxidation of synapic acid (SA), a lignin-related phenol of industrial interest [91]. This time the authors started using an experimentally protocol, iterative saturation mutagenesis [92], for screening more than 15,000 clones, from which two variants (C14F12 and CA32F1) and parent 3A4 were purified and characterized. Then, they rationalized the effect of mutations using computational methods. To do so, the 3D structures for 3A4 (mutation V162A), C14F12 (mutations V162A and F392N) and CA32F1 (mutations V162R, T164E and F392N) were built from coordinates of PM1L (PDB code: 5ANH), which shares 98% sequence identity with the chimeric variant 3A4. The 3D models were energy minimized with a force field and, in case of triple mutant variant, a 20 ns MD simulation was necessary to assure the structural stability of the model. With respect to SA, it was modeled in its protonated and deprotonated form, SAH and SA⁻, respectively, due to its pK_a = 4.9 and MD simulations and experiments were performed at pH = 5.0 and 3.0. Then the binding modes of SA, and other substrates, on the T1 site of chimeric proteins were sampled with PELE [56]. Finally, and after filtering structures based on their energy from PELE simulations, the atomic Mulliken spin density of the substrate was computed on randomly selected structures that were minimized with QM/MM methodology including in the QM part the substrate, the T1 copper atom, and residues H394, C450, H455, I452 and F460; meanwhile the rest of structure was modeled classically with OPLS force field. The main results from computational simulation, and their relation with experimental findings, in this research work can be resumed as follow: (1) In 3A4, SAH binds with high affinity inside T1 site interacting with residues P393 and N207 through its carboxylic acid moiety and with residues A162, P163 and the side chain of S264 through its phenolic group. Its average spin density was about 88%. On the other hand, SA⁻ showed fewer structures populating the T1 site and instead a second site was explored. In that site substrate bound in a different orientation and further away from the metal center due to the favorable interaction between the negatively charged carboxylate group, the backbone of G391 and F392 and the hydroxyl group of S387. The SA⁻ averaged spin densities in those sites were 87 and 61%, respectively. That roughly means that 3A4 laccase oxidized, preferentially, SAH instead of SA⁻. (2) In C14F12, SA⁻ binds with more affinity to T1 site than SAH showing the same interactions of its carboxylic group with residues G391, N392 and S387’s side chain. Notably, the mutation F392N created a favorable binding site for the negatively charged substrate with the N392 side chain; therefore substrate’s phenolic group showed two main orientations (one inside T1 site and other pointing to the solvent) with improved spin densities of 93 and 48%, respectively. This, in turn, was in agreement with improvement of catalytic constants for C14F12 (k_cat = 251 s⁻¹) relative to 3A4 (k_cat = 156 s⁻¹) that could be due to improved oxidation of SA⁻ (that is abundant at this pH than at pH = 3). (3) In CA32F1, the SAH and SA⁻ species were also bound with high affinity to T1 site, showing spin densities of 99 and 86%, respectively. However the relative orientation of the side chains from R162, E164 and D205, and specially the new acid–basic pair formed by R162 and E164, affected the substrate oxidation. Overall, the authors concluded that QM and MD calculations were performed to rationalize the effect of the selected mutations, revealing the critical role of the residues of the enzyme pocket in order to provide the precise binding of the substrate that enabled an efficient electron transfer to the T1 Cu. Additionally, they stressed the usefulness of atomic simulations for unveiling the molecular determinants for the efficient oxidation of a target molecule by these kind of enzymes.

4. Perspectives and conclusions

This chapter aimed to collect and highlight the increasing evidence of computational modeling methods as supporting and validation tools in the study of oxidoreductases as biocatalysts and the characterization of their biotransformation products. We decided to focus on five enzymes that have been investigated intensively during last years due to their broad application as biocatalysts in chemical, agriculture, environmental and biotechnological applications. Those enzymes are LiP, MnP and VP, which are heme peroxidases in the peroxidase-catalase superfamily, DyP as a recently discovered heme peroxidase in the CDE superfamily (including chlorite dismutases, DyPs and the EfeB proteins of E. coli), and laccase as the most studied member of copper-containing oxidoreductases. The relevant literature was intensively searched for reports were computational methods—mainly MD simulations, QM and QM/MM—were presented as auxiliary tools to explain the experimental evidence in aspects like substrate binding and oxidation of lignin or derivatives against different binding sites reported for these enzymes, improvement of catalytic activity of a given oxidoreductase for a specific substrate, rational protein design and more recently computed-aided protein engineering. However, we believe that, despite the recent advances in hardware, software and algorithms, the use of computational modeling methods in the understanding of substrate binding and enzymatic catalytic processes between oxidoreductases and their substrates is yet scarce. Nevertheless, the two recently successful and relevant applications, that were discussed in last section, about the interplay between experimental and computational protocols applied to directed evolution of improved laccases as biocatalysts are very encouraging. This accumulated knowledge in rational protein design from presented case studies on laccase paves the way to other studies of the same kind in enzymes from oxidoreductases.

Undoubtedly, during next years the research in this area will be benefited from advance in computational software and hardware as well as from more accurate and predictive computational tools. Additionally, it is expected that more efficient biocatalysts are successfully and quickly generated from conjunction between these experimental/computational techniques presented here with the aim to improve catalytic processes that are urgently demanded in environmental, food, chemical, pharmaceutical and energy industries.

Acknowledgments

J. A-M. and I. F-V. thank financial support from project FONDECYT No. 1140618. J.R. and R.R. acknowledge support from CONICYT to perform their doctoral studies through scholarships No. 21160905 and No. 21130949, respectively.