Nitrocefin hydrolysis by metallo-beta-lactamase is an important model chemical reaction mimicking cephalosporin antibiotic inactivation. Due to the specific chromogenic properties of the nitrocefin, transient kinetic data for this reaction is available. Despite its importance in the understanding of the reaction mechanism, these data can be utilised to verify benchmark calculations. This reaction is complicated from the computational viewpoint as the active site carries two double charged cations and therefore is highly polarised; nucleophilic attack and formation of the electrophilic site should be properly described. We calculate Gibbs-free energy profiles of three chemical steps comprising the entire reaction at the QM(DFT)/MM molecular dynamics level. We compare results obtained with three hybrid functionals differing in the contribution of the exact Hartree-Fock exchange, B3LYP-D3, PBE0-D3 and BHHLYP-D3. Among them, only calculations performed at the QM(PBE0-D3)/MM level were able to properly describe the intermediate accumulation and the limiting step. Laplacian of electron density maps clarify the influence of the computational protocol on the electrophilic site formation and covalent bond polarisation.

The Long Interspersed Element-1 (L1) retrotransposon is an ancient genetic parasite that comprises a significant part of the human genome. ORF2p is a multifunctional enzyme with endonuclease (EN) and reverse transcriptase (RT) activities that mediate target-primed reverse transcription of RNA into DNA. Structural studies of LINE-1 ORF2p consistently show a single Mg2+ cation in the reverse transcriptase active site, conflicting with the common DNA polymerase mechanism which involves two divalent cations. We explored a reaction pathway of the DNA elongation based on the recent high-resolution ternary complex structure of the ORF2p. The combined quantum and molecular mechanics approach at the QM (PBE0-D3/6-31G**)/MM (CHARMM) level is employed for biased umbrella sampling molecular dynamics simulations followed by umbrella integration utilized to obtain the free energy profile. The nucleotidyl transfer reaction proceeds in a single step with a free energy barrier of 15.1 ± 0.8 kcal/mol, and 7.8 ± 1.2 kcal/mol product stabilization relative to reagents. Concerted nucleophilic attack by DNA O3′ and proton transfer to Asp703 occur without a second catalytic metal ion. Estimated rate constant ∼60 s−1 aligns with RT kinetics, while analysis of the Laplacian of the electron density along the cleaving P-O bond identifies a dissociative mechanism.

The active sites of enzymes are able to activate substrates and perform chemical reactions that cannot occur in solutions. We focus on the hydrolysis reactions catalyzed by enzymes and initiated by the nucleophilic attack of the substrate’s carbonyl carbon atom. From an electronic structure standpoint, substrate activation can be characterized in terms of the Laplacian of the electron density. This is a simple and easily visible imaging technique that allows one to “visualize” the electrophilic site on the carbonyl carbon atom, which occurs only in the activated species. The efficiency of substrate activation by the enzymes can be quantified from the ratio of reactive and nonreactive states derived from the molecular dynamics trajectories executed with quantum mechanics/molecular mechanics potentials. We propose a neural network that assigns the species to reactive and nonreactive ones using the Laplacian of electron density maps. The neural network is trained on the cysteine protease enzyme-substrate complexes, and successfully validated on the zinc-containing hydrolase, thus showing a wide range of applications using the proposed approach.

