Protonation states of methionine aminopeptidase and their relevance for inhibitor binding and catalytic activity

We have performed a computational study of different protomeric states of the methionine aminopeptidase active site using a combined quantum-mechanical/molecular mechanical simulation approach. The aim of this study was to clarify the native protonation state of the enzyme, which is needed for the development of novel irreversible inhibitors that can possibly be used as antiangiogenic and antibiotic drugs by virtual screening and other drug design methods. The results of the simulations indicated that two protonation states are possible without disturbing the overall geometry of the active site. We then verified experimentally the presence of the two protonation states by studying the substrate hydrolysis and inhibitor binding reactions at different pH values and come to the conclusion that one of the protomeric states is relevant for inhibitor binding, whereas the other is relevant for substrate hydrolysis. This result has implications for the development of other inhibitors of this class of enzymes and adds a new perspective to the pharmacological properties of the antiangiogenic drug fumagillin, which is an irreversible inhibitor of the human methionine aminopeptidase type II.     

Comparative molecular dynamics simulations of HIV-1 integrase and the T66I/M154I mutant: binding modes and drug resistance to a diketo acid inhibitor

HIV-1 IN is an essential enzyme for viral replication and an interesting target for the design of new pharmaceuticals for use in multidrug therapy of AIDS. L-731,988 is one of the most active molecules of the class of beta-diketo acids. Individual and combined mutations of HIV-1 IN at residues T66, S153, and M154 confer important degrees of resistance to one or more inhibitors belonging to this class. In an effort to understand the molecular mechanism of the resistance of T66I/M154I IN to the inhibitor L-731,988 and its specific binding modes, we have carried out docking studies, explicit solvent MD simulations, and binding free energy calculations. The inhibitor was docked against different protein conformations chosen from prior MD trajectories, resulting in 2 major orientations within the active site. MD simulations have been carried out for the T66I/M154I DM IN, DM IN in complex with L-731,988 in 2 different orientations, and 1QS4 IN in complex with L-731,988. The results of these simulations show a similar dynamical behavior between T66I/M154I IN alone and in complex with L-731,988, while significant differences are observed in the mobility of the IN catalytic loop (residues 138-149). Water molecules bridging the inhibitor to residues from the active site have been identified, and residue Gln62 has been found to play an important role in the interactions between the inhibitor and the protein. This work provides information about the binding modes of L-731,988, as well as insight into the mechanism of inhibitor-resistance in HIV-1 integrase.     

Determinants of reactivity and selectivity in soluble epoxide hydrolase from quantum mechanics/molecular mechanics modeling

Soluble epoxide hydrolase (sEH) is an enzyme involved in drug metabolism that catalyzes the hydrolysis of epoxides to form their corresponding diols. sEH has a broad substrate range and shows high regio- and enantioselectivity for nucleophilic ring opening by Asp333. Epoxide hydrolases therefore have potential synthetic applications. We have used combined quantum mechanics/molecular mechanics (QM/MM) umbrella sampling molecular dynamics (MD) simulations (at the AM1/CHARMM22 level) and high-level ab initio (SCS-MP2) QM/MM calculations to analyze the reactions, and determinants of selectivity, for two substrates: trans-stilbene oxide (t-SO) and trans-diphenylpropene oxide (t-DPPO). The calculated free energy barriers from the QM/MM (AM1/CHARMM22) umbrella sampling MD simulations show a lower barrier for phenyl attack in t-DPPO, compared with that for benzylic attack, in agreement with experiment. Activation barriers in agreement with experimental rate constants are obtained only with the highest level of QM theory (SCS-MP2) used. Our results show that the selectivity of the ring-opening reaction is influenced by several factors, including proximity to the nucleophile, electronic stabilization of the transition state, and hydrogen bonding to two active site tyrosine residues. The protonation state of His523 during nucleophilic attack has also been investigated, and our results show that the protonated form is most consistent with experimental findings. The work presented here illustrates how determinants of selectivity can be identified from QM/MM simulations. These insights may also provide useful information for the design of novel catalysts for use in the synthesis of enantiopure compounds.     

Investigation of the acylation mechanism of class C beta-lactamase: pKa calculation, molecular dynamics simulation and quantum mechanical calculation

