Methods, systems, and media of quantum mechanics-based drug long-time course kinetics simulation

Abstract

This invention provides a method, system, and medium for long-term drug dynamics simulation based on quantum mechanics. The method includes: acquiring and preprocessing target protein crystal structure coordinate data and drug molecule structure data; obtaining a stabilized initial atomic coordinate and topological parameter dataset and performing atomic-level partitioning; then obtaining an enhanced sampling dynamics-driven dataset for long-term simulation; obtaining a quality-controlled, qualified long-term atomic trajectory time-series dataset and processing it in conjunction with the enhanced sampling dynamics-driven dataset; finally obtaining a long-term drug dynamics analytical dataset. Through the construction and data preprocessing of the initial complex system, QM / MM partitioning and quantum mechanical force data generation, enhanced sampling and collective variable data-driven simulation, long-term trajectory data acquisition and quality control, and free energy surface reconstruction and conformational data extraction, high-precision sampling of rare conformational events of protein-drug complexes is achieved, improving the accuracy and completeness of drug action mechanism analysis.

Description

Technical Field

[0001] This invention relates to the fields of computational drug design and computational biophysics, and more specifically, to a method, system, and medium for long-term drug dynamics simulation based on quantum mechanics. Background Technology

[0002] Traditional molecular dynamics simulations are limited to nanosecond to submicrosecond timescales, making it difficult to capture key biophysical processes at the microsecond to millisecond scale, such as protein allosteric regulation, binding pocket reconstruction, and drug-induced conformational changes. Conventional classical molecular dynamics cannot describe electronic effects such as covalent bonding, proton transfer, and charge redistribution. Neither enhanced sampling nor QM / MM alone can simultaneously achieve high precision, long-term simulation, and efficient sampling. Existing technologies, published literature and patents only disclose separate classical molecular dynamics long-term simulations, separate QM / MM calculations, and separate meta-dynamics enhanced sampling. No publicly available schemes have been found that deeply couple real-time QM / MM, meta-dynamics, and tens of microsecond long-term simulations to form a complete data processing chain. This results in a trade-off between simulation accuracy, sampling efficiency, and timescale, failing to meet the demands of modern drug development for atomic-level dynamic mechanism analysis. Summary of the Invention

[0003] The purpose of this invention is to provide a method, system, and medium for long-term drug dynamics simulation based on quantum mechanics, which can achieve efficient and high-precision sampling of rare conformational events of protein-drug complexes, thereby improving the accuracy and completeness of drug action mechanism analysis.

[0004] This invention also provides a long-term drug dynamics simulation method based on quantum mechanics, comprising the following steps: The target protein crystal structure coordinate file data and drug molecule structure data are acquired and preprocessed to obtain initial complex conformation data and construct a periodic simulation box. Then, the coordinate data of the complete solvation and ion addition system are obtained and hierarchical energy minimization processing is performed to obtain a stable initial atomic coordinate and topological parameter dataset. Atomic partitioning is performed based on the stabilized initial atomic coordinates and topology parameter dataset. Then, the partition index dataset is obtained and calculated and coupled. Finally, the total energy of QM / MM coupling and the full atomic coupling force dataset are obtained and encapsulated to obtain the frame-by-frame force source dataset. Based on the frame-by-frame force source dataset, normalization and scale calibration are performed to obtain the spatial features of collective variables and the adaptive parameter dataset. An adaptive Gaussian bias potential is then applied, and the enhanced sampling dynamics driving dataset is obtained. Long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain a time-series trajectory dataset and perform quality control monitoring, thereby obtaining a quality-controlled qualified long-term atomic trajectory time-series dataset. The data is processed based on the quality control qualified long-term atomic trajectory time series dataset and the enhanced sampling dynamics driven dataset to obtain the free energy extremum and steady-state feature dataset. Principal component analysis and cluster analysis are then performed to obtain the transition path and transition state dataset and perform interaction time series analysis, ultimately obtaining the drug long-term dynamics analysis dataset.

[0005] Optionally, the step of acquiring target protein crystal structure coordinate file data and drug molecule structure data and preprocessing them to obtain initial complex conformation data and construct a periodic simulation box, thereby obtaining coordinate data of the complete solvation and ionization system and performing hierarchical energy minimization processing to obtain a stable initial atomic coordinate and topological parameter dataset, includes: The target protein crystal structure coordinate file data is obtained and preprocessed to obtain a pure atomic coordinate dataset of the protein. Acquire drug molecular structure data and optimize it to obtain a low-energy three-dimensional conformation; Preprocessing is performed based on the low-energy three-dimensional conformation to obtain a three-dimensional coordinate dataset of drug ligands; Molecular docking was performed based on the protein pure atomic coordinate dataset and the drug ligand three-dimensional coordinate dataset to obtain initial complex conformation data. A periodic simulation box is constructed based on the initial complex conformation data to obtain an all-atom solvation coordinate dataset; Based on the aforementioned all-atom solvation coordinate dataset, the coordinate data of the complete solvation and ion addition system are obtained through calculation. Based on the coordinate data of the complete solvation-ion addition system, a hierarchical energy minimization process is performed to obtain a stable initial atomic coordinate and topological parameter dataset.

[0006] Optionally, the step of performing atomic-level partitioning based on the stabilized initial atomic coordinates and topology parameter dataset, then obtaining a partition index dataset and performing calculations and coupling, and then obtaining the QM / MM coupled total energy and the all-atomic coupled force dataset and encapsulating it, to obtain a frame-by-frame force source dataset, includes: Atomic partitioning is performed based on the stabilized initial atomic coordinates and topological parameter dataset to obtain a partition index dataset; The partitioned index dataset includes the QM atomic index table, the MM atomic index table, and the QM / MM boundary atomic index table; The boundary correction parameter dataset is obtained by calculating the QM / MM boundary atom index table using a preset boundary processing model. The QM atom index table is calculated using a pre-defined theoretical model to obtain QM coupling data, including QM energy, QM atom forces, and QM region electron density and charge datasets. The MM atom index table is calculated using a preset molecular force field model to obtain MM coupling data, including MM energy, MM atom force and non-bonded interaction parameter datasets. The QM coupling data, MM coupling data and boundary correction parameter dataset are coupled by a preset coupling model to obtain the total QM / MM coupling energy and the all-atom coupling force dataset. The QM / MM coupled total energy and the all-atom coupled force dataset are encapsulated to obtain a frame-by-frame force source dataset, which includes atom number, coordinates, force vector, potential energy components and boundary correction terms.

[0007] Optionally, the frame-by-frame force source dataset is normalized and scaled to obtain the spatial features of the collective variables and the adaptive parameter dataset, and an adaptive Gaussian bias potential is applied to obtain the enhanced sampling dynamics-driven dataset, including: Based on the frame-by-frame force source dataset, normalization and scale calibration are performed to obtain collective variable definitions and physical characterization data; Pre-simulation is performed based on the definition of collective variables and physical representation data to obtain the spatial characteristics of collective variables and adaptive parameter dataset; An adaptive Gaussian bias potential is applied based on the spatial characteristics of the collective variables and the adaptive parameter dataset to obtain the bias potential; The bias force and free energy increment are calculated based on the bias potential and collective variables. The enhanced sampling dynamics-driven dataset is obtained by superimposing the bias force and the frame-by-frame force source dataset.

[0008] Optionally, long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain a time-series trajectory dataset and perform quality control monitoring, thereby obtaining a quality-controlled, qualified long-term atomic trajectory time-series dataset, including: Long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain dynamic process parameters. Based on the dynamic process parameters, a time-series trajectory dataset is obtained; Quality control monitoring is performed based on the time-series trajectory dataset to obtain a quality-controlled, qualified long-term atomic trajectory time-series dataset.

[0009] Optionally, the quality control qualified long-term atomic trajectory time series dataset is processed in conjunction with the enhanced sampling kinetics-driven dataset to obtain a free energy extremum and steady-state feature dataset, and principal component analysis and cluster analysis are performed to obtain a transition path and transition state dataset, and interaction time series analysis is performed to finally obtain a drug long-term kinetic analysis dataset, including: Preprocessing is performed on the quality-controlled qualified long-term atomic trajectory time series dataset to obtain standardized trajectory and collective variable time series datasets; The dataset is reweighted based on the enhanced sampling dynamics driving dataset to obtain a two-dimensional free energy surface dataset. Based on the standardized trajectory and collective variable time series dataset combined with the two-dimensional free energy surface dataset, local extremum search and stable conformation state identification are performed to obtain free energy extremum and steady-state feature datasets. Principal component analysis and cluster analysis were performed based on the free energy extreme value and steady-state characteristic dataset to obtain conformational clusters and representative structure datasets; Based on the conformational clusters and representative structure datasets, conformational transition paths and transition states are extracted to obtain a dataset of transition paths and transition states. Based on the interaction time series analysis of the transition path and the transition state dataset, a long-term drug kinetic analysis dataset is obtained.

[0010] Secondly, embodiments of the present invention provide a long-term drug dynamics simulation system based on quantum mechanics. This system includes a memory and a processor. The memory includes a program for a long-term drug dynamics simulation method based on quantum mechanics. When the program for the long-term drug dynamics simulation method based on quantum mechanics is executed by the processor, it performs the following steps: The target protein crystal structure coordinate file data and drug molecule structure data are acquired and preprocessed to obtain initial complex conformation data and construct a periodic simulation box. Then, the coordinate data of the complete solvation and ion addition system are obtained and hierarchical energy minimization processing is performed to obtain a stable initial atomic coordinate and topological parameter dataset. Atomic partitioning is performed based on the stabilized initial atomic coordinates and topology parameter dataset. Then, the partition index dataset is obtained and calculated and coupled. Finally, the total energy of QM / MM coupling and the full atomic coupling force dataset are obtained and encapsulated to obtain the frame-by-frame force source dataset. Based on the frame-by-frame force source dataset, normalization and scale calibration are performed to obtain the spatial features of collective variables and the adaptive parameter dataset. An adaptive Gaussian bias potential is then applied, and the enhanced sampling dynamics driving dataset is obtained. Long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain a time-series trajectory dataset and perform quality control monitoring, thereby obtaining a quality-controlled qualified long-term atomic trajectory time-series dataset. The data is processed based on the quality control qualified long-term atomic trajectory time series dataset and the enhanced sampling dynamics driven dataset to obtain the free energy extremum and steady-state feature dataset. Principal component analysis and cluster analysis are then performed to obtain the transition path and transition state dataset and perform interaction time series analysis, ultimately obtaining the drug long-term dynamics analysis dataset.