ORF2p (open reading frame 2 protein) is a multifunctional multidomain enzyme that demonstrates both reverse transcriptase and endonuclease activities and is associated with the pathophysiology of cancer. The 3D structure of the entire seven-domain ORF2p complex was revealed with the recent achievements in structural studies. The different arrangements of the CTD (carboxy-terminal domain) and tower domains were identified as the "closed-ring" and "open-ring" conformations, which differed by the hairpin position of the tower domain, but the structural diversity of these complexes has the potential to be more extensive. To study this, we performed sub-microsecond all-atom molecular dynamics simulations of the entire ORF2p complex with different starting configurations. The obtained molecular dynamic trajectories frames were assigned to several clusters following the dimension reduction to three principal components of the 1275 distances feature matrix. Five and six clusters were obtained for the "open" and "closed" ring models, respectively. While the fingers-palm-thumb core retains its rigid configuration during the MD (molecular dynamics) simulations, all other domains display the complicated dynamic behavior not observed in the experimental structures. The EN (endonuclease) and CTD domains display significant translations and rotations while their internal structures stay rigid. The CTD domain can either form strong contacts with the tower or be far apart from it for both formal "open" and "closed" ring states because the tower hairpin position is not the only determining factor of the protein complex configuration. While only the "thumb up" conformation is observed in all the trajectories, the active site can be obstructed by the movement of the CTD domain. Thus, molecular modeling and machine learning techniques provide valuable insights into the dynamical behavior of the ORF2p complex, which is hard to uncover with experimental methods, given the complexity and size of the object.

Firefly bioluminescence is a product of chemical reactions that involve luciferin chromophore oxidation in the active site of luciferase proteins. We perform a series of classical molecular dynamic simulations and combined quantum and molecular mechanics (QM/MM) calculations to expose the molecular mechanism of C4 carbon atom deprotonation in luciferyl adenylate molecule. QM/MM calculations confirm that ND-protonated His245 residue is a suitable proton acceptor in the wild type Photinus pyralis luciferase. Classical molecular dynamic simulations reveal oxygen binding cavities inside the protein including the one located close to the C4 atom of luciferin. In mutant forms that lack direct interactions with the His245 side chain, a proton wire comprising water molecules stimulates either protonation of the phosphate group of the luciferyl adenylate forming an unstable intermediate or keto-enol tautomerization of luciferin. The comparison of the ionization potentials of molecular systems with the ‘deprotonated’ C4 carbon atom revealed that the ionization energy for the enol tautomeric form is close to the system with the His245 proton acceptor. Thus, existence of the keto-enol tautomerization channel might explain bioluminescence in case of absence of the amino acid proton acceptor.

N-acetylaspartilglutamate is the most common dipeptide in brain cells, which is synthesized using the enzyme N-acetylaspartilglutamate synthase. Herein we utilize bioinformatics methods to predict the protein structure from the primary sequence of the coding gene, classical molecular dynamics to obtain a stable protein complex with N-acetylaspartate and glutamate ligands within the trajectory, as well as machine learning methods to analyze, describe and select potential reactive and non-reactive conformations of the model system describing the enzyme-substrate complex. Molecular dynamics simulations with combined quantum mechanics / molecular mechanics potentials were performed for a set of selected conformations and the potential reaction mechanism were characterized.

The combined quantum mechanics/molecular mechanics method is most often used to describe the molecular mechanisms of enzymatic reactions. The review discusses the main methodological issues, gives practical recommendations, and also illustrates the progress of the method over the past 20 years using an important example of the reaction of guanosine triphosphate hydrolysis by a protein complex.

The high-order Rayleigh-Schr & ouml;dinger perturbation theory (RSPT) can be applied for studying anharmonic vibrational problem formulated with the isomorphic Hougen Hamiltonian, but the resulting series usually possess slowly convergent or even divergent behavior. This flaw can be overcome by the resummation of such series with the multi-valued Hermite-Pad & eacute; approximant (HPA) that accurately reproduces variational matrix eigenvalues provided the basis set is the same. Besides, the state-to-state juxtaposition of HPA branch points can provide an accurate quantitative description of resonance phenomena. Such resummation was earlier proven to be efficient for three- and four-atomic asymmetric top molecules, as well as for linear molecules CO2 2 and C2H2. 2 H 2 . In the present work, this technique was systematically applied for studying vibrational resonances of practically important isotopologues of the linear carbonyl sulfide molecule ( 16 O 12 C 32 S, 16O12C34S, O 12 C 34 S, 16O13C32S, O 13 C 32 S, 18O12C32S, O 12 C 32 S, 16O12C33S, O 12 C 33 S, 16 O 13 C 34 S). The isomorphic Hamiltonians were constructed using the ab initio equilibrium geometry and quartic PES calculated at the CCSD(T)/cc-pV(Q+d)Z level. The analysis of HPA common branch points of 125 vibrational states for each isotopologue predicted comprehensive resonance pictures. These resonances reproduced most of experimentally observed couplings and indicated a possible break down of the known polyad formula P = 2 v 1 + v2 2 + 4 v 3 . The demonstrated efficiency of this purely ab initio approach opens a perspective of further studies of resonances phenomena of new and hardly accessible molecules.