β-Lactamases are bacterial enzymes that act as a bacterial defense system against β-lactam antibiotics. β-Lactamase cleaves the β-lactam ring of the antibiotic by a two step mechanism involving acylation and deacylation steps. Although class C β-lactamases have been investigated extensively, the details of their mechanism of action are not well understood at the molecular level. In this study, we investigated the mechanism of the acylation step of class C β-lactamase using pKa calculations, molecular dynamics (MD) simulations and quantum mechanical (QM) calculations. Serine64 (Ser64) is an active site residue that attacks the β-lactam ring. In this study, we considered three possible scenarios for activation of the nucleophile Ser64, where the activation base is (1) Tyrosine150 (Tyr150), (2) Lysine67 (Lys67), or (3) substrate. From the pKa calculation, we found that Tyr150 and Lys67 are likely to remain in their protonated states in the pre-covalent complex between the enzyme and substrate, although their role as activator would require them to be in the deprotonated state. It was found that the carboxylate group of the substrate remained close to Ser64 for most of the simulation. The energy barrier for hydrogen abstraction from Ser64 by the substrate was calculated quantum mechanically using a large truncated model of the enzyme active site and found to be close to the experimental energy barrier, which suggests that the substrate can initiate the acylation mechanism in class C β-lactamase.     

Elucidation of the protonation states of the catalytic residues in mtKasA: implications for inhibitor design

KasA (β-ketoacyl ACP synthase I) is involved in the biosynthetic pathway of mycolic acids, an essential component of the cell wall in Mycobacterium tuberculosis. It was shown that KasA is essential for the survival of the pathogen and thus could serve as a new drug target for the treatment of tuberculosis. The active site of KasA was previously characterized by X-ray crystallography. However, questions regarding the protonation state of specific amino acids, the orientation of the histidine groups within the active site, and additional conformers being accessible at ambient temperatures remain open and have to be addressed prior to the design of new inhibitors. We investigate the active site of KasA in this work by means of structural motifs and relative energies. Molecular dynamics (MD) simulations, free energy perturbation computations, and calculations employing the hybrid quantum mechanics/molecular mechanics (QM/MM) method made it possible to determine the protonation status and reveal important details about the catalytic mechanism of KasA. Additionally, we can rationalize the molecular basis for the acyl-transfer activity in the H311A mutant. Our data strongly suggest that inhibitors should be able to inhibit different protonation states because the enzyme can switch easily between a zwitterionic and neutral state.     

Molecular dynamics simulations of arachidonic acid complexes with COX-1 and COX-2: insights into equilibrium behavior

The cyclooxygenase (COX) enzymes are responsible for the committed step in prostaglandin biosynthesis, the generation of prostaglandin H(2). As a result, these enzymes are pharmacologically important targets for nonsteroidal antiinflammatory drugs, such as aspirin and newer COX-2 selective inhibitors. The cyclooxygenases are functional homodimers, and each subunit contains both a cyclooxygenase and a peroxidase active site. These enzymes are quite interesting mechanistically, as the conversion of arachidonic acid to prostaglandin H(2) requires two oxygenation and two cyclization reactions, resulting in the formation of five new chiral centers with nearly absolute regio- and stereochemical fidelity. We have used molecular dynamics (MD) simulations to investigate the equilibrium behavior of both COX-1 and COX-2 enzyme isoforms with bound arachidonate. These simulations were compared with reference simulations of arachidonate in solution to explore the effect of enzyme on substrate conformation and positioning in the active site. The simulations suggest that the substrate has greater conformational freedom in the COX-2 active site, consistent with the larger COX-2 active site volume observed in X-ray crystal structures. The simulations reveal different conformational behavior for arachidonate in each subunit over the course of extended equilibrium MD simulations. The simulations also provide detailed information for several protein channels that might be important for oxygen and water transport to or from active sites or for intermediate trafficking between the cyclooxygenase and peroxidase active sites. The detailed comparisons for COX-1 versus COX-2 active site structural fluctuations may also provide useful information for design of new isozyme-selective inhibitors.     

In silico analyses of substrate interactions with human serum paraoxonase 1

