A method for predicting hERG cardiotoxicity of compounds based on multidimensional molecular fingerprints

Abstract

A method for predicting compound hERG cardiotoxicity based on multidimensional molecular fingerprints, the method comprising: the first step, performing molecular dynamics simulation on compounds in a cardiotoxicity data set; the second step, analyzing the molecular dynamics simulation results and calculating molecular fingerprints based on molecular dynamics simulation; the third step, calculating the molecular fingerprint of the compound by molecular fingerprint software, and forming a multidimensional molecular fingerprint with the molecular dynamics simulation fingerprint; the fourth step, using feature selection methods to find the optimal feature subset from the multidimensional molecular fingerprint; the fifth step, using multiple machine learning methods to construct a compound hERG cardiotoxicity prediction model. Experimental results show that the model established using multidimensional molecular fingerprints has higher prediction performance. This method can achieve higher compound hERG cardiotoxicity prediction accuracy and can be used for early toxicity screening in drug design.

Description

Technical Field

[0001] The present invention relates to the field of computer-aided drug screening, and in particular to a method for predicting the hERG cardiotoxicity of a compound based on multidimensional molecular fingerprints. Background Art

[0002] Currently, drug toxicity has become a serious medical and economic burden. acLQTS (acquired long QT syndrome) is an acquired cardiac ion channel disease caused by drugs blocking the hERG channel. Acquired long QT syndrome (acLQTS) is one of the most important diseases caused by drug-induced toxicity. It is a potentially life-threatening arrhythmia disease that also increases the risk of syncope, sudden cardiac death, and epileptic seizures.

[0003] The hERG protein is a tetrameric potassium ion channel that is primarily associated with cardiotoxicity and acLQTS. Therefore, it is essential to avoid drug cardiotoxicity during drug screening, and the use of computer models can quickly predict the ability of a compound to bind to hERG and assess the compound's cardiotoxicity.

[0004] Molecular fingerprint is a molecular structure encoding method. It can generally calculate the physical and chemical properties, two-dimensional structural characteristics and three-dimensional structural characteristics of compounds, and can be used to extract the characteristics of compounds in compound property prediction models. Molecular dynamics fingerprint (MDFP) is a fingerprint based on the trajectory of computational molecular dynamics simulation. After adding the time dimension, MDFP can reflect the characteristics of molecules in the time dimension. Studies have found that in the study of free energy prediction of heterogeneous fusion models constructed using MDFP, MDFP showed good performance. Therefore, MDFP can be used as an alternative to traditional molecular fingerprints and has great application potential in the prediction of hERG cardiotoxicity. Summary of the invention

[0005] The present invention provides a method for predicting the hERG cardiotoxicity of a compound based on multidimensional molecular fingerprints.

[0006] In order to achieve the above purpose, the present invention adopts the following technical solutions:

[0007] Step 1: Perform molecular dynamics simulation on the compounds in the cardiotoxicity dataset. The specific process is as follows:

[0008] First, generate the three-dimensional structure file of the compound molecule, and then use the AMBER all-atom force field to generate the topological structure of the compound. Then, place the compound in a water box with a diameter of 1.0-3.0nm, minimize the energy of the system, and then balance the system under NVT and NPT systems respectively. Finally, perform molecular dynamics simulation on the system and collect the simulation trajectory.

[0009] Step 2: Analyze the molecular dynamics simulation results obtained in step 1 and calculate the molecular fingerprint based on the molecular dynamics simulation. The specific process is as follows:

[0010] The molecular dynamics simulation trajectory was analyzed by the GROMACS toolkit, and the gyration radius rg, root mean square deviation rmsd, solvent accessible surface area sasa, total energy tenergy, number of hydrogen bonds hbond, and kinetic energy kinetic of the small molecule in the simulation were calculated as the four-dimensional molecular fingerprint of the compound.

[0011] Step 3: Calculate the traditional molecular fingerprint of the compound through molecular fingerprint software, and combine it with the molecular dynamics simulation fingerprint to form a multi-dimensional molecular fingerprint;

[0012] The traditional molecular fingerprints include one-dimensional, two-dimensional and three-dimensional molecular fingerprints of compounds such as PropertyFP, Baseline2D and ECFP.

[0013] Step 4: Use feature selection method to find the optimal feature subset from the multidimensional molecular fingerprint;

[0014] Firstly, the zero-variation and near-zero-variation features in all molecular fingerprints were deleted, and then recursive feature elimination (RFE) was performed to select the optimal feature subset.

[0015] Step 5: Use multiple machine learning methods to construct a compound hERG cardiotoxicity prediction model.

[0016] Machine learning algorithms such as random forest, support vector machine, gradient boosting, and partial least squares regression were used to construct a prediction model for hERG cardiotoxicity with the optimal feature subset as training data.

[0017] Step 6: For the constructed compound hERG cardiotoxicity prediction model, the model was evaluated using the cross-validation method, and the root mean square error (RMSE), mean unsigned error (MUE) and R2 were used as indicators to evaluate the prediction performance of the model.