We demonstrate that machine learning models trained on a set of features obtained from QM/MM molecular dynamic trajectories of fluorescent proteins can be used to predict the chromophore dipole moment variation upon excitation, the quantity related to the electronic excitation energy. Linear regression, gradient boosting, and artificial neural network- based models were considered using cross-validation on the training dataset. Gradient boosting approach proved to be the most accurate for both internal (R2 = 0.77) and external (R2 = 0.7) test sets.

Molecular dynamic simulations using QM/MM potentials are performed for the enzyme-substrate complex of adenosine triphosphate (ATP) with the motor protein myosin. Machine learning methods are applied to a dataset consisting of the geometry parameters of the active site in the enzyme-substrate complex to predict the Laplacian of electron density at the bond critical point of the PG-O3B bond being broken in ATP. Using a gradient boosting machine learning model, a mean absolute error of 0.01 a.u. and an R 2 score of 0.99 are achieved, and it is found that the PG-O3B bond length is the most important feature, contributing 2/3, while other geometry features contribute 1/3.

DNA aptamers are oligonucleotides that specifically bind to target molecules, similar to how antibodies bind to antigens. We identified an aptamer named MEZ that is highly specific to the receptor-binding domain, RBD, of the SARS-CoV-2 spike protein from the Wuhan-Hu-1 strain. The SELEX procedure was utilized to enrich the initial 31-mer oligonucleotide library with the target aptamer. The aptamer identification was performed using the novel protocol based on nanopore sequencing developed in this study. The MEZ aptamer was chemically synthesized and tested for binding with the SARS-CoV-2 RBD of the spike protein from different strains. The Kd is 6.5 nM for the complex with the RBD from the Wuhan-Hu-1 strain, which is comparable with known aptamers; the advantage is that the MEZ aptamer is smaller than known analogs. The proposed aptamer is highly selective for the RBD protein from the Wuhan-Hu-1 strain and does not form complexes with the RBD from Beta, Delta and Omicron strains. Experimental and theoretical studies together revealed the molecular mechanism of aptamer binding. The aptamer occupies the same binding site as ACE2 when bound to the RBD. The 3 '-end of the MEZ aptamer is important for complex formation and is responsible for the discrimination of the RBD protein from a specific strain. The 5 '-end is responsible for the formation of a loop in the 3D structure of the aptamer, which is important for proper binding. MEZ is a 31-mer aptamer that is highly specific to the RBD from the SARS-CoV-2 Wuhan-Hu-1 strain with Kd = 6.5 nM.