Human paraoxonase (HuPON1) is a serum enzyme that exhibits a broad spectrum of hydrolytic activities, including the hydrolysis of various organophosphates, esters, and recently identified lactone substrates. Despite intensive site-directed mutagenesis and other biological studies, the structural basis for the specificity of substrate interactions of HuPON1 remains elusive. In this study, we apply homology modeling, docking, and molecular dynamic (MD) simulations to probe the binding interactions of HuPON1 with representative substrates. The results suggest that the active site of HuPON1 is characterized by two distinct binding regions: the hydrophobic binding site for arylesters/lactones, and the paraoxon binding site for phosphotriesters. The unique binding modes proposed for each type of substrate reveal a number of key residues governing substrate specificity. The polymorphic residue R/Q192 interacts with the leaving group of paraoxon, suggesting it plays an important role in the proper positioning of this substrate in the active site. MD simulations of the optimal binding complexes show that residue Y71 undergoes an "open-closed" conformational change upon ligand binding, and forms strong interactions with substrates. Further binding free energy calculations and residual decomposition give a more refined molecular view of the energetics and origin of HuPON1/substrate interactions. These studies provide a theoretical model of substrate binding and specificity associated with wild type and mutant forms of HuPON1, which can be applied in the rational design of HuPON1 variants as bioscavengers with enhanced catalytic activity.     

Insights into the structure and dynamics of the dinuclear zinc beta-lactamase site from Bacteroides fragilis

Herein, we report quantum chemical calculations and molecular dynamics (MD) simulations of the dinuclear form of the Bacteroides fragilis zinc beta-lactamase. We studied four different configurations which differ in the protonation state of the Asp103 residue and in the presence or absence of a Zn1-OH-Zn2 bridge. The flexibility of the Zn1-OH-Zn2 bridge was studied by means of quantum mechanical (QM) calculations on cluster models while the relative stabilities of the different configurations were estimated from QM linear scaling calculations on the enzyme. Contacts between important residues (Cys104, Asp69, Lys185, etc.), the solvation of the zinc ions, and the conformation of the active site beta-hairpin loop were characterized by the MD analyses. The influence of the buried sodium ion close to the Zn2 position was investigated by carrying out a secondary simulation where the sodium ion was replaced with an internal water molecule. The comparative structural analyses among the different MD trajectories augmented with energetic calculations have demonstrated that the B. fragilis protein efficiently binds the internal Na(+) ion observed crystallographically. Moreover, we found that when Asp103 is unprotonated, a rigid Zn1-OH-Zn2 bridge results, while for neutral Asp103, a fluctuating Zn1-Zn2 distance was possible via the breaking and formation of the Zn1-OH-Zn2 bridge. The mechanistic implications of these observations are discussed in detail.     

Cofactor assisted gating mechanism in the active site of NADH oxidase from Thermus thermophilus

NADH oxidase (NOX) from Thermus thermophilus is a member of a structurally homologous flavoprotein family of nitroreductases and flavin reductases. The importance of local conformational dynamics in the active site of NOX has been recently demonstrated. The enzyme activity was increased by 250% in the presence of 1 M urea with no apparent perturbation of the native structure of the protein. The present in silico results correlate with the in vitro data and suggest the possible explanation about the effect of urea on NOX activity at the molecular level. Both, X-ray structure and molecular dynamics (MD) simulations, show open conformation of the active site represented by approximately 0.9 nm distance between the indole ring of Trp47 and the isoalloxazine ring of FMN412. In this conformation, the substrate molecule can bind in the active site without sterical restraints. MD simulations also indicate more stable conformation of the active site called "closed" conformation. In this conformation, Trp47 and the isoalloxazine ring of FMN412 are so close to each other (approximately 0.5 nm) that the substrate molecule is unable to bind between them without perturbing this conformation. The open/close transition of the active site between Trp47 and the flavin ring is accompanied by release of the "tightly" bound water molecule from the active site--cofactor assisted gating mechanism. The presence of urea in aqueous solutions of NOX prohibits closing of the active site and even unlocks the closed active site because of the concomitant binding of a urea molecule in the active site cavity. The binding of urea in the active site is stabilized by formation of one/two persistent hydrogen bonds involving the carbonyl group of the urea molecule. Our report represents the first MD study of an enzyme from the novel flavoprotein family of nitroreductases and flavin reductases. The common occurrence of aromatic residues covering the active sites in homologous enzymes suggests the possibility of a general gating mechanism and the importance of local dynamics within this flavoprotein family.     

Distinguishing two groups of flavin reductases by analyzing the protonation state of an active site carboxylic acid

