Eugen Hruska

Despite widespread use of molecular dynamics (MD) to refine docking poses, no large-scale statistical validation exists for metalloproteins like Cytochrome P450 (CYP450). We present the first statistically rigorous benchmark of docking engines and MD refinement on CYP450 systems, evaluating AutoDock Vina, GNINA, GalaxyDock2-HEME, and RoseTTAFold All-Atom (RFAA) across crossdocking scenarios. Our analyses revealed that short MD simulations provide no statistically significant geometric refinement toward experimental structures. Using automated interaction comparison that categorises protein-ligand contacts relative to crystal references, we found this behaviour showed no significant dependence on ligand size or flexibility. Among docking methods, RFAA achieved geometric drift approximately three times lower than AutoDock Vina. These findings provide clarity regarding the efficacy of MD refinement for organometallic complexes. The automated framework developed for this study (DOCKTOPUS/CYCLOPS packages) enables the high-throughput simulations and statistical testing necessary to rigorously validate computational methods and reduce publication bias inherent in stochastic simulations.

Molecular docking is a standard tool for predicting ligand binding, yet refining these poses via molecular dynamics (MD) remains challenging to automate for metalloproteins like Cytochrome P450 (CYP450). Here, we introduce DOCKTOPUS (Dock to P450) and CYCLOPS (CYtochrome Complex Ligand Optimization with Protein Simulation) Python packages designed to automate the docking and molecular dynamics setup for CYP450-ligand complexes. These tools integrate seamlessly with representative predictive engines AutoDock Vina, GNINA, GalaxyDock2-HEME, and RoseTTAFold All-Atom (RFAA) by automatically handling the complex parameterisation of the covalently bound heme-thiolate cofactor and ligand topologies, while also providing capabilities to filter out unphysical ligand geometries generated by RFAA. We validated the pipeline through a medium-sized case study to enable a rigorous statistical analysis. In crossdocking scenarios RFAA achieved a mean geometric drift relative to experiments approximately three times lower than that of AutoDock Vina. Our structural and pharmacophoric analyses revealed that subsequent MD simulations provide no statistically significant geometric refinement towards the experimental crystal structure. This behaviour was remarkably consistent: we observed no significant differences by specific interaction type or ligand flexibility. These findings suggest that short MD refinement with contemporary classical force fields for organometallic systems might have limited effects beyond pose validity. Furthermore the DOCKTOPUS and CYCLOPS packages offer a robust platform for statistical testing of docking and molecular dynamics effects.

Docking and molecular dynamics are commonly used methods in drug development for the prediction of interactions between the drug in question and the target protein. When used together, docking studies can provide starting structures for molecular dynamics, which can improve the relatively low accuracy of docking but at greater computational cost. When the structure contains a cofactor with covalently bound metal atom, few automated methods are capable of setting up the system to perform molecular dynamics simulations. Another limitation that can be encountered is the parametrization of the ligand, which is not always easy to set up and automate using popular suites, such as AMBER or GROMACS. Here, we propose an automated workflow for setting up and running the simulation, providing a seamless integration for both flexible and rigid protein docking engines. The automatic setup works for around 90% of all ligands for rigid docking engines, with this number decreasing for RoseTTAFold All-Atom (flexible docking engine) to around 30% due to its known problem with hallucination of the ligand conformation.

In this work, we benchmark 4 selected open source docking engines for use in the cytochrome P450 protein family. The key enzymes family of phase I metabolism is characterized by a wide variety of accepted substrates due to flexible active site. This work is a benchmark study which aims to evaluate the capabilities of current rigid and induced-fit docking methods for prediction of correct heme-ligand orientation. To asses it, we use two unique distances to heme iron and a SuCOS score to quantify reconstruction of orientation and chemical features. We selected three rigid protein docking engines: GNINA, AutoDock VINA, GalaxyDock2 HEME and a flexible docking model, RosettaFold-All-Atoms to test them on a dataset of 128 CYP-binding ligands. We report mean absolute error for RosetttaFold-All-Atom on key distance, to the atom closest to heme iron in experimental reference structure, 3 times lower than AutoDock VINA engine in the same simulation. Our results indicate that induced fit method is a significant improvement over rigid methods for flexible active site, but still offer limited predictivity. During crossdocking, RosettaFold-All-Atoms was able to recreate over a quarter of distances up to 20 percent difference from experiment. Further analysis indicates a low overlap in the distribution of ligand chemical features, based on a SuCOS score, which suggests a space for further improvement.