The search for efficient inhibitors of the SARS-CoV-2 enzymes is ongoing due to the continuing COVID-19 pandemic. We report the results of computational modeling of the reactions of the SARS-CoV-2 main protease (MPro ) with four potential covalent inhibitors. Two of them, carmofur and nirmatrelvir, have been shown experimentally the ability to inhibit MPro . Two other compounds, X77A and X77C, were designed computationally in this work, derived from the structure of X77, a non-covalent inhibitor forming a tight surface complex with MPro . We modified the X77 structure by introducing warheads capable of efficient chemical reactions with the catalytic cysteine residue in the M Pro active site. The reaction mechanisms of the four molecules with M Pro were investigated by quantum mechanics/molecular mechanics (QM/MM) calculations using large quantum subsystems. First, at the QM/MM level, we optimized structures of stationary points on the potential energy surfaces corresponding to the reactants, products, intermediates, and transition states along the hypothesized reaction coordinates. Analysis of these structures has informed the selection of collective variables for the subsequent calculations of the Gibbs energy profiles using molecular dynamics simulations with QM/MM potentials (QM/MM MD). In these simulations, the QM part was treated by DFT with the PBE0 functional. The results show that all four compounds form covalent adducts with the catalytic cysteine Cys 145 of MPro . From the chemical perspective, the reactions of these four compounds with M Pro follow three distinct mechanisms. In all cases, the reaction is initiated by a nucleophilic attack of the thiolate group of the deprotonated cysteine residue from the catalytic dyad Cys145-His41 of MPro . In the case of carmofur and X77A, the covalent binding of the thiolate to the ligand is accompanied by the formation of the fluoro-uracil leaving group. The reaction with X77C follows the nucleophilic aromatic substitution SN Ar mechanism. The reaction of M Pro with nirmatrelvir, which has a reactive nitrile group, leads to the formation of the covalent thioimidate adduct with the thiolate of the Cys145 residue in the enzyme active site.

Modern quantum-based methods are employed to model interaction of the flavin-dependent enzyme RutA with the uracil and oxygen molecules. This complex presents the structure of reactants for the chain of chemical reactions of monooxygenation in the enzyme active site, which is important in drug metabolism. In this case, application of quantum-based approaches is an essential issue, unlike conventional modeling of protein-ligand interaction with force fields using molecular mechanics and classical molecular dynamics methods. We focus on two difficult problems to characterize the structure of reactants in the RutA-FMN-O2 -uracil complex, where FMN stands for the flavin mononucleotide species. First, location of a small O2 molecule in the triplet spin state in the protein cavities is required. Second, positions of both ligands, O2 and uracil, must be specified in the active site with a comparable accuracy. We show that the methods of molecular dynamics with the interaction potentials of quantum mechanics/molecular mechanics theory (QM/MM MD) allow us to characterize this complex and, in addition, to surmise possible reaction mechanism of uracil oxygenation by RutA.

We report the results of a computational study of the mechanism of the light-induced chemical reaction of chromophore hydration in the fluorescent protein Dreiklang, responsible for its switching from the fluorescent ON-state to the dark OFF-state. We explore the relief of the charge-transfer excited-state potential energy surface in the ON-state to locate minimum energy conical intersection points with the ground-state energy surface. Simulations of the further evolution of model systems allow us to characterize the ground-state reaction intermediate tentatively suggested in the femtosecond studies of the light-induced dynamics in Dreiklang and finally to arrive at the reaction product. The obtained results clarify the details of the photoswitching mechanism in Dreiklang, which is governed by the chemical modification of its chromophore.

Interaction of molecular oxygen 3O2 with the flavin-dependent protein miniSOG after light illumination results in creation of singlet oxygen 1O2 and superoxide O2●−. Despite the recently resolved crystal structures of miniSOG variants, oxygen-binding sites near the flavin chromophore are poorly characterized. We report the results of computational studies of the protein−oxygen systems using molecular dynamics (MD) simulations with force-field interaction potentials and quantum mechanics/molecular mechanics (QM/MM) potentials for the original miniSOG and the mutated protein. We found several oxygen-binding pockets and pointed out possible tunnels bridging the bulk solvent and the isoalloxazine ring of the chromophore. These findings provide an essential step toward understanding photophysical properties of miniSOG—an important singlet oxygen photosensitizer.