[0011] Optionally, the step of acquiring target protein crystal structure coordinate file data and drug molecule structure data and preprocessing them to obtain initial complex conformation data and construct a periodic simulation box, thereby obtaining coordinate data of the complete solvation and ionization system and performing hierarchical energy minimization processing to obtain a stable initial atomic coordinate and topological parameter dataset, includes: The target protein crystal structure coordinate file data is obtained and preprocessed to obtain a pure atomic coordinate dataset of the protein. Acquire drug molecular structure data and optimize it to obtain a low-energy three-dimensional conformation; Preprocessing is performed based on the low-energy three-dimensional conformation to obtain a three-dimensional coordinate dataset of drug ligands; Molecular docking was performed based on the protein pure atomic coordinate dataset and the drug ligand three-dimensional coordinate dataset to obtain initial complex conformation data. A periodic simulation box is constructed based on the initial complex conformation data to obtain an all-atom solvation coordinate dataset; Based on the aforementioned all-atom solvation coordinate dataset, the coordinate data of the complete solvation and ion addition system are obtained through calculation. Based on the coordinate data of the complete solvation-ion addition system, a hierarchical energy minimization process is performed to obtain a stable initial atomic coordinate and topological parameter dataset.

[0012] Optionally, the step of performing atomic-level partitioning based on the stabilized initial atomic coordinates and topology parameter dataset, then obtaining a partition index dataset and performing calculations and coupling, and then obtaining the QM / MM coupled total energy and the all-atomic coupled force dataset and encapsulating it, to obtain a frame-by-frame force source dataset, includes: Atomic partitioning is performed based on the stabilized initial atomic coordinates and topological parameter dataset to obtain a partition index dataset; The partitioned index dataset includes the QM atomic index table, the MM atomic index table, and the QM / MM boundary atomic index table; The boundary correction parameter dataset is obtained by calculating the QM / MM boundary atom index table using a preset boundary processing model. The QM atom index table is calculated using a pre-defined theoretical model to obtain QM coupling data, including QM energy, QM atom forces, and QM region electron density and charge datasets. The MM atom index table is calculated using a preset molecular force field model to obtain MM coupling data, including MM energy, MM atom force and non-bonded interaction parameter datasets. The QM coupling data, MM coupling data and boundary correction parameter dataset are coupled by a preset coupling model to obtain the total QM / MM coupling energy and the all-atom coupling force dataset. The QM / MM coupled total energy and the all-atom coupled force dataset are encapsulated to obtain a frame-by-frame force source dataset, which includes atom number, coordinates, force vector, potential energy components and boundary correction terms.

[0013] Thirdly, embodiments of the present invention also provide a computer-readable storage medium, the computer-readable storage medium including a quantum mechanics-based long-term drug kinetic simulation method program, which, when executed by a processor, implements the steps of the quantum mechanics-based long-term drug kinetic simulation method as described in any of the preceding claims.

[0014] According to the technical solution of the present invention, through the construction and data preprocessing of the initial complex system, QM / MM partitioning and quantum mechanical force data generation, enhanced sampling and collective variable data driving, parallel acquisition and quality control of long-term trajectory data, and free energy surface reconstruction and conformation data extraction, efficient and high-precision sampling of rare conformational events of protein and drug complexes is achieved, thereby improving the accuracy and completeness of drug action mechanism analysis.

[0015] Other features and advantages of the invention will be set forth in the following description, and will be apparent in part from the description, or may be learned by practicing embodiments of the invention. The objects and other advantages of the invention may be realized and obtained by means of the structures particularly pointed out in the written description and the accompanying drawings. Attached Figure Description

[0016] To more clearly illustrate the technical solution of the present invention, the accompanying drawings used in the present invention will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained from these drawings without creative effort.

[0017] Figure 1 A flowchart illustrating a quantum mechanics-based long-term drug kinetics simulation method provided in an embodiment of the present invention; Figure 2 A flowchart illustrating the construction and data preprocessing of the initial complex system for a quantum mechanics-based long-term drug kinetics simulation method provided in this embodiment of the invention. Figure 3A high-level flowchart of a quantum mechanics-based long-term drug kinetics simulation method provided in an embodiment of the present invention. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0019] It should be noted that similar reference numerals and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures. Furthermore, in the description of this invention, terms such as "first," "second," etc., are used only to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0020] In this invention, substances participating in the reaction, such as proteins, drug molecules, and water, are placed in a virtual simulation box to simulate the reaction. The core region where the drug and protein bind is designated as a quantum mechanical (QM) region, and the remaining portion as a molecular mechanical (MM) region. The former requires precise calculations of the movement of atomic nuclei and electrons, while the latter uses a coarse force field to describe the atomic connections. By adding biases, the drug-protein complex is made to attempt different conformations more quickly. Under constant temperature and pressure, the positions and velocities of all atoms are updated step-by-step in extremely short time steps (femtoseconds), and the electronic structure of the QM region is recalculated at each step, thereby simulating atomic trajectory data. Based on this trajectory data, key conformational change events can be extracted. The technical solution of this invention is described in detail below.

[0021] Please refer to Figure 1 , Figure 1 This is a flowchart illustrating a quantum mechanics-based long-term drug kinetics simulation according to an embodiment of the present invention. This quantum mechanics-based long-term drug kinetics simulation method is used in terminal devices, such as mobile phones and computers. The quantum mechanics-based long-term drug kinetics simulation method includes the following steps: S11. Obtain the target protein crystal structure coordinate file data and drug molecule structure data and perform preprocessing to obtain initial complex conformation data and construct a periodic simulation box, thereby obtaining coordinate data of the complete solvation and ion addition system and performing hierarchical energy minimization processing to obtain a stable initial atomic coordinate and topological parameter dataset. S12. Perform atomic-level partitioning based on the stabilized initial atomic coordinates and topology parameter dataset, then obtain the partition index dataset and perform calculation and coupling, then obtain the QM / MM coupling total energy and the full atomic coupling force dataset and encapsulate it to obtain the frame-by-frame force source dataset. S13. Normalize and scale the frame-by-frame force source dataset to obtain the collective variable spatial features and adaptive parameter dataset, apply an adaptive Gaussian bias potential, and then obtain the enhanced sampling dynamics driving dataset. S14. Perform long-term simulation based on the enhanced sampling dynamics driven dataset to obtain a time-series trajectory dataset and perform quality control monitoring, thereby obtaining a quality-controlled qualified long-term atomic trajectory time-series dataset. S15. The quality control qualified long-term atomic trajectory time series dataset is processed in combination with the enhanced sampling dynamics driven dataset to obtain the free energy extreme value and steady-state feature dataset, and principal component analysis and cluster analysis are performed to obtain the transition path and transition state dataset and perform interaction time series analysis, and finally obtain the drug long-term dynamics analysis dataset.

[0022] Based on the above processing flow, through the construction and data preprocessing of the initial complex system, QM / MM partitioning and quantum mechanical force data generation, enhanced sampling and collective variable data-driven processing, parallel acquisition and quality control of long-term trajectory data, and free energy surface reconstruction and conformation data extraction, efficient and high-precision sampling of rare conformational events of protein-drug complexes is achieved, thereby improving the accuracy and completeness of drug action mechanism analysis.

[0023] Please refer to Figure 2 , Figure 2 This is a flowchart illustrating the construction and data preprocessing of the initial complex system in the quantum mechanics-based long-term drug kinetics simulation method according to an embodiment of the present invention. According to this embodiment, the process of acquiring target protein crystal structure coordinate file data and drug molecule structure data and performing preprocessing to obtain initial complex conformation data and constructing a periodic simulation box, thereby obtaining coordinate data of the complete solvation-ionization system and performing hierarchical energy minimization processing to obtain a stabilized initial atomic coordinate and topological parameter dataset, specifically: The target protein crystal structure coordinate file data is obtained and preprocessed to obtain a pure atomic coordinate dataset of the protein. Acquire drug molecular structure data and optimize it to obtain a low-energy three-dimensional conformation; Preprocessing is performed based on the low-energy three-dimensional conformation to obtain a three-dimensional coordinate dataset of drug ligands; Molecular docking was performed based on the protein pure atomic coordinate dataset and the drug ligand three-dimensional coordinate dataset to obtain initial complex conformation data. A periodic simulation box is constructed based on the initial complex conformation data to obtain an all-atom solvation coordinate dataset; Based on the aforementioned all-atom solvation coordinate dataset, the coordinate data of the complete solvation and ion addition system are obtained through calculation. Based on the coordinate data of the complete solvation-ion addition system, a hierarchical energy minimization process is performed to obtain a stable initial atomic coordinate and topological parameter dataset.