[0018] Beneficial effects of the present invention: The present invention adopts the above-mentioned scheme, and the four-dimensional molecular fingerprint based on molecular dynamics simulation is used to overcome the limitations of previous methods, and the obtained prediction results are more accurate. The present invention retrieved 203 compounds that have the ability to inhibit hERG protein activity from the literature, and used them as training data sets, and established a prediction model using four machine learning algorithms based on multidimensional molecular fingerprints. Experimental results show that this model can achieve higher prediction accuracy when combined with traditional molecular fingerprints. The present invention proposes a four-dimensional molecular fingerprint (MDFP) based on molecular dynamics simulation and uses it to predict the cardiotoxicity of compounds. For the first time, molecular dynamics simulation fingerprints and other dimensional molecular fingerprints are combined to improve the hERG cardiotoxicity prediction accuracy of the model, and an integrated method is used to integrate multiple machine learning models to make the prediction results better than a single model. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 The figure is a flow chart of the method of the present invention.

[0020] Figure 2 This is a schematic diagram of MDFP variable properties.

[0021] Including kinetic, LJ-14, LJ-SR, tenergy, rgyr, hbond, sasa, rmsd. Each attribute is expressed by avr (average), std (standard deviation) and mid (median).

[0022] Figure 3a Schematic diagram of selecting important features from MDFP fingerprints using recursive feature elimination.

[0023] Figure 3b Schematic diagram of selecting important features from Baseline2D fingerprints using recursive feature elimination.

[0024] Figure 3c Schematic diagram of selecting important features from ECFP4 fingerprints using recursive feature elimination.

[0025] Figure 3d Schematic diagram of selecting important features from PropertyFP fingerprints using recursive feature elimination.

[0026] Figure 4a Schematic diagram of the correlation between the integrated predicted values ​​and experimental values ​​of each model.

[0027] Figure 4b Schematic diagram of the correlation between GBM's predicted values ​​and experimental values.

[0028] Figure 4c Schematic diagram of the correlation between the predicted values ​​of PLS ​​and the experimental values.

[0029] Figure 4d Schematic diagram of the correlation between the predicted and experimental values ​​of RF.

[0030] Figure 4e Schematic diagram of the correlation between the SVM predicted value and the experimental value. DETAILED DESCRIPTION

[0031] The technical solutions in the embodiments of the present invention will be described clearly and completely below. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all of the embodiments.

[0032] Embodiment 1:

[0033] The compound hERG cardiotoxicity prediction method based on multidimensional molecular fingerprints comprises the following steps:

[0034] The first step is data preprocessing, collecting high-quality hERG inhibitor data sets from the literature. The half-inhibitory concentration value (IC50) represents the biochemical half-maximum inhibition concentration and is used to indicate the inhibitory ability of the compounds in this data set on hERG. Delete the toxicity data with the same name and half-inhibitory concentration value. If the name of the molecule is repeated and the difference in IC50 value is less than one order of magnitude, the IC50 value of this molecule is averaged. Finally, 203 molecules with specific hERG half-inhibitory concentration values ​​are obtained. The obtained data are grouped, with 80% for training set and 20% for validation set. The training set uses five-fold cross-validation, and the validation set is used to test the predictive performance of the new compound model. And the IC50 value is converted to pIC50 (negative logarithm of IC50 value) value.

[0035] In the second step, molecular dynamics simulation of the compound was performed based on GROMACS. The 3D structure file of the molecule was obtained through the SMILES file in Zinc15. The topological structure of the compound was generated using the AMBER14 force field in ACPYPE. The compound was placed in a dodecahedral box with a center size of 1.0 nm and solvated using the TIP3P water model. Then, the system was energy minimized for 100 ps and balanced for 1 ns under the NVT and NPT (300 K, 1 bar) ensembles, respectively. Finally, the system ran a 5 ns molecular dynamics simulation, writing coordinates every 10 ps and energies every 1 ps.

[0036] The third step is to obtain three 2D molecular fingerprints. Use RDKit to obtain the Baseline2D molecular fingerprint. Use the Descriptastorus package of RDKit to obtain the PropertyFP molecular fingerprint. Use the Morgan algorithm of RDKit (verctor length is 2048 and radius is 2) to obtain the ECFP4 molecular fingerprint.

[0037] The fourth step is to calculate the molecular fingerprint based on molecular dynamics simulation. The molecular dynamics trajectory is obtained by analyzing the GROMACS toolkit, and the gyration radius (rg), root mean square deviation (rmsd), solvent accessible surface area (sasa), total energy (tenergy), number of hydrogen bonds (hbond), kinetic energy (kinetic) LJ-SR, LJ14 feature information is used as the four-dimensional molecular fingerprint of the compound. Then the R language package is used to obtain the mean, median and standard deviation of the features.