Flavin-containing reductases are involved in a wide variety of physiological reactions such as photosynthesis, nitric oxide synthesis, and detoxification of foreign compounds, including therapeutic drugs. Ferredoxin-NADP(H)-reductase (FNR) is the prototypical enzyme of this family. The fold of this protein is highly conserved and occurs as one domain of several multidomain enzymes such as the members of the diflavin reductase family. The enzymes of this family have emerged as fusion of a FNR and a flavodoxin. Although the active sites of these enzymes are very similar, different enzymes function in opposite directions, that is, some reduce oxidized nicotinamide adenine dinucleotide phosphate (NADP(+)) and some oxidize reduced nicotinamide adenine dinucleotide phosphate (NADPH). In this work, we analyze the protonation behavior of titratable residues of these enzymes through electrostatic calculations. We find that a highly conserved carboxylic acid in the active site shows a different titration behavior in different flavin reductases. This residue is deprotonated in flavin reductases present in plastids, but protonated in bacterial counterparts and in diflavin reductases. The protonation state of the carboxylic acid may also influence substrate binding. The physiological substrate for plastidic enzymes is NADP(+), but it is NADPH for the other mentioned reductases. In this article, we discuss the relevance of the environment of this residue for its protonation and its importance in catalysis. Our results allow to reinterpret and explain experimental data.     

Molecular dynamics simulations of human butyrylcholinesterase

Herein, we present results from molecular dynamics (MD) simulations of the human butyrylcholinesterase (BuChE) enzyme in aqueous solution. Two configurations of the unbound form of BuChE differing in the presence or absence of a sodium ion inside the protein gorge were simulated for 10 and 5 ns, respectively. Besides complementing the structural information provided by X-ray data, the MD simulations give insight into the structure of the native BuChE enzyme. For example, it is shown that: the nucleophilic Ser(198) residue and the various binding subsites in the BuChE catalytic cavity are readily accessible from the exterior of the protein; the presence of the sodium ion dynamically explores two different binding sites in the gorge leading to the active site and stabilizes the productive conformation of the Glu(325)/His(438)/Ser(198) catalytic triad; several long-lived water bridges are fully integrated into the architecture of the active site; the positions of the residues at the rim of the gorge region display large deviations with respect to the crystal structure; and two side doors, constituted by residues situated at the tip of the acyl- and Omega-loops, respectively, open wide enough to allow the passage of water molecules. In conclusion, we compare our theoretical results with those from previous work on mouse acetylcholinesterase and discuss their implications for substrate binding and catalysis in BuChE.     

Crystal structure of carbapenam synthetase (CarA)

Carbapenam synthetase (CarA) is an ATP/Mg2+-dependent enzyme that catalyzes formation of the beta-lactam ring in (5R)-carbapenem-3-carboxylic acid biosynthesis. CarA is homologous to beta-lactam synthetase (beta-LS), which is involved in clavulanic acid biosynthesis. The catalytic cycles of CarA and beta-LS mediate substrate adenylation followed by beta-lactamization via a tetrahedral intermediate or transition state. Another member of this family of ATP/Mg2+-dependent enzymes, asparagine synthetase (AS-B), catalyzes intermolecular, rather than intramolecular, amide bond formation in asparagine biosynthesis. The crystal structures of apo-CarA and CarA complexed with the substrate (2S,5S)-5-carboxymethylproline (CMPr), ATP analog alpha,beta-methyleneadenosine 5'-triphosphate (AMP-CPP), and a single Mg2+ ion have been determined. CarA forms a tetramer. Each monomer resembles beta-LS and AS-B in overall fold, but key differences are observed. The N-terminal domain lacks the glutaminase active site found in AS-B, and an extended loop region not observed in beta-LS or AS-B is present. Comparison of the C-terminal synthetase active site to that in beta-LS reveals that the ATP binding site is highly conserved. By contrast, variations in the substrate binding pocket reflect the different substrates of the two enzymes. The Mg2+ coordination is also different. Several key residues in the active site are conserved between CarA and beta-LS, supporting proposed roles in beta-lactam formation. These data provide further insight into the structures of this class of enzymes and suggest that CarA might be a versatile target for protein engineering experiments aimed at developing improved production methods and new carbapenem antibiotics.     

Functional and structural insights revealed by molecular dynamics simulations of an essential RNA editing ligase in Trypanosoma brucei