[0024] The process involves preprocessing target protein crystal structure coordinate files obtained from PDB databases or homologous sequence construction. This includes removing redundant water of crystallization, non-functional small molecules, redundant chains, and repeating atoms; repairing missing main and side chain residues; filling in missing atoms; correcting unreasonable bond lengths and angles; labeling active sites and binding pocket regions; and obtaining a cleaned and repaired pure atomic coordinate dataset of the protein. SMILES, SDF, or mol2 format structural data of the drug molecule are then acquired and optimized to obtain a low-energy three-dimensional conformation. This low-energy three-dimensional conformation is then preprocessed, and Gasteiger or AM1-BCC charges are assigned. By supplementing polar hydrogen atoms and defining docking characteristic parameters such as rotatable bonds, hydrogen bond donors and acceptors, and hydrophobic centers, a three-dimensional coordinate dataset of drug ligands with charge and flexibility parameters is obtained. Molecular docking is performed based on the protein pure atomic coordinate dataset and the drug ligand three-dimensional coordinate dataset. A docking grid is defined with the protein active site and binding pocket as the center, and the grid size covers the pocket and extends outward by 0.5–1.0 nm. The ligand coordinates are input into the docking program for global search and local refinement. The conformation with the lowest binding energy is retained, and poses with spatial conflicts greater than the threshold are removed. The optimal conformation is retained and merged with the protein coordinates to obtain the initial conformational data of the protein-drug complex.

[0025] Based on the initial conformational data of the protein-drug complex, a cubic or truncated octahedral periodic simulation box was constructed. The initial complex was placed at the center of the simulation system, with the closest distance between the box boundary and the complex surface ≥1.2 nm to ensure no periodic mirror self-interactions. An explicit water molecule model was filled into the box, covering all voids and the protein and drug surfaces, obtaining a full-atom solvation coordinate dataset containing solvent. The total charge of the entire system was calculated based on this dataset. Na⁺ and Cl⁻ ion coordinates were added according to the principle of electroneutrality, further supplementing ions to a physiological concentration of 0.15 mol / L, uniformly distributed in the solvent region and far from the binding pocket. The minimum distance between ions and proteins / drugs was ≥0.5 nm to avoid initial strong electrostatic conflicts, obtaining complete solvation and ion coordinate data. Based on this complete solvation and ion coordinate data, hierarchical energy minimization processing was performed to assign compatible force fields to the protein, drug ligand, water, and ions. CHARMM36 and AMBER were used for the protein. Protein force fields such as ff19SB or CHARMM22 are used. The drug ligands are topologically and in terms of parameters generated using Anttecamber and ParamChem. Water is modeled using TIP3P and SPC / E, and ions are modeled using non-bonded parameters. A dataset of topological parameters, bond parameters, non-bonded parameters, and atom types for the entire system is generated. In the first stage, the protein backbone and heavy atoms of the drug molecule are constrained, and only the solvent water molecules, ions, and amino acid side chains are optimized. The iteration continues until energy convergence is achieved, eliminating local spatial collisions. In the second stage, the energy of the entire system is minimized without constraints. The steepest descent method and the conjugate gradient method are combined until the maximum force is less than 1000 kJ / (mol·nm). All constraints are then released, eliminating all unreasonable atomic contacts and internal tensions, and obtaining a stable initial atomic coordinate and topological parameter dataset that is free of structural conflicts and thermodynamically stable.

[0026] According to an embodiment of the present invention, the step of performing atomic-level partitioning based on the stabilized initial atomic coordinates and topological parameter dataset, then obtaining a partition index dataset and performing calculations and coupling, and then obtaining the QM / MM coupling total energy and the all-atomic coupling force dataset and encapsulating them to obtain a frame-by-frame force source dataset, specifically involves: Atomic partitioning is performed based on the stabilized initial atomic coordinates and topological parameter dataset to obtain a partition index dataset; The partitioned index dataset includes the QM atomic index table, the MM atomic index table, and the QM / MM boundary atomic index table; The boundary correction parameter dataset is obtained by calculating the QM / MM boundary atom index table using a preset boundary processing model. The QM atom index table is calculated using a pre-defined theoretical model to obtain QM coupling data, including QM energy, QM atom forces, and QM region electron density and charge datasets. The MM atom index table is calculated using a preset molecular force field model to obtain MM coupling data, including MM energy, MM atom force and non-bonded interaction parameter datasets. The QM coupling data, MM coupling data and boundary correction parameter dataset are coupled by a preset coupling model to obtain the total QM / MM coupling energy and the all-atom coupling force dataset. The QM / MM coupled total energy and the all-atom coupled force dataset are encapsulated to obtain a frame-by-frame force source dataset, which includes atom number, coordinates, force vector, potential energy components and boundary correction terms.

[0027] Specifically, atomic-level partitioning was performed based on the stabilized initial atomic coordinates and topological parameter dataset. All heavy atoms and polar hydrogen atoms of the drug molecule were included in the quantum mechanical (QM) region, preserving the complete chemical bonds and electronic structure. All atoms of the side chains of key residues at active sites and allosteric sites were included in the QM region, while the main chain atoms were retained in the molecular mechanical (MM) region. Interfacial water molecules within 0.3 nm of the drug's active group that participate in hydrogen bonding networks or proton transfer were included in the QM region. The protein backbone, distal residues, bulk solvent, and background ions were all included in the MM region. This yielded a partitioned index dataset, including a QM atomic index table, a MM atomic index table, and a QM / MM boundary atomic index table.

[0028] The QM / MM boundary atom index table is calculated using a pre-defined boundary treatment model, such as the embedded linked atom method. Virtual hydrogen atoms are inserted at the covalent bond truncation points in the QM and MM regions to replace boundary carbon atoms, saturating the QM region valence bonds. Boundary geometric constraints are applied to maintain the consistency of boundary bond lengths, bond angles, and dihedral angles with the natural configuration, avoiding boundary distortion. Boundary charge correction terms are calculated to eliminate spurious charges and spurious dipoles introduced by truncation, obtaining a boundary correction parameter dataset. Real-time single-point energy and gradient calculations of the QM atom index table are performed using a pre-defined theoretical model, such as density functional theory. The functionals used are B3LYP, M06-2X, or ωB97X-D, and the basis sets are 6-31G(d), 6-31+G(d,p), or def2-SVP. Implicit solvent polarization effects are considered, and the SMD / PCM solvation model is enabled. Self-consistent field iterations are performed in each simulation step, with a convergence threshold ≤10⁻. 6Hartree obtains QM coupling data for each step, including QM energy, QM atomic forces, and QM region electron density and charge datasets. It calculates the MM atom index table using a pre-defined molecular force field model, such as one compatible with classical molecular force fields. For proteins, CHARMM36, AMBERff19SB, or CHARMM22 are used; for solvents and ions, TIP3P, SPC / E, or water model-matched ion parameters are used. Bond stretching, angular bending, dihedral torsion, van der Waals interactions, and electrostatic interactions are calculated. This yields MM coupling data for each step, including MM energy, MM atomic forces, and non-bonded interaction parameter datasets. A pre-defined coupling model, such as an embedded QM / MM Hamiltonian coupling model, couples the QM coupling data, MM coupling data, and boundary correction parameter datasets. The total energy is calculated. The force is the sum of QM energy, MM energy, and QM / MM cross-boundary energy. The M / MM cross-boundary energy includes van der Waals interactions, electrostatic interactions, boundary constraint energies, and charge mapping correction terms between QM and MM. The atomic force is the vector sum of QM force, MM force, and QM / MM cross-boundary force. Each simulation step is strictly synchronously calculated without delay, interpolation, or approximate replacement. This yields the total QM / MM coupling energy and all-atom coupling force dataset for each step. The QM / MM coupling energy and all-atom coupling force for each step are encapsulated according to the dynamics engine format to obtain a frame-by-frame force source dataset that can directly drive molecular dynamics. This dataset includes atom indices, coordinates, force vectors, potential energy components, and boundary correction terms. This dataset is the sole input force source data for subsequent meta-dynamics enhancement sampling and does not introduce other external forces or correction terms.

[0029] According to an embodiment of the present invention, the step of normalizing and scaling the frame-by-frame force source dataset to obtain the spatial features of collective variables and the adaptive parameter dataset, applying an adaptive Gaussian bias potential, and then obtaining the enhanced sampling dynamics-driven dataset specifically involves: Based on the frame-by-frame force source dataset, normalization and scale calibration are performed to obtain collective variable definitions and physical characterization data; Pre-simulation is performed based on the definition of collective variables and physical representation data to obtain the spatial characteristics of collective variables and adaptive parameter dataset; An adaptive Gaussian bias potential is applied based on the spatial characteristics of the collective variables and the adaptive parameter dataset to obtain the bias potential; The bias force and free energy increment are calculated based on the bias potential and collective variables. The enhanced sampling dynamics-driven dataset is obtained by superimposing the bias force and the frame-by-frame force source dataset.

[0030] The process involves normalizing and scaling the frame-by-frame force source dataset, loading system topology, atom indices, QM / MM partitioning information, and boundary constraint parameters, initializing the bias potential history library, free energy surface cache, and collective variable statistical array, establishing a meta-dynamics runtime context, and completing real-time data integration between the enhanced sampling module and the QM / MM calculation module. Based on the conformational transition mechanism of the protein-drug complex, at least one preferred two-dimensional coupled collective variable CV is selected from the frame-by-frame force source dataset. The first collective variable CV1 is the Euclidean distance between the centroid of the drug molecule and the Cα atom of the key protein residue, representing the binding depth and pocket opening degree. The second collective variable CV2 is a weighted combination of the dihedral angles φ / ψ of the main chain of the protein functional loop region, representing conformational inversion and skeletal deformation. The collective variable CV is normalized and scaled to obtain the definition and physical characterization data of the collective variable CV, including the CV definition parameters, atom indices, weight coefficients, and value range dataset.