[0038] The fifth step is feature selection. First, the nearZeroVar function of the caret package in R language is used to remove zero-variation and near-zero-variation features in all molecular fingerprints, and then the rfe function (ten five-fold cross-validation settings) is used to perform recursive feature elimination (RFE) to select the optimal feature subset. In the RFE process, all features are first ranked according to the feature importance values ​​obtained by the random forest (RF) algorithm, and then the gradually reduced features are iteratively trained according to the ranking to evaluate the performance of the feature subset.

[0039] The sixth step is to construct a machine learning model. We use machine learning algorithms such as random forest, support vector machine, gradient boosting, partial least squares regression, etc. to build a prediction model for hERG cardiotoxicity with the optimal feature subset as training data:

[0040] Construction of RF: In this study, the random forest function was used to construct the random forest classifier. The number of randomly selected classification trees and variables for each node overflow was set to ntree=500, and mtry was optimized from one-third of the number of features minus 10 to plus 15.

[0041] Construction of SVM: This study uses radial basis function (RBF) kernel to construct SVM classifier. At the same time, a random search method is used to optimize specific SVM parameters using regularization parameter C and σ parameter. C is calculated from e -2 To e 6 , σ from e -7 to e, with a step length of e 0.5 .

[0042] GBM construction: In this study, the total number of trees (n.tree) and the maximum depth of each tree (interaction.depth) were optimized using the caret package from 1 to 3000 and 1 to 10, respectively. shrinkage and n.minobsinnode were set to 0.005 and 10, respectively.

[0043] Construction of PLS: In this study, n_components of PLS ​​was optimized from 1 to the maximum eigenvalue or sample size.

[0044] Step 6. Model Evaluation

[0045] (1) The model was evaluated using a 10-fold 5-fold cross validation method. The original dataset was randomly divided into 5 equal subsets, 4 of which were used for training and the other for testing. The 5-fold cross validation was then repeated 10 times to reduce randomness. The cross validation was performed 10 times with different random seeds. Then, the average value was calculated to evaluate the prediction performance of the model.

[0046] Use root mean square error (RMSE), mean unsigned error (MUE), and R 2 To evaluate the prediction performance of the model. The calculation formulas for these indicators are as follows:

[0047]

[0048]

[0049]

[0050] Among them, P, E,n represent the predicted value, the average of experimental values, the experimental value, and the number of molecules, respectively.

[0051] The results show that the RMSE of MDFP is 0.967-0.139, second only to 0.860-0.960 of PropertyFP. ECFP4 and Baseline2D are relatively poor. R2 and MUE also show the same trend (see the table below).

[0052]

Claims

1. A method for predicting the hERG cardiotoxicity of a compound based on multidimensional molecular fingerprints, characterized in that The following steps are involved: Step 1: molecular dynamics simulation of the compounds in the cardiotoxicity dataset; Step 2: Analyze the molecular dynamics simulation results obtained in step 1 and calculate the four-dimensional molecular fingerprint based on the molecular dynamics simulation; Step 3: Calculate the traditional one-dimensional, two-dimensional and three-dimensional molecular fingerprints of the compound through molecular fingerprint software, and combine them with the molecular dynamics simulation fingerprint to form a multi-dimensional molecular fingerprint; Step 4: Use feature selection method to find the optimal feature subset from the multidimensional molecular fingerprint; Step 5: Use multiple machine learning methods to build a compound hERG cardiotoxicity prediction model; The specific process of step 1 is as follows: first, generate a three-dimensional structure file of the compound molecule, then use the AMBER all-atom force field to generate the topological structure of the compound, then place the compound in a water box with a diameter of 1.0-3.0 nm, minimize the energy of the system, then balance the system under NVT and NPT ensembles respectively, finally perform molecular dynamics simulation on the system, and collect simulation trajectories; The specific process of step 2 is as follows: using the GROMACS toolkit to analyze the molecular dynamics simulation trajectory, calculate the radius of gyration rg, root mean square deviation rmsd, solvent accessible surface area sasa, total energy tenergy, number of hydrogen bonds hbond, and kinetic energy kinetic of the small molecule in the simulation as the four-dimensional molecular fingerprint of the compound; The traditional molecular fingerprints described in step 3 include one-dimensional, two-dimensional, and three-dimensional molecular fingerprints of PropertyFP, Baseline2D, and ECFP compounds; The specific process of step 4 is as follows: firstly, the zero-variation and near-zero-variation features in all molecular fingerprints are deleted, and then recursive feature elimination (RFE) is performed to select the optimal feature subset; The multiple machine learning methods in step five include: random forest, support vector machine, gradient boosting, partial least squares regression machine learning algorithms, and the optimal feature subset is used as training data to construct a prediction model for compound hERG cardiotoxicity.

2. The method for predicting the hERG cardiotoxicity of a compound based on multidimensional molecular fingerprint according to claim 1, characterized in that: The sixth step is to evaluate the constructed compound hERG cardiotoxicity prediction model using a cross-validation method, using the root mean square error RMSE, mean unsigned error MUE and R 2 The predictive performance of the model is evaluated as an indicator.

Images