We report the results of computational studies of the guanosine triphosphate (GTP) hydrolysis in the active site of the KRas-NF1 protein complex, where KRas stands for the K-isoform of the Ras (ras sarcoma) protein and NF1 (neurofbromin-1) is the activating protein. The model system was constructed using coordinates of heavy atoms from the crystal structure PDB ID 6OB2 with the GTP analog GMPPNP. Large-scale classical molecular dynamics (MD) calculations were performed to analyze conformations of the enzyme-substrate complexes. The Gibbs energy profiles for the hydrolysis reaction were computed using MD simulations with quantum mechanics/molecular mechanics (QM/MM) interaction potentials. The density functional theory DFT(ωB97X-D3/6-31G**) approach was applied in QM and the CHARMM36 force field parameters in MM. The most likely scenario of the chemical step of the GTP hydrolysis in KRas-NF1 corresponds to the water-assisted mechanism of the formation of the inorganic phosphate coupled with the dissociation of GTP to GDP.

The dynamic properties of enzyme-substrate complexes in the reaction of the hydrolysis of guanosine triphosphate (GTP) are compared by the native Ras enzyme and its variant G12VRas with oncogenic point substitution using molecular dynamic methods with potentials constructed according to the theory of quantum mechanics/molecular mechanics (QM/MM). Molecular models of systems including the GTP-binding protein Ras, the GAP accelerator protein, the GTP substrate, and the catalytic water molecule are constructed based on the atomic coordinates of the complex from the database of protein structures. An analysis of the obtained distributions along molecular dynamics trajectories for the distances between the oxygen atom of the catalytic water molecule in the active center of the complex and the phosphorus atom of the γ-phosphate group of GTP showed that the G12V substitution leads to the growth of populations with large distances between the reactants, thus reducing the proportion of the fraction of reactive conformations for the hydrolysis reaction.

The mechanism of the reactivation reaction of the double mutant butyrylcholinesterase (BChE) double mutantAsn322Glu/Glu325Gly, inhibited by the organophosphorus compound (OPC) echothiophat is studied by molecular modeling methods. The ability of this mutant to spontaneously reactivate itself was previously shown experimentally. The energy profile of the reaction path calculated by the method of quantum mechanics/molecular mechanics (QM/MM) confirms that such a mechanism is possible. Molecular dynamics calculations with QM/MM potentials are used to investigate the proton transfer pathways and determine the protonated state of the glutamic acids in the active site of the double mutant. Possible ways of increasing its conformational stability are investigated by the methods of classical molecular dynamics.

This work explores the level of transparency in reporting the details of computational protocols that is required for practical reproducibility of quantum mechanics/molecular mechanics (QM/MM) simulations. Using the reaction of an essential SARS-CoV-2 enzyme (the main protease) with a covalent inhibitor (carmofur) as a test case of chemical reactions in biomolecules, we carried out QM/MM calculations to determine the structures and energies of the reactants, the product, and the transition state/intermediate using analogous QM/MM models implemented in two software packages, NWChem and Q-Chem. Our main benchmarking goal was to reproduce the key energetics computed with the two packages. Our results indicate that quantitative agreement (within the numerical thresholds used in calculations) is difficult to achieve. We show that rather minor details of QM/ MM simulations must be reported in order to ensure the reproducibility of the results and offer suggestions toward developing practical guidelines for reporting the results of biosimulations.

Supercomputer molecular modeling methods are applied to characterize structure and dynamics of the flavin-dependent enzyme RutA in the complex with molecular oxygen. Following construction of a model protein system, molecular dynamics (MD) simulations were carried out using either classical force field interaction potentials or the quantum mechanics/molecular mechanics (QM/MM) potentials. Several oxygen-binding pockets in the protein cavities were located in these simulations. The QM/MM-based MD calculations rely on the interface between the quantum chemistry package TeraChem and the MD package NAMD. The results show a stable localization of the oxygen molecule in the enzyme active site. Static QM/MM calculations carried out with two different packages, NWChem and TURBOMOLE, allowed us to establish the structure of the RutA-O2 complex. Biochemical perspectives of the hallmark reaction of incorporating oxygen into organic compounds emerged from these simulations are formulated.