RNA editing ligase 1 (TbREL1) is required for the survival of both the insect and bloodstream forms of Trypanosoma brucei, the parasite responsible for the devastating tropical disease African sleeping sickness. The type of RNA editing that TbREL1 is involved in is unique to the trypanosomes, and no close human homolog is known to exist. In addition, the high-resolution crystal structure revealed several unique features of the active site, making this enzyme a promising target for structure-based drug design. In this work, two 20 ns atomistic molecular dynamics (MD) simulations are employed to investigate the dynamics of TbREL1, both with and without the ATP substrate present. The flexibility of the active site, dynamics of conserved residues and crystallized water molecules, and the interactions between TbREL1 and the ATP substrate are investigated and discussed in the context of TbREL1's function. Differences in local and global motion upon ATP binding suggest that two peripheral loops, unique to the trypanosomes, may be involved in interdomain signaling events. Notably, a significant structural rearrangement of the enzyme's active site occurs during the apo simulations, opening an additional cavity adjacent to the ATP binding site that could be exploited in the development of effective inhibitors directed against this protozoan parasite. Finally, ensemble averaged electrostatics calculations over the MD simulations reveal a novel putative RNA binding site, a discovery that has previously eluded scientists. Ultimately, we use the insights gained through the MD simulations to make several predictions and recommendations, which we anticipate will help direct future experimental studies and structure-based drug discovery efforts against this vital enzyme.     

Atomic resolution crystal structures and quantum chemistry meet to reveal subtleties of hydroxynitrile lyase catalysis

Hydroxynitrile lyases are versatile enzymes that enantiospecifically cope with cyanohydrins, important intermediates in the production of various agrochemicals or pharmaceuticals. We determined four atomic resolution crystal structures of hydroxynitrile lyase from Hevea brasiliensis: one native and three complexes with acetone, isopropyl alcohol, and thiocyanate. We observed distinct distance changes among the active site residues related to proton shifts upon substrate binding. The combined use of crystallography and ab initio quantum chemical calculations allowed the determination of the protonation states in the enzyme active site. We show that His(235) of the catalytic triad must be protonated in order for catalysis to proceed, and we could reproduce the cyanohydrin synthesis in ab initio calculations. We also found evidence for the considerable pK(a) shifts that had been hypothesized earlier. We envision that this knowledge can be used to enhance the catalytic properties and the stability of the enzyme for industrial production of enantiomerically pure cyanohydrins.     

pH-dependent conformational flexibility of the SARS-CoV main proteinase (M(pro)) dimer: molecular dynamics simulations and multiple X-ray structure analyses

The SARS coronavirus main proteinase (M(pro)) is a key enzyme in the processing of the viral polyproteins and thus an attractive target for the discovery of drugs directed against SARS. The enzyme has been shown by X-ray crystallography to undergo significant pH-dependent conformational changes. Here, we assess the conformational flexibility of the M(pro) by analysis of multiple crystal structures (including two new crystal forms) and by molecular dynamics (MD) calculations. The MD simulations take into account the different protonation states of two histidine residues in the substrate-binding site and explain the pH-activity profile of the enzyme. The low enzymatic activity of the M(pro) monomer and the need for dimerization are also discussed.     

Molecular dynamics simulations of apo,holo, and inactivator bound GABA‐at reveal the role of active site residues in PLP dependent enzymes

The pyridoxal 5‐phosphate (PLP) cofactor is a significant organic molecule in medicinal chemistry. It is often found covalently bound to lysine residues in proteins to form PLP dependent enzymes. An example of this family of PLP dependent enzymes is γ‐aminobutyric acid aminotransferase (GABA‐AT) which is responsible for the degradation of the neurotransmitter GABA. Its inhibition or inactivation can be used to prevent the reduction of GABA concentration in brain which is the source of several neurological disorders. As a test case for PLP dependent enzymes, we have performed molecular dynamics simulations of GABA‐AT to reveal the roles of the protein residues and its cofactor. Three different states have been considered: the apoenzyme, the holoenzyme, and the inactive state obtained after the suicide inhibition by vigabatrin. Different protonation states have also been considered for PLP and two key active site residues: Asp298 and His190. Together, 24 independent molecular dynamics trajectories have been simulated for a cumulative total of 2.88 µs. Our results indicate that, unlike in aqueous solution, the PLP pyridine moiety is protonated in GABA‐AT. This is a consequence of a pK_a shift triggered by a strong charge–charge interaction with an ionic “diad” formed by Asp298 and His190 that would help the activation of the first half‐reaction of the catalytic mechanism in GABA‐AT: the conversion of PLP to free pyridoxamine phosphate (PMP). In addition, our MD simulations exhibit additional strong hydrogen bond networks between the protein and PLP: the phosphate group is held in place by the donation of at least three hydrogen bonds while the carbonyl oxygen of the pyridine ring interacts with Gln301; Phe181 forms a π–π stacking interaction with the pyridine ring and works as a gate keeper with the assistance of Val300. All these interactions are hypothesized to help maintain free PMP in place inside the protein active site to facilitate the second half‐reaction in GABA‐AT: the regeneration of PLP‐bound GABA‐AT (i.e., the holoenzyme). Proteins 2016; 84:875–891. © 2016 Wiley Periodicals, Inc.     