[0031] Based on the definition of collective variables and physical characterization data, short-duration QM / MM pre-simulations are performed. CV time-series data are collected, and the mean, standard deviation, sampling range, and distribution density of CV are statistically analyzed. The width, height, deposition interval, and effective exploration range of the Gaussian bias potential are determined to avoid conformational distortion due to an excessively strong bias potential or insufficient sampling due to an excessively weak bias potential. The spatial characteristics of the collective variable CV and the adaptive parameter dataset are obtained. Based on the spatial characteristics of the collective variable CV and the adaptive parameter dataset, an adaptive Gaussian bias potential is applied. Within the low-dimensional space spanned by the collective variable CV, a Gaussian bias potential is deposited at fixed time intervals, in the form V(s,t)=∑h·exp[-(ss) k[2σ²], height h is 1.0 kJ / mol, which can adaptively decay with sampling density, width σ is 0.2 times the characteristic length of collective variable CV to ensure smoothness without jagged edges, deposition interval is applied once every 100 simulation steps, strictly synchronized with QM / MM force calculation, real-time recording of the deposited bias potential, including position, intensity and timestamp, updating historical potential energy surface, calculation based on bias potential and collective variable, according to the atomic contribution coefficient and vector direction of collective variable CV, the bias force is distributed to all atoms participating in CV calculation, the bias force acting on the atom is obtained by numerically differentiating the bias potential with respect to the collective variable, synchronously calculate the free energy increment of the current conformation for subsequent free energy surface reconstruction, the bias force only acts on conformation driving, does not change the electronic structure of QM region and QM / MM coupling relationship, according to the bias force and frame-by-frame force source dataset superimposed, the actual The atomic-level bias force calculated in real time is vector-superimposed with the QM / MM coupled all-atom forces in the frame-by-frame force source dataset, F_total(r,t)=F_QM / MM(r,t)+F_bias(CV,t). The superposition process maintains a strict one-to-one correspondence between the direction, magnitude, and atomic index of the force, without introducing additional constraints, modifying the energy term, or interpolating approximations, thus ensuring dynamic consistency. The superimposed total force, total potential energy, real-time value of the collective variable CV, bias potential state, timestamp, and other information are uniformly encapsulated to obtain an enhanced sampling dynamic driving dataset, including the total force vector of all atoms, the total potential energy of the system (QM / MM potential energy + bias potential), real-time CV value and free energy increment, time step, and sampling state identifier. This dataset is the sole driving data for subsequent long-term trajectory parallel acquisition and is used to achieve efficient and accelerated exploration of conformation space under the high-precision QM / MM framework.

[0032] According to an embodiment of the present invention, the step of performing long-term simulation based on the enhanced sampling dynamics driven dataset to obtain a time-series trajectory dataset and performing quality control monitoring, thereby obtaining a quality-controlled qualified long-term atomic trajectory time-series dataset, specifically includes: Long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain dynamic process parameters. Based on the dynamic process parameters, a time-series trajectory dataset is obtained; Quality control monitoring is performed based on the time-series trajectory dataset to obtain a quality-controlled, qualified long-term atomic trajectory time-series dataset.

[0033] The simulation involved long-term time-history simulations based on the enhanced sampling dynamics-driven dataset. This included loading the entire system's topology file, atomic coordinates, QM / MM partition index, boundary constraint parameters, elementary dynamic bias potential history, and collective variable states. Synchronous initialization of the simulation engine, QM computational kernel, and enhanced sampling module was completed, establishing a unified timestamp and iteration counter to ensure no data misalignment. The simulation was conducted using an NPT isothermal and isobaric ensemble at a temperature of 300K, utilizing Velocity. A rescaling or Langevin temperature controller with a coupling time constant τ_T = 0.1–0.5 ps and a pressure of 1 bar is used. A Parrinello–Rahman or Berendsen pressure bath with a coupling time constant τ_P = 1.0–5.0 ps is employed. The integration algorithm uses a leapfrog integrator with an integration step size of 0.5 fs (femtosecond level), adapted to the high-frequency electronic structure calculation and chemical bond vibration stability requirements of the QM region. The LINCS algorithm is used to constrain hydrogen bond lengths, allowing a balance between step size stability and sampling efficiency. The evolution is gradual and synchronous, with QM calculation and MM evolution strictly framed. Within each time step, execution is carried out in a fixed order, without omission, desynchronization, or delay. Based on the current coordinates, the QM kernel is called to complete the real-time calculation of single-point energy and atomic forces in the QM region. Based on the force field and non-bonded interactions, the energy and force calculations in the MM region are completed. The elemental dynamic bias force is superimposed to obtain the total atomic forces. Atomic velocities and positions are updated according to Newton's equations of motion, completing one dynamic iteration to obtain dynamic process parameters. The collective variable CV is updated synchronously. The bias potential library and free energy cache ensure high-precision synchronous calculations at every step, without using interpolation forces, approximation forces, or cache forces as substitutes. Calculations are performed based on dynamic process parameters, employing a CPU / GPU heterogeneous parallel architecture to perform large-scale calculations with a parallel scale of ≥256 CPU cores. It supports hybrid parallelism of MPI multi-process + OpenMP multi-thread, decomposes parallelism by atomic domain, uses spatial decomposition for the MM region, and uses independent computing core binding for the QM region. Asynchronous buffering is used for trajectory writing to avoid disk I / O blocking the simulation process. The checkpoint file is automatically saved and restarted every 10ns, and seamless continuation of calculations is supported after abnormal interruptions, ensuring stable and continuous operation for a long time of 50 microseconds. Trajectory data is collected according to fixed rules to obtain a time-series trajectory dataset. The trajectory is saved every 10ps, and the saved content includes atomic number, three-dimensional coordinates, system potential energy, temperature, pressure, real-time CV value, and timestamp. The format is XTC / TRR format compatible with standard MD analysis tools, balancing accuracy and storage space. The cumulative duration is 50 microseconds of continuous operation, with a total number of trajectory frames ≥5000 frames, forming a complete long-term evolution record.

[0034] Quality control monitoring was performed based on the time-series trajectory dataset. The simulation was conducted online, monitoring three core indicators throughout: energy conservation (total potential and kinetic energy drift <0.5%, exceeding which indicates numerical instability); charge drift (total charge fluctuation in the QM region <±0.01e, preventing abnormal electronic structure calculations); and temperature fluctuation (system temperature fluctuation range <±5K, ensuring normal temperature control; exceeding these limits immediately triggered alarms and automatic corrections; automatically rolling back to the previous checkpoint, reducing the step size, or enhancing SCF convergence before restarting). After simulation completion, global quality control verification was performed on all trajectory frames, removing abnormal frames and supplementing valid segments to obtain a quality-controlled long-term atomic trajectory time-series dataset. This dataset included a complete 50-microsecond time-series coordinate trajectory, energy, temperature, pressure, and CV time-series curves, a quality control report, and simulation logs. This data served as the sole input for subsequent free energy surface reconstruction and key conformational event extraction.

[0035] According to an embodiment of the present invention, the process of combining the quality control qualified long-term atomic trajectory time series dataset with the enhanced sampling kinetics-driven dataset to obtain a free energy extremum and steady-state feature dataset, and performing principal component analysis and cluster analysis, thereby obtaining a transition path and transition state dataset and performing interaction time series analysis, ultimately obtaining a long-term drug kinetics analytical dataset, specifically: Preprocessing is performed on the quality-controlled qualified long-term atomic trajectory time series dataset to obtain standardized trajectory and collective variable time series datasets; The dataset is reweighted based on the enhanced sampling dynamics driving dataset to obtain a two-dimensional free energy surface dataset. Based on the standardized trajectory and collective variable time series dataset combined with the two-dimensional free energy surface dataset, local extremum search and stable conformation state identification are performed to obtain free energy extremum and steady-state feature datasets. Principal component analysis and cluster analysis were performed based on the free energy extreme value and steady-state characteristic dataset to obtain conformational clusters and representative structure datasets; Based on the conformational clusters and representative structure datasets, conformational transition paths and transition states are extracted to obtain a dataset of transition paths and transition states. Based on the interaction time series analysis of the transition path and the transition state dataset, a long-term drug kinetic analysis dataset is obtained.

[0036] The process involves preprocessing a quality-controlled, qualified long-term atomic trajectory time-series dataset. Valid frames retained due to quality control filtering are removed. The time axis and CV sequence are aligned, and NaN values, anomalous jump points, and non-physical mutation frames are removed to obtain standardized trajectory and collective variable (CV) time-series datasets. These datasets are then combined with an augmented sampling dynamics-driven dataset for reweighting. The accumulated Gaussian bias potential history and CV access records from the augmented sampling dynamics-driven dataset are read and subjected to meta-dynamic standard reweighting. The space composed of the two-dimensional collective variables (CV1, CV2) is meshed, with the mesh resolution set according to the characteristic scale of the CV variables. An exponential function is applied based on the accumulated bias potential value of each frame. The bias field introduced by the elementary dynamics is reweighted and corrected. The probability density of the conformational distribution after correction is statistically analyzed. The free energy is calculated according to the formula G(CV)=-kBT·ln[P(CV)], with the unit uniformly set to kJ / mol or kcal / mol. The grid noise is eliminated by smoothing, and a two-dimensional free energy surface dataset is obtained. Local extremum search and stable conformational state identification are performed based on the two-dimensional free energy surface data. The free energy minimum point (energy trap) is located, corresponding to the stable conformational state. The free energy maximum point (energy barrier) is located, corresponding to the conformational transition state. The energy barrier height, relative free energy difference and occupancy probability between each steady state are calculated. The coordinate range and time window of the key state in CV space are marked, and a free energy extremum and steady state feature dataset is obtained.

[0037] Principal component analysis and cluster analysis were performed on the free energy extremum and steady-state characteristic dataset. Principal component analysis was performed on proteins binding pockets, drug molecules, and functional loop atoms. The first two to three principal components were extracted to characterize global conformational changes. K-means clustering or DBSCAN density clustering was used, and the number of clusters was determined based on the number of free energy minima. The central structure of each cluster was taken as the representative conformation. The root mean square deviation (RMSD), occupancy ratio, and residence time within each cluster were calculated to obtain a dataset of conformational clusters and representative structures, including cluster labels, representative structures, cluster center coordinates, and cluster statistical characteristics. Conformational transition paths and transition states were extracted based on the conformational cluster data. Conformational evolution was tracked along the time axis and CV space. Continuous transition segments from one free energy minima to another were selected, and keyframe structures and atomic coordinates of the transition initiation state, transition state, and final state were extracted. The time point, duration, energy barrier crossing rate, and cooperative motion sequence of the transition were determined to obtain a dataset of transition paths and transition states, including conformational transition paths. The transition state structure and transition time series dataset is used for interaction time series analysis based on the transition pathway and transition state dataset. Interactions between representative conformations and transition pathways are analyzed, and the occupancy and bond length / bond angle distribution of hydrogen bonds, salt bridges, hydrophobic interactions, π-stacking, and covalent bonds between the drug and protein are statistically analyzed. Changes in the dihedral angle, pocket volume, and hydrophobic surface area of ​​key residue side chains over time are calculated. Leading events, signal transduction pathways, and allosteric coupling residues triggering conformational transitions are located. The dataset outputs key interaction time series, occupancy, bond parameters, and allosteric transduction pathways, and integrates free energy surface, conformational clusters, transition pathways, and interaction data to complete the analysis of the drug action mechanism. This clarifies conformational transition rules, drug binding modes, stable state distribution, and energy barrier regulation mechanisms, resulting in a long-term drug kinetic analysis dataset, including two-dimensional free energy surface data, free energy extrema, representative conformational structures, transition pathways, key interaction time series, and quality control reports. This dataset is the final output of this simulation and can be directly used for drug molecule optimization and mechanism of action elucidation.

[0038] Please refer to Figure 3 , Figure 3 This is a high-level flowchart of the long-term drug dynamics simulation method based on quantum mechanics in the embodiments of the present invention.

[0039] This invention also discloses a quantum mechanics-based long-term drug kinetics simulation system, including a memory and a processor. The memory includes a quantum mechanics-based long-term drug kinetics simulation method program. When the processor executes the quantum mechanics-based long-term drug kinetics simulation method program, it performs the following steps: The target protein crystal structure coordinate file data and drug molecule structure data are acquired and preprocessed to obtain initial complex conformation data and construct a periodic simulation box. Then, the coordinate data of the complete solvation and ion addition system are obtained and hierarchical energy minimization processing is performed to obtain a stable initial atomic coordinate and topological parameter dataset. Atomic partitioning is performed based on the stabilized initial atomic coordinates and topology parameter dataset. Then, the partition index dataset is obtained and calculated and coupled. Finally, the total energy of QM / MM coupling and the full atomic coupling force dataset are obtained and encapsulated to obtain the frame-by-frame force source dataset. Based on the frame-by-frame force source dataset, normalization and scale calibration are performed to obtain the spatial features of collective variables and the adaptive parameter dataset. An adaptive Gaussian bias potential is then applied, and the enhanced sampling dynamics driving dataset is obtained. Long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain a time-series trajectory dataset and perform quality control monitoring, thereby obtaining a quality-controlled qualified long-term atomic trajectory time-series dataset. The data is processed based on the quality control qualified long-term atomic trajectory time series dataset and the enhanced sampling dynamics driven dataset to obtain the free energy extremum and steady-state feature dataset. Principal component analysis and cluster analysis are then performed to obtain the transition path and transition state dataset and perform interaction time series analysis, ultimately obtaining the drug long-term dynamics analysis dataset.

[0040] Based on the above processing flow, through the construction and data preprocessing of the initial complex system, QM / MM partitioning and quantum mechanical force data generation, enhanced sampling and collective variable data-driven processing, parallel acquisition and quality control of long-term trajectory data, and free energy surface reconstruction and conformation data extraction, efficient and high-precision sampling of rare conformational events of protein-drug complexes is achieved, thereby improving the accuracy and completeness of drug action mechanism analysis.

[0041] According to an embodiment of the present invention, the acquisition of target protein crystal structure coordinate file data and drug molecule structure data, followed by preprocessing to obtain initial complex conformation data and construct a periodic simulation box, and then obtaining coordinate data of the complete solvation-ion addition system and performing hierarchical energy minimization processing to obtain a stable initial atomic coordinate and topological parameter dataset, specifically: The target protein crystal structure coordinate file data is obtained and preprocessed to obtain a pure atomic coordinate dataset of the protein. Acquire drug molecular structure data and optimize it to obtain a low-energy three-dimensional conformation; Preprocessing is performed based on the low-energy three-dimensional conformation to obtain a three-dimensional coordinate dataset of drug ligands; Molecular docking was performed based on the protein pure atomic coordinate dataset and the drug ligand three-dimensional coordinate dataset to obtain initial complex conformation data. A periodic simulation box is constructed based on the initial complex conformation data to obtain an all-atom solvation coordinate dataset; Based on the aforementioned all-atom solvation coordinate dataset, the coordinate data of the complete solvation and ion addition system are obtained through calculation. Based on the coordinate data of the complete solvation-ion addition system, a hierarchical energy minimization process is performed to obtain a stable initial atomic coordinate and topological parameter dataset.

[0042] The process involves preprocessing target protein crystal structure coordinate files obtained from PDB databases or homologous sequence construction. This includes removing redundant water of crystallization, non-functional small molecules, redundant chains, and repeating atoms; repairing missing main and side chain residues; filling in missing atoms; correcting unreasonable bond lengths and angles; labeling active sites and binding pocket regions; and obtaining a cleaned and repaired pure atomic coordinate dataset of the protein. SMILES, SDF, or mol2 format structural data of the drug molecule are then acquired and optimized to obtain a low-energy three-dimensional conformation. This low-energy three-dimensional conformation is then preprocessed, and Gasteiger or AM1-BCC charges are assigned. By supplementing polar hydrogen atoms and defining docking characteristic parameters such as rotatable bonds, hydrogen bond donors and acceptors, and hydrophobic centers, a three-dimensional coordinate dataset of drug ligands with charge and flexibility parameters is obtained. Molecular docking is performed based on the protein pure atomic coordinate dataset and the drug ligand three-dimensional coordinate dataset. A docking grid is defined with the protein active site and binding pocket as the center, and the grid size covers the pocket and extends outward by 0.5–1.0 nm. The ligand coordinates are input into the docking program for global search and local refinement. The conformation with the lowest binding energy is retained, and poses with spatial conflicts greater than the threshold are removed. The optimal conformation is retained and merged with the protein coordinates to obtain the initial conformational data of the protein-drug complex.

[0043] Based on the initial conformational data of the protein-drug complex, a cubic or truncated octahedral periodic simulation box was constructed. The initial complex was placed at the center of the simulation system, with the closest distance between the box boundary and the complex surface ≥1.2 nm to ensure no periodic mirror self-interactions. An explicit water molecule model was filled into the box, covering all voids and the protein and drug surfaces, obtaining a full-atom solvation coordinate dataset containing solvent. The total charge of the entire system was calculated based on this dataset. Na⁺ and Cl⁻ ion coordinates were added according to the principle of electroneutrality, further supplementing ions to a physiological concentration of 0.15 mol / L, uniformly distributed in the solvent region and far from the binding pocket. The minimum distance between ions and proteins / drugs was ≥0.5 nm to avoid initial strong electrostatic conflicts, obtaining complete solvation and ion coordinate data. Based on this complete solvation and ion coordinate data, hierarchical energy minimization processing was performed to assign compatible force fields to the protein, drug ligand, water, and ions. CHARMM36 and AMBER were used for the protein. Protein force fields such as ff19SB or CHARMM22 are used. The drug ligands are topologically and in terms of parameters generated using Anttecamber and ParamChem. Water is modeled using TIP3P and SPC / E, and ions are modeled using non-bonded parameters. A dataset of topological parameters, bond parameters, non-bonded parameters, and atom types for the entire system is generated. In the first stage, the protein backbone and heavy atoms of the drug molecule are constrained, and only the solvent water molecules, ions, and amino acid side chains are optimized. The iteration continues until energy convergence is achieved, eliminating local spatial collisions. In the second stage, the energy of the entire system is minimized without constraints. The steepest descent method and the conjugate gradient method are combined until the maximum force is less than 1000 kJ / (mol·nm). All constraints are then released, eliminating all unreasonable atomic contacts and internal tensions, and obtaining a stable initial atomic coordinate and topological parameter dataset that is free of structural conflicts and thermodynamically stable.

[0044] According to an embodiment of the present invention, the step of performing atomic-level partitioning based on the stabilized initial atomic coordinates and topological parameter dataset, then obtaining a partition index dataset and performing calculations and coupling, and then obtaining the QM / MM coupling total energy and the all-atomic coupling force dataset and encapsulating them to obtain a frame-by-frame force source dataset, specifically involves: Atomic partitioning is performed based on the stabilized initial atomic coordinates and topological parameter dataset to obtain a partition index dataset; The partitioned index dataset includes the QM atomic index table, the MM atomic index table, and the QM / MM boundary atomic index table; The boundary correction parameter dataset is obtained by calculating the QM / MM boundary atom index table using a preset boundary processing model. The QM atom index table is calculated using a pre-defined theoretical model to obtain QM coupling data, including QM energy, QM atom forces, and QM region electron density and charge datasets. The MM atom index table is calculated using a preset molecular force field model to obtain MM coupling data, including MM energy, MM atom force and non-bonded interaction parameter datasets. The QM coupling data, MM coupling data and boundary correction parameter dataset are coupled by a preset coupling model to obtain the total QM / MM coupling energy and the all-atom coupling force dataset. The QM / MM coupled total energy and the all-atom coupled force dataset are encapsulated to obtain a frame-by-frame force source dataset, which includes atom number, coordinates, force vector, potential energy components and boundary correction terms.

[0045] Specifically, atomic-level partitioning was performed based on the stabilized initial atomic coordinates and topological parameter dataset. All heavy atoms and polar hydrogen atoms of the drug molecule were included in the quantum mechanical (QM) region, preserving the complete chemical bonds and electronic structure. All atoms of the side chains of key residues at active sites and allosteric sites were included in the QM region, while the main chain atoms were retained in the molecular mechanical (MM) region. Interfacial water molecules within 0.3 nm of the drug's active group that participate in hydrogen bonding networks or proton transfer were included in the QM region. The protein backbone, distal residues, bulk solvent, and background ions were all included in the MM region. This yielded a partitioned index dataset, including a QM atomic index table, a MM atomic index table, and a QM / MM boundary atomic index table.

[0046] The QM / MM boundary atom index table is calculated using a pre-defined boundary treatment model, such as the embedded linked atom method. Virtual hydrogen atoms are inserted at the covalent bond truncation points in the QM and MM regions to replace boundary carbon atoms, saturating the QM region valence bonds. Boundary geometric constraints are applied to maintain the consistency of boundary bond lengths, bond angles, and dihedral angles with the natural configuration, avoiding boundary distortion. Boundary charge correction terms are calculated to eliminate spurious charges and spurious dipoles introduced by truncation, obtaining a boundary correction parameter dataset. Real-time single-point energy and gradient calculations of the QM atom index table are performed using a pre-defined theoretical model, such as density functional theory. The functionals used are B3LYP, M06-2X, or ωB97X-D, and the basis sets are 6-31G(d), 6-31+G(d,p), or def2-SVP. Implicit solvent polarization effects are considered, and the SMD / PCM solvation model is enabled. Self-consistent field iterations are performed in each simulation step, with a convergence threshold ≤10⁻. 6Hartree obtains QM coupling data for each step, including QM energy, QM atomic forces, and QM region electron density and charge datasets. It calculates the MM atom index table using a pre-defined molecular force field model, such as one compatible with classical molecular force fields. For proteins, CHARMM36, AMBERff19SB, or CHARMM22 are used; for solvents and ions, TIP3P, SPC / E, or water model-matched ion parameters are used. Bond stretching, angular bending, dihedral torsion, van der Waals interactions, and electrostatic interactions are calculated. This yields MM coupling data for each step, including MM energy, MM atomic forces, and non-bonded interaction parameter datasets. A pre-defined coupling model, such as an embedded QM / MM Hamiltonian coupling model, couples the QM coupling data, MM coupling data, and boundary correction parameter datasets. The total energy is calculated. The force is the sum of QM energy, MM energy, and QM / MM cross-boundary energy. The M / MM cross-boundary energy includes van der Waals interactions, electrostatic interactions, boundary constraint energies, and charge mapping correction terms between QM and MM. The atomic force is the vector sum of QM force, MM force, and QM / MM cross-boundary force. Each simulation step is strictly synchronously calculated without delay, interpolation, or approximate replacement. This yields the total QM / MM coupling energy and all-atom coupling force dataset for each step. The QM / MM coupling energy and all-atom coupling force for each step are encapsulated according to the dynamics engine format to obtain a frame-by-frame force source dataset that can directly drive molecular dynamics. This dataset includes atom indices, coordinates, force vectors, potential energy components, and boundary correction terms. This dataset is the sole input force source data for subsequent meta-dynamics enhancement sampling and does not introduce other external forces or correction terms.

[0047] According to an embodiment of the present invention, the step of normalizing and scaling the frame-by-frame force source dataset to obtain the spatial features of collective variables and the adaptive parameter dataset, applying an adaptive Gaussian bias potential, and then obtaining the enhanced sampling dynamics-driven dataset specifically involves: Based on the frame-by-frame force source dataset, normalization and scale calibration are performed to obtain collective variable definitions and physical characterization data; Pre-simulation is performed based on the definition of collective variables and physical representation data to obtain the spatial characteristics of collective variables and adaptive parameter dataset; An adaptive Gaussian bias potential is applied based on the spatial characteristics of the collective variables and the adaptive parameter dataset to obtain the bias potential; The bias force and free energy increment are calculated based on the bias potential and collective variables. The enhanced sampling dynamics-driven dataset is obtained by superimposing the bias force and the frame-by-frame force source dataset.

[0048] The process involves normalizing and scaling the frame-by-frame force source dataset, loading system topology, atom indices, QM / MM partitioning information, and boundary constraint parameters, initializing the bias potential history library, free energy surface cache, and collective variable statistical array, establishing a meta-dynamics runtime context, and completing real-time data integration between the enhanced sampling module and the QM / MM calculation module. Based on the conformational transition mechanism of the protein-drug complex, at least one preferred two-dimensional coupled collective variable CV is selected from the frame-by-frame force source dataset. The first collective variable CV1 is the Euclidean distance between the centroid of the drug molecule and the Cα atom of the key protein residue, representing the binding depth and pocket opening degree. The second collective variable CV2 is a weighted combination of the dihedral angles φ / ψ of the main chain of the protein functional loop region, representing conformational inversion and skeletal deformation. The collective variable CV is normalized and scaled to obtain the definition and physical characterization data of the collective variable CV, including the CV definition parameters, atom indices, weight coefficients, and value range dataset.

[0049] Based on the definition of collective variables and physical characterization data, short-duration QM / MM pre-simulations are performed. CV time-series data are collected, and the mean, standard deviation, sampling range, and distribution density of CV are statistically analyzed. The width, height, deposition interval, and effective exploration range of the Gaussian bias potential are determined to avoid conformational distortion due to an excessively strong bias potential or insufficient sampling due to an excessively weak bias potential. The spatial characteristics of the collective variable CV and the adaptive parameter dataset are obtained. Based on the spatial characteristics of the collective variable CV and the adaptive parameter dataset, an adaptive Gaussian bias potential is applied. Within the low-dimensional space spanned by the collective variable CV, a Gaussian bias potential is deposited at fixed time intervals, in the form V(s,t)=∑h·exp[-(ss) k[2σ²], height h is 1.0 kJ / mol, which can adaptively decay with sampling density, width σ is 0.2 times the characteristic length of collective variable CV to ensure smoothness without jagged edges, deposition interval is applied once every 100 simulation steps, strictly synchronized with QM / MM force calculation, real-time recording of the deposited bias potential, including position, intensity and timestamp, updating historical potential energy surface, calculation based on bias potential and collective variable, according to the atomic contribution coefficient and vector direction of collective variable CV, the bias force is distributed to all atoms participating in CV calculation, the bias force acting on the atom is obtained by numerically differentiating the bias potential with respect to the collective variable, synchronously calculate the free energy increment of the current conformation for subsequent free energy surface reconstruction, the bias force only acts on conformation driving, does not change the electronic structure of QM region and QM / MM coupling relationship, according to the bias force and frame-by-frame force source dataset superimposed, the actual The atomic-level bias force calculated in real time is vector-superimposed with the QM / MM coupled all-atom forces in the frame-by-frame force source dataset, F_total(r,t)=F_QM / MM(r,t)+F_bias(CV,t). The superposition process maintains a strict one-to-one correspondence between the direction, magnitude, and atomic index of the force, without introducing additional constraints, modifying the energy term, or interpolating approximations, thus ensuring dynamic consistency. The superimposed total force, total potential energy, real-time value of the collective variable CV, bias potential state, timestamp, and other information are uniformly encapsulated to obtain an enhanced sampling dynamic driving dataset, including the total force vector of all atoms, the total potential energy of the system (QM / MM potential energy + bias potential), real-time CV value and free energy increment, time step, and sampling state identifier. This dataset is the sole driving data for subsequent long-term trajectory parallel acquisition and is used to achieve efficient and accelerated exploration of conformation space under the high-precision QM / MM framework.

[0050] According to an embodiment of the present invention, the step of performing long-term simulation based on the enhanced sampling dynamics driven dataset to obtain a time-series trajectory dataset and performing quality control monitoring, thereby obtaining a quality-controlled qualified long-term atomic trajectory time-series dataset, specifically includes: Long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain dynamic process parameters. Based on the dynamic process parameters, a time-series trajectory dataset is obtained; Quality control monitoring is performed based on the time-series trajectory dataset to obtain a quality-controlled, qualified long-term atomic trajectory time-series dataset.

[0051] The simulation involved long-term time-history simulations based on the enhanced sampling dynamics-driven dataset. This included loading the entire system's topology file, atomic coordinates, QM / MM partition index, boundary constraint parameters, elementary dynamic bias potential history, and collective variable states. Synchronous initialization of the simulation engine, QM computational kernel, and enhanced sampling module was completed, establishing a unified timestamp and iteration counter to ensure no data misalignment. The simulation was conducted using an NPT isothermal and isobaric ensemble at a temperature of 300K, utilizing Velocity. A rescaling or Langevin temperature controller with a coupling time constant τ_T = 0.1–0.5 ps and a pressure of 1 bar is used. A Parrinello–Rahman or Berendsen pressure bath with a coupling time constant τ_P = 1.0–5.0 ps is employed. The integration algorithm uses a leapfrog integrator with an integration step size of 0.5 fs (femtosecond level), adapted to the high-frequency electronic structure calculation and chemical bond vibration stability requirements of the QM region. The LINCS algorithm is used to constrain hydrogen bond lengths, allowing a balance between step size stability and sampling efficiency. The evolution is gradual and synchronous, with QM calculation and MM evolution strictly framed. Within each time step, execution is carried out in a fixed order, without omission, desynchronization, or delay. Based on the current coordinates, the QM kernel is called to complete the real-time calculation of single-point energy and atomic forces in the QM region. Based on the force field and non-bonded interactions, the energy and force calculations in the MM region are completed. The elemental dynamic bias force is superimposed to obtain the total atomic forces. Atomic velocities and positions are updated according to Newton's equations of motion, completing one dynamic iteration to obtain dynamic process parameters. The collective variable CV is updated synchronously. The bias potential library and free energy cache ensure high-precision synchronous calculations at every step, without using interpolation forces, approximation forces, or cache forces as substitutes. Calculations are performed based on dynamic process parameters, employing a CPU / GPU heterogeneous parallel architecture to perform large-scale calculations with a parallel scale of ≥256 CPU cores. It supports hybrid parallelism of MPI multi-process + OpenMP multi-thread, decomposes parallelism by atomic domain, uses spatial decomposition for the MM region, and uses independent computing core binding for the QM region. Asynchronous buffering is used for trajectory writing to avoid disk I / O blocking the simulation process. The checkpoint file is automatically saved and restarted every 10ns, and seamless continuation of calculations is supported after abnormal interruptions, ensuring stable and continuous operation for a long time of 50 microseconds. Trajectory data is collected according to fixed rules to obtain a time-series trajectory dataset. The trajectory is saved every 10ps, and the saved content includes atomic number, three-dimensional coordinates, system potential energy, temperature, pressure, real-time CV value, and timestamp. The format is XTC / TRR format compatible with standard MD analysis tools, balancing accuracy and storage space. The cumulative duration is 50 microseconds of continuous operation, with a total number of trajectory frames ≥5000 frames, forming a complete long-term evolution record.

[0052] Quality control monitoring was performed based on the time-series trajectory dataset. The simulation was conducted online, monitoring three core indicators throughout: energy conservation (total potential and kinetic energy drift <0.5%, exceeding which indicates numerical instability); charge drift (total charge fluctuation in the QM region <±0.01e, preventing abnormal electronic structure calculations); and temperature fluctuation (system temperature fluctuation range <±5K, ensuring normal temperature control; exceeding these limits immediately triggered alarms and automatic corrections; automatically rolling back to the previous checkpoint, reducing the step size, or enhancing SCF convergence before restarting). After simulation completion, global quality control verification was performed on all trajectory frames, removing abnormal frames and supplementing valid segments to obtain a quality-controlled long-term atomic trajectory time-series dataset. This dataset included a complete 50-microsecond time-series coordinate trajectory, energy, temperature, pressure, and CV time-series curves, a quality control report, and simulation logs. This data served as the sole input for subsequent free energy surface reconstruction and key conformational event extraction.

[0053] According to an embodiment of the present invention, the process of combining the quality control qualified long-term atomic trajectory time series dataset with the enhanced sampling kinetics-driven dataset to obtain a free energy extremum and steady-state feature dataset, and performing principal component analysis and cluster analysis, thereby obtaining a transition path and transition state dataset and performing interaction time series analysis, ultimately obtaining a long-term drug kinetics analytical dataset, specifically: Preprocessing is performed on the quality-controlled qualified long-term atomic trajectory time series dataset to obtain standardized trajectory and collective variable time series datasets; The dataset is reweighted based on the enhanced sampling dynamics driving dataset to obtain a two-dimensional free energy surface dataset. Based on the standardized trajectory and collective variable time series dataset combined with the two-dimensional free energy surface dataset, local extremum search and stable conformation state identification are performed to obtain free energy extremum and steady-state feature datasets. Principal component analysis and cluster analysis were performed based on the free energy extreme value and steady-state characteristic dataset to obtain conformational clusters and representative structure datasets; Based on the conformational clusters and representative structure datasets, conformational transition paths and transition states are extracted to obtain a dataset of transition paths and transition states. Based on the interaction time series analysis of the transition path and the transition state dataset, a long-term drug kinetic analysis dataset is obtained.

[0054] The process involves preprocessing a quality-controlled, qualified long-term atomic trajectory time-series dataset. Valid frames retained due to quality control filtering are removed. The time axis and CV sequence are aligned, and NaN values, anomalous jump points, and non-physical mutation frames are removed to obtain standardized trajectory and collective variable (CV) time-series datasets. These datasets are then combined with an augmented sampling dynamics-driven dataset for reweighting. The accumulated Gaussian bias potential history and CV access records from the augmented sampling dynamics-driven dataset are read and subjected to meta-dynamic standard reweighting. The space composed of the two-dimensional collective variables (CV1, CV2) is meshed, with the mesh resolution set according to the characteristic scale of the CV variables. An exponential function is applied based on the accumulated bias potential value of each frame. The bias field introduced by the elementary dynamics is reweighted and corrected. The probability density of the conformational distribution after correction is statistically analyzed. The free energy is calculated according to the formula G(CV)=-kBT·ln[P(CV)], with the unit uniformly set to kJ / mol or kcal / mol. The grid noise is eliminated by smoothing, and a two-dimensional free energy surface dataset is obtained. Local extremum search and stable conformational state identification are performed based on the two-dimensional free energy surface data. The free energy minimum point (energy trap) is located, corresponding to the stable conformational state. The free energy maximum point (energy barrier) is located, corresponding to the conformational transition state. The energy barrier height, relative free energy difference and occupancy probability between each steady state are calculated. The coordinate range and time window of the key state in CV space are marked, and a free energy extremum and steady state feature dataset is obtained.

[0055] Principal component analysis and cluster analysis were performed on the free energy extremum and steady-state characteristic dataset. Principal component analysis was performed on proteins binding pockets, drug molecules, and functional loop atoms. The first two to three principal components were extracted to characterize global conformational changes. K-means clustering or DBSCAN density clustering was used, and the number of clusters was determined based on the number of free energy minima. The central structure of each cluster was taken as the representative conformation. The root mean square deviation (RMSD), occupancy ratio, and residence time within each cluster were calculated to obtain a dataset of conformational clusters and representative structures, including cluster labels, representative structures, cluster center coordinates, and cluster statistical characteristics. Conformational transition paths and transition states were extracted based on the conformational cluster data. Conformational evolution was tracked along the time axis and CV space. Continuous transition segments from one free energy minima to another were selected, and keyframe structures and atomic coordinates of the transition initiation state, transition state, and final state were extracted. The time point, duration, energy barrier crossing rate, and cooperative motion sequence of the transition were determined to obtain a dataset of transition paths and transition states, including conformational transition paths. The transition state structure and transition time series dataset is used for interaction time series analysis based on the transition pathway and transition state dataset. Interactions between representative conformations and transition pathways are analyzed, and the occupancy and bond length / bond angle distribution of hydrogen bonds, salt bridges, hydrophobic interactions, π-stacking, and covalent bonds between the drug and protein are statistically analyzed. Changes in the dihedral angle, pocket volume, and hydrophobic surface area of ​​key residue side chains over time are calculated. Leading events, signal transduction pathways, and allosteric coupling residues triggering conformational transitions are located. The dataset outputs key interaction time series, occupancy, bond parameters, and allosteric transduction pathways, and integrates free energy surface, conformational clusters, transition pathways, and interaction data to complete the analysis of the drug action mechanism. This clarifies conformational transition rules, drug binding modes, stable state distribution, and energy barrier regulation mechanisms, resulting in a long-term drug kinetic analysis dataset, including two-dimensional free energy surface data, free energy extrema, representative conformational structures, transition pathways, key interaction time series, and quality control reports. This dataset is the final output of this simulation and can be directly used for drug molecule optimization and mechanism of action elucidation.

[0056] A third aspect of the present invention provides a readable storage medium comprising a quantum mechanics-based long-term drug kinetic simulation method program, wherein when the quantum mechanics-based long-term drug kinetic simulation method program is executed by a processor, it implements the steps of the quantum mechanics-based long-term drug kinetic simulation method as described in any of the preceding claims.

[0057] According to the technical solution of the present invention, through the construction and data preprocessing of the initial complex system, QM / MM partitioning and quantum mechanical force data generation, enhanced sampling and collective variable data driving, parallel acquisition and quality control of long-term trajectory data, and free energy surface reconstruction and conformation data extraction, efficient and high-precision sampling of rare conformational events of protein and drug complexes is achieved, thereby improving the accuracy and completeness of drug action mechanism analysis.

Claims

1. A long-term drug kinetic simulation method based on quantum mechanics, characterized in that, Includes the following steps: The target protein crystal structure coordinate file data and drug molecule structure data are acquired and preprocessed to obtain initial complex conformation data and construct a periodic simulation box. Then, the coordinate data of the complete solvation and ion addition system are obtained and hierarchical energy minimization processing is performed to obtain a stable initial atomic coordinate and topological parameter dataset. Atomic partitioning is performed based on the stabilized initial atomic coordinates and topology parameter dataset. Then, the partition index dataset is obtained and calculated and coupled. Finally, the total energy of QM / MM coupling and the full atomic coupling force dataset are obtained and encapsulated to obtain the frame-by-frame force source dataset. Based on the frame-by-frame force source dataset, normalization and scale calibration are performed to obtain the spatial features of collective variables and the adaptive parameter dataset. An adaptive Gaussian bias potential is then applied, and the enhanced sampling dynamics driving dataset is obtained. Long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain a time-series trajectory dataset and perform quality control monitoring, thereby obtaining a quality-controlled qualified long-term atomic trajectory time-series dataset. The data is processed based on the quality control qualified long-term atomic trajectory time series dataset and the enhanced sampling dynamics driven dataset to obtain the free energy extremum and steady-state feature dataset. Principal component analysis and cluster analysis are then performed to obtain the transition path and transition state dataset and perform interaction time series analysis, ultimately obtaining the drug long-term dynamics analysis dataset.

2. The long-term drug dynamics simulation method based on quantum mechanics according to claim 1, characterized in that, The process involves acquiring and preprocessing the target protein crystal structure coordinates and drug molecule structure data to obtain initial complex conformation data and constructing a periodic simulation box. This leads to the acquisition of coordinate data for the complete solvation and ionization system, followed by hierarchical energy minimization to obtain a stable initial atomic coordinates and topological parameter dataset. The target protein crystal structure coordinate file data is obtained and preprocessed to obtain a pure atomic coordinate dataset of the protein. Acquire drug molecular structure data and optimize it to obtain a low-energy three-dimensional conformation; Preprocessing is performed based on the low-energy three-dimensional conformation to obtain a three-dimensional coordinate dataset of drug ligands; Molecular docking was performed based on the protein pure atomic coordinate dataset and the drug ligand three-dimensional coordinate dataset to obtain initial complex conformation data. A periodic simulation box is constructed based on the initial complex conformation data to obtain an all-atom solvation coordinate dataset; Based on the aforementioned all-atom solvation coordinate dataset, the coordinate data of the complete solvation and ion addition system are obtained through calculation. Based on the coordinate data of the complete solvation-ion addition system, a hierarchical energy minimization process is performed to obtain a stable initial atomic coordinate and topological parameter dataset.

3. The long-term drug dynamics simulation method based on quantum mechanics according to claim 2, characterized in that, The process involves atomic-level partitioning based on the stabilized initial atomic coordinates and topology parameter dataset, obtaining a partition index dataset, performing calculations and coupling, and then obtaining the QM / MM coupled total energy and the all-atomic coupled force dataset, which are then encapsulated to obtain a frame-by-frame force source dataset, including: Atomic partitioning is performed based on the stabilized initial atomic coordinates and topological parameter dataset to obtain a partition index dataset; The partitioned index dataset includes the QM atomic index table, the MM atomic index table, and the QM / MM boundary atomic index table; The boundary correction parameter dataset is obtained by calculating the QM / MM boundary atom index table using a preset boundary processing model. The QM atom index table is calculated using a pre-defined theoretical model to obtain QM coupling data, including QM energy, QM atom forces, and QM region electron density and charge datasets. The MM atom index table is calculated using a preset molecular force field model to obtain MM coupling data, including MM energy, MM atom force and non-bonded interaction parameter datasets. The QM coupling data, MM coupling data and boundary correction parameter dataset are coupled by a preset coupling model to obtain the total QM / MM coupling energy and the all-atom coupling force dataset. The QM / MM coupled total energy and the all-atom coupled force dataset are encapsulated to obtain a frame-by-frame force source dataset, which includes atom number, coordinates, force vector, potential energy components and boundary correction terms.

4. The long-term drug dynamics simulation method based on quantum mechanics according to claim 1, characterized in that, The process involves normalizing and scaling the frame-by-frame force source dataset to obtain the spatial features of the collective variables and the adaptive parameter dataset, applying an adaptive Gaussian bias potential, and then obtaining the enhanced sampling dynamics-driven dataset, including: Based on the frame-by-frame force source dataset, normalization and scale calibration are performed to obtain collective variable definitions and physical characterization data; Pre-simulation is performed based on the definition of collective variables and physical representation data to obtain the spatial characteristics of collective variables and adaptive parameter dataset; An adaptive Gaussian bias potential is applied based on the spatial characteristics of the collective variables and the adaptive parameter dataset to obtain the bias potential; The bias force and free energy increment are calculated based on the bias potential and collective variables. The enhanced sampling dynamics-driven dataset is obtained by superimposing the bias force and the frame-by-frame force source dataset.

5. The long-term drug dynamics simulation method based on quantum mechanics according to claim 1, characterized in that, The step of performing long-term simulation based on the enhanced sampling dynamics-driven dataset to obtain a time-series trajectory dataset and performing quality control monitoring, thereby obtaining a quality-controlled qualified long-term atomic trajectory time-series dataset, includes: Long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain dynamic process parameters. Based on the dynamic process parameters, a time-series trajectory dataset is obtained; Quality control monitoring is performed based on the time-series trajectory dataset to obtain a quality-controlled, qualified long-term atomic trajectory time-series dataset.

6. The long-term drug dynamics simulation method based on quantum mechanics according to claim 1, characterized in that, The process involves combining the quality-controlled qualified long-term atomic trajectory time-series dataset with the enhanced sampling kinetics-driven dataset to obtain a free energy extremum and steady-state characteristic dataset. Principal component analysis and cluster analysis are then performed to obtain a transition path and transition state dataset, followed by interaction time-series analysis. Finally, a long-term drug kinetics analytical dataset is obtained, including: Preprocessing is performed on the quality-controlled qualified long-term atomic trajectory time series dataset to obtain standardized trajectory and collective variable time series datasets; The dataset is reweighted based on the enhanced sampling dynamics driving dataset to obtain a two-dimensional free energy surface dataset. Based on the standardized trajectory and collective variable time series dataset combined with the two-dimensional free energy surface dataset, local extremum search and stable conformation state identification are performed to obtain free energy extremum and steady-state feature datasets. Principal component analysis and cluster analysis were performed based on the free energy extreme value and steady-state characteristic dataset to obtain conformational clusters and representative structure datasets; Based on the conformational clusters and representative structure datasets, conformational transition paths and transition states are extracted to obtain a dataset of transition paths and transition states. Based on the interaction time series analysis of the transition path and the transition state dataset, a long-term drug kinetic analysis dataset is obtained.

7. A long-term drug dynamics simulation system based on quantum mechanics, characterized in that, The system includes a memory and a processor. The memory contains a program for a quantum mechanics-based long-term drug kinetic simulation method. When the processor executes the program for the quantum mechanics-based long-term drug kinetic simulation method, it performs the following steps: The target protein crystal structure coordinate file data and drug molecule structure data are acquired and preprocessed to obtain initial complex conformation data and construct a periodic simulation box. Then, the coordinate data of the complete solvation and ion addition system are obtained and hierarchical energy minimization processing is performed to obtain a stable initial atomic coordinate and topological parameter dataset. Atomic partitioning is performed based on the stabilized initial atomic coordinates and topology parameter dataset. Then, the partition index dataset is obtained and calculated and coupled. Finally, the total energy of QM / MM coupling and the full atomic coupling force dataset are obtained and encapsulated to obtain the frame-by-frame force source dataset. Based on the frame-by-frame force source dataset, normalization and scale calibration are performed to obtain the spatial features of collective variables and the adaptive parameter dataset. An adaptive Gaussian bias potential is then applied, and the enhanced sampling dynamics driving dataset is obtained. Long-term simulations are performed on the enhanced sampling dynamics-driven dataset to obtain a time-series trajectory dataset and perform quality control monitoring, thereby obtaining a quality-controlled qualified long-term atomic trajectory time-series dataset. The data is processed based on the quality control qualified long-term atomic trajectory time series dataset and the enhanced sampling dynamics driven dataset to obtain the free energy extremum and steady-state feature dataset. Principal component analysis and cluster analysis are then performed to obtain the transition path and transition state dataset and perform interaction time series analysis, ultimately obtaining the drug long-term dynamics analysis dataset.

8. The long-term drug dynamics simulation system based on quantum mechanics according to claim 7, characterized in that, The process involves acquiring and preprocessing the target protein crystal structure coordinates and drug molecule structure data to obtain initial complex conformation data and constructing a periodic simulation box. This leads to the acquisition of coordinate data for the complete solvation and ionization system, followed by hierarchical energy minimization to obtain a stable initial atomic coordinates and topological parameter dataset. The target protein crystal structure coordinate file data is obtained and preprocessed to obtain a pure atomic coordinate dataset of the protein. Acquire drug molecular structure data and optimize it to obtain a low-energy three-dimensional conformation; Preprocessing is performed based on the low-energy three-dimensional conformation to obtain a three-dimensional coordinate dataset of drug ligands; Molecular docking was performed based on the protein pure atomic coordinate dataset and the drug ligand three-dimensional coordinate dataset to obtain initial complex conformation data. A periodic simulation box is constructed based on the initial complex conformation data to obtain an all-atom solvation coordinate dataset; Based on the aforementioned all-atom solvation coordinate dataset, the coordinate data of the complete solvation and ion addition system are obtained through calculation. Based on the coordinate data of the complete solvation-ion addition system, a hierarchical energy minimization process is performed to obtain a stable initial atomic coordinate and topological parameter dataset.

9. The long-term drug dynamics simulation system based on quantum mechanics according to claim 8, characterized in that, The process involves atomic-level partitioning based on the stabilized initial atomic coordinates and topology parameter dataset, obtaining a partition index dataset, performing calculations and coupling, and then obtaining the QM / MM coupled total energy and the all-atomic coupled force dataset, which are then encapsulated to obtain a frame-by-frame force source dataset, including: Atomic partitioning is performed based on the stabilized initial atomic coordinates and topological parameter dataset to obtain a partition index dataset; The partitioned index dataset includes the QM atomic index table, the MM atomic index table, and the QM / MM boundary atomic index table; The boundary correction parameter dataset is obtained by calculating the QM / MM boundary atom index table using a preset boundary processing model. The QM atom index table is calculated using a pre-defined theoretical model to obtain QM coupling data, including QM energy, QM atom forces, and QM region electron density and charge datasets. The MM atom index table is calculated using a preset molecular force field model to obtain MM coupling data, including MM energy, MM atom force and non-bonded interaction parameter datasets. The QM coupling data, MM coupling data and boundary correction parameter dataset are coupled by a preset coupling model to obtain the total QM / MM coupling energy and the all-atom coupling force dataset. The QM / MM coupled total energy and the all-atom coupled force dataset are encapsulated to obtain a frame-by-frame force source dataset, which includes atom number, coordinates, force vector, potential energy components and boundary correction terms.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a quantum mechanics-based long-term drug kinetic simulation method program, which, when executed by a processor, implements the steps of the quantum mechanics-based long-term drug kinetic simulation method as described in any one of claims 1 to 6.

Images