We investigated the various factors impacting the performance of Δ-machine learning (Δ-ML) solution phase molecular properties.

Due to the limitation of solvent models, quantum chemistry calculated solution-phase molecular properties often deviates from experimental measurements. Recently, ∆-machine learning (∆-ML) was shown to be a promising approach to correcting errors in the quantum chemistry calculation of solvated molecules. However, this approach's applicability to different molecular properties and its performance in various use cases are still unknown. In this work, we tested the performance of Δ-ML in correcting redox potential and absorption energy calculations using four types of input descriptors and various ML methods. We sought to understand the dependence of ∆-ML performance on the property to predict, the quantum chemistry method, the data set distribution/size, the type of input features, and the feature selection techniques. We found that ∆-ML can effectively correct the errors in redox potentials calculated by density functional theory (DFT) and absorption energies calculated by time-dependent DFT. For both properties, the ∆-ML corrected results showed less sensitivity to the DFT functional choice than the raw results. The optimal input descriptor depends on the property, regardless of the specific ML method used. The solvent-solute descriptor (SS) is the best for redox potential, whereas the combined molecular fingerprint (cFP) is the best for absorption energy. A detailed analysis of the feature space and the physical foundation of different descriptors well explained these observations. Feature selection did not further improve the Δ-ML performance. Finally, we analyzed the limitation of our Δ-ML solvent effects approach in data sets with molecules of varying degrees of electronic structure errors.

Interpretation of chemistry on an atomic scale improves with explainable artificial intelligence (XAI). The parts of the molecule with the most significant influence on the chemical property of interest can be visualized with atomwise and bondwise attributions. Nonetheless, the attributions from different XAI methods regularly disagree substantially, causing uncertainty about which explanation is correct. To determine a ground truth for attributions, we define chemical operations which avoid alchemical steps or approximations and allow extracting one attribution per atom or bond from existing datasets of chemical properties. This general procedure allows for generating large datasets of ground truth attributions. The approach allowed us to create a ground truth explanation dataset with more than 5 million data points for the HOMO-LUMO gap chemical property. This open-source dataset of atomistic ground truth explanations may serve as a reference for XAI approaches.

The availability of large, high-quality datasets is crucial for artificial intelligence design and discovery in chemistry. Despite the essential roles of solvents in chemistry, the rapid computational dataset generation of solution-phase molecular properties at the quantum mechanical level of theory was previously hampered by the complicated simulation procedure. Software toolkits that can automate the procedure to set up high-throughput explicit-solvent quantum chemistry (QC) calculations for arbitrary solutes and solvents in an open-source framework are still lacking. We developed AutoSolvate, an open-source toolkit, to streamline the workflow for QC calculation of explicitly solvated molecules. It automates the solvated-structure generation, force field fitting, configuration sampling, and the final extraction of microsolvated cluster structures that QC packages can readily use to predict molecular properties of interest. AutoSolvate is available through both a command line interface and a graphical user interface, making it accessible to the broader scientific community. To improve the quality of the initial structures generated by AutoSolvate, we investigated the dependence of solute–solvent closeness on solute/solvent identities and trained a machine learning model to predict the closeness and guide initial structure generation. Finally, we tested the capability of AutoSolvate for rapid dataset curation by calculating the outer-sphere reorganization energy of a large dataset of 166 redox couples, which demonstrated the promise of the AutoSolvate package for chemical discovery efforts.

Prediction of redox potentials is essential for catalysis and energy storage. Although density functional theory (DFT) calculations have enabled rapid redox potential predictions for numerous compounds, prominent errors persist compared to experimental measurements. In this work, we develop machine learning (ML) models to reduce the errors of redox potential calculations in both implicit and explicit solvent models. Training and testing of the ML correction models are based on the diverse ROP313 data set with experimental redox potentials measured for organic and organometallic compounds in a variety of solvents. For the implicit solvent approach, our ML models can reduce both the systematic bias and the number of outliers. ML corrected redox potentials also demonstrate less sensitivity to DFT functional choice. For the explicit solvent approach, we significantly reduce the computational costs by embedding the microsolvated cluster in implicit bulk solvent, obtaining converged redox potential results with a smaller solvation shell. This combined implicit–explicit solvent model, together with GPU-accelerated quantum chemistry methods, enabled rapid generation of a large data set of explicit-solvent-calculated redox potentials for 165 organic compounds, allowing detailed investigation of the error sources in explicit solvent redox potential calculations.