Investigation into mechanism of orotidine 5′-monophosphate decarboxylase enzyme by MM-PBSA/MM-GBSA and molecular docking

In this work, we studied the binding affinity of orotidine 5′-monophosphate (OMP) and 6-hydroxy-UMP (BMP) for Saccharomyces cerevisiae orotidine 5′-monophosphate decarboxylase (OMPDC) enzyme by using Molecular Mechanics-Poisson–Boltzmann Surface Area (MM-PBSA) and the Molecular Mechanics-Generalised Born Surface Area (MM-GBSA) calculations. In all simulations, Asp91, which is an important residue in the enzyme active site, was considered in both anionic (present in the native form of the enzyme) and neutral states. A series of 10-ns molecular dynamics simulations were performed for the four OMPDC–ligand complexes, two ligand-free enzymes and two free ligands, followed by MM-PBSA/MM-GBSA calculations on the collected snapshots, and molecular docking calculations using the free enzymes and ligands. The results of MM-PBSA/MM-GBSA calculations indicate that all of the OMPDC–ligand complexes form favourable systems in water, which is in agreement with corresponding experimental data. The results of the MM-PBSA and molecular docking methods also showed that OMPDC–BMP complexes, transition state analogue and inhibitor of the OMPDC enzyme have the highest binding affinities. The fact that in the native anionic state BMP shows a higher binding affinity compared with the substrate suggests the contribution of a transition state stabilisation mechanism in the debatable catalytic mechanism of the OMPDC enzyme.     

Effects of organic solvents and substrate binding on trypsin in acetonitrile and hexane media

In this work, we used molecular dynamic (MD) simulation to study trypsin with and without a six-amino-acid peptide bound in three different solvents (water, acetonitrile and hexane) in order to provide molecular information for well understanding the structure and function of enzymes in non-aqueous media. The results show that the enzyme is more compact and less native-like in hexane than in the other two polar solvents. The substrate could stabilize the native protein structure in the two polar media, but not in the non-polar hexane. There are no significant differences in the conformation of the S1 pocket upon the substrate binding in water and acetonitrile media while a reverse behavior is observed in hexane media, implying a possible induced fit binding mechanism in the non-polar media. The substrate binding enhances the stability of catalytic H-bond network since it could expel the solvent molecules from the active site. The enzyme and the substrate appear to be more appropriate to the reactive conformation in the organic solvents compared with aqueous solution. There is much greater substrate binding strength in hexane media than the water and acetonitrile ones since the polar solvent significantly weakens electrostatic interactions, which are observed to be the main driving force to the binding. In addition, some residues of the S1 pocket could remain favorable contribution to the binding despite the solvent change, but with differences in the contribution extent, the number and the type of residues between the three media.

Molecular dynamics simulations of lignin peroxidase in solution

The dynamical and structural properties of lignin peroxidase and its Trp171Ala mutant have been investigated in aqueous solution using molecular dynamics (MD) simulations. In both cases, the enzyme retained its overall backbone structure and all its noncovalent interactions in the course of the MD simulations. Very interestingly, the analysis of the MD trajectories showed the presence of large fluctuations in correspondence of the residues forming the heme access channel; these movements enlarge the opening and facilitate the access of substrates to the enzyme active site. Moreover, steered molecular dynamics docking simulations have shown that lignin peroxidase natural substrate (veratryl alcohol) can easily approach the heme edge through the access channel.