A drug virtual screening strategy and system based on molecular docking and dynamics simulation

Abstract

This invention discloses a drug virtual screening strategy and system based on molecular docking and kinetic simulation, specifically relating to the field of computer-aided design technology. It involves acquiring the molecular dynamics simulation trajectory of target proteins and constructing receptor conformational contact maps to identify receptor conformational pocket opening and closing transition nodes. The system then screens and evaluates the distinguishing ability of each candidate receptor conformation in molecular docking, thereby obtaining complementary target receptor conformation combinations and determining the screening order accordingly. Ultimately, this achieves efficient virtual screening of a candidate compound library. This invention improves the ability of different receptor conformations to distinguish candidate compounds, effectively enhancing the accuracy and reliability of drug virtual screening.

Description

Technical Field

[0001] This invention relates to the field of computer-aided design technology, and more specifically, to a drug virtual screening strategy and system based on molecular docking and kinetic simulation. Background Technology

[0002] Molecular dynamics simulations are widely used to obtain different conformational states of protein receptors and to perform virtual drug screening in order to improve the success rate of drug discovery.

[0003] Existing methods typically use multiple receptor conformations obtained directly as input for virtual screening, lacking means to distinguish and evaluate different receptor conformations. This makes it impossible to reasonably assess the contribution and distinguishing ability of each receptor conformation in the drug molecule screening process, resulting in low overall efficiency of virtual screening and difficulty in fully utilizing the advantages brought by the diversity of protein receptor conformations. Summary of the Invention

[0004] To overcome the aforementioned deficiencies of the prior art, embodiments of the present invention provide a drug virtual screening strategy and system based on molecular docking and kinetic simulation to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] A virtual drug screening strategy based on molecular docking and kinetic simulation includes the following steps:

[0007] S1: Obtain the molecular dynamics simulation trajectory of the target protein and extract the receptor conformation set, and establish a receptor conformation contact map based on the contact relationship of binding pocket residues;

[0008] S2: Identify pocket opening-closing transition nodes according to the receptor conformational contact map, and select candidate receptor conformations from both sides of the pocket opening-closing transition nodes;

[0009] S3: Dock the baseline active molecule set and the decoy molecule set to the candidate receptor conformation to generate a set of binding postures and a set of scores and rankings for the candidate receptor conformation;

[0010] S4: For each candidate receptor conformation, calculate the differences in burial depth distribution and key contact rearrangement times in the binding posture set, and combine them with the scoring and ranking set to form a receptor conformation discrimination feature set;

[0011] S5: Identify the target receptor conformational combinations complementary to misclassified molecules based on the receptor conformational discrimination feature set, and determine the screening order according to the coverage order of the pocket state by the target receptor conformational combinations.

[0012] S6: Connect the candidate compound library to the target receptor conformation combination in the order of screening, and output the candidate compound screening results.

[0013] In a preferred embodiment, S1 specifically refers to:

[0014] The trajectory frames in the molecular dynamics simulation trajectory of the target protein are aligned with coordinates, and trajectory frames missing the coordinates of binding pocket residues are removed. The receptor conformation set is extracted according to the sampling interval.

[0015] The contact relationships of binding pocket residues for each receptor conformation are obtained based on the atomic distance relationships between binding pocket residues.

[0016] By associating the contact relationships of the pocket residues according to the order of the receptor conformation origin trajectory frames, a receptor conformation contact map is constructed.

[0017] In a preferred embodiment, S2 specifically refers to:

[0018] Read the receptor conformation contact map in the order of the receptor conformation source trajectory frame, compare the contact relationship of the binding pocket residues corresponding to adjacent receptor conformations, and obtain the number of newly added contacts and the number of lost contacts.

[0019] The source trajectory frame where the contact change direction reverses is determined based on the number of new contacts and the number of contacts disappearing, and the source trajectory frame where the contact change direction reverses is used as the pocket opening and closing transition node.

[0020] Receptor conformations that maintain continuous contact relationships between binding pocket residues are selected on both sides of the pocket opening-closing transition node to form candidate receptor conformations.

[0021] In a preferred embodiment, S3 specifically refers to:

[0022] The docking search region is determined by the contact relationship of the binding pocket residues corresponding to the candidate receptor conformation, and the conformation and charge state of the benchmark active molecule set and the decoy molecule set are standardized.

[0023] The standardized baseline active molecule set and decoy molecule set are docked to each candidate receptor conformation, respectively, preserving the binding posture that satisfies the docking search region constraints;

[0024] Record the molecular origin category, binding posture, and docking score corresponding to each candidate receptor conformation, and generate a set of binding postures and a set of scores for candidate receptor conformations according to the docking scores.

[0025] In a preferred embodiment, S4 specifically refers to:

[0026] The binding posture set of each candidate receptor conformation is read according to the molecular origin category. The burial depth is formed based on the distance from the heavy atom of the molecule in the binding posture to the boundary of the docking search region, and the difference in the burial depth distribution between the baseline active molecule set and the decoy molecule set is obtained.

[0027] An orientational contact relationship is formed based on the atomic distance between the heavy atoms of the molecule and the residues in the binding pocket;

[0028] By comparing the postureal contact relationships with the contact relationships of binding pocket residues, the differences in the number of key contact rearrangements formed by contact addition, contact disappearance, and contact replacement were statistically analyzed.

[0029] By correlating the differences in burial depth distribution, the differences in the number of key contact rearrangements, and the scoring and ranking sets, a set of receptor conformation discrimination features is formed.

[0030] In a preferred embodiment, S5 specifically refers to:

[0031] Based on the scoring and ranking set and the molecular origin category, molecules ranked lower than the baseline active molecule set and molecules ranked higher than the baseline active molecule set are labeled to form a misclassified molecule set for each candidate receptor conformation.

[0032] Based on the receptor conformation discrimination feature set, the overlap relationship of misclassified molecule sets, the direction of difference in burial depth distribution, and the direction of difference in the number of key contact rearrangements are compared, and candidate receptor conformations that are complementary to the misclassified molecule sets are selected as the target receptor conformation combination.

[0033] The screening order is determined according to the sequence of source trajectory frames on both sides of the pocket opening and closing transition node corresponding to the target receptor conformation combination.

[0034] In a preferred embodiment, S6 specifically refers to:

[0035] The candidate compound library is standardized in terms of conformation and charge state to form a standardized candidate compound library;

[0036] The standardized candidate compound library is sequentially docked to the target receptor conformation in the target receptor conformation combination according to the screening order, and the binding posture that satisfies the docking search region constraint corresponding to the target receptor conformation is retained.

[0037] Record candidate compounds, target receptor conformations, binding postures, and docking scores. Based on the docking scores and binding posture retention of the same candidate compounds in the target receptor conformation combinations, form the candidate compound screening results.

[0038] On the other hand, the present invention provides a drug virtual screening system based on molecular docking and kinetic simulation, comprising:

[0039] Receptor mapping module: Obtains the molecular dynamics simulation trajectory of the target protein and extracts the receptor conformation set, and establishes a receptor conformation contact map based on the contact relationship of binding pocket residues;

[0040] The transition node identification module identifies pocket-opening and closing transition nodes according to the receptor conformation contact map and selects candidate receptor conformations from both sides of the pocket-opening and closing transition nodes.

[0041] Molecular docking evaluation module: docks the baseline active molecule set and the decoy molecule set to the candidate receptor conformation, generating a set of binding postures and a set of scores and rankings for the candidate receptor conformation;

[0042] Conformation discrimination analysis module: For each candidate receptor conformation, calculate the differences in burial depth distribution and key contact rearrangement times in the binding posture set, and combine them with the scoring and ranking set to form a receptor conformation discrimination feature set;

[0043] Combination sorting determination module: Based on the receptor conformation discrimination feature set, identify the target receptor conformation combination complementary to the misclassified molecules, and determine the screening order according to the coverage order of the pocket state of the target receptor conformation combination.

[0044] Virtual screening output module: Connects the candidate compound library to the target receptor conformation combination in the screening order and outputs the candidate compound screening results.

[0045] The technical effects and advantages of the drug virtual screening strategy and system based on molecular docking and kinetic simulation proposed in this invention are as follows:

[0046] By acquiring the molecular dynamics simulation trajectory of target proteins and establishing receptor conformational contact maps, the dynamic changes in receptor conformation can be characterized by the contact relationships of binding pocket residues. By identifying pocket opening and closing transition nodes and selecting candidate receptor conformations on both sides of the nodes, the indiscriminate use of conformations obtained from molecular dynamics simulations for virtual screening is avoided. By evaluating the docking of candidate receptor conformations using a set of baseline active molecules and a set of decoy molecules, the binding posture set and scoring ranking set corresponding to different candidate receptor conformations can be obtained. By forming a set of receptor conformation discrimination features based on differences in burial depth distribution, differences in the number of key contact rearrangements, and the scoring ranking set, the screening discrimination ability of candidate receptor conformations can be quantitatively characterized. By identifying target receptor conformation combinations complementary to misclassified molecules and determining the screening order, receptor conformations in different pocket states can synergistically participate in candidate compound screening. Finally, the candidate compound library is sequentially docked to the target receptor conformation combinations, which helps to improve the reliability and specificity of candidate compound screening results. Attached Figure Description

[0047] Figure 1 This is a schematic diagram of a drug virtual screening strategy based on molecular docking and kinetic simulation according to the present invention.

[0048] Figure 2 This is a schematic diagram of the structure of a drug virtual screening system based on molecular docking and kinetic simulation according to the present invention. Detailed Implementation

[0049] 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. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0050] Example 1

[0051] Figure 1 This invention presents a virtual drug screening strategy based on molecular docking and kinetic simulation, which includes the following steps:

[0052] S1: Obtain the molecular dynamics simulation trajectory of the target protein and extract the receptor conformation set, and establish a receptor conformation contact map based on the contact relationship of binding pocket residues;

[0053] S2: Identify pocket opening-closing transition nodes according to the receptor conformational contact map, and select candidate receptor conformations from both sides of the pocket opening-closing transition nodes;

[0054] S3: Dock the baseline active molecule set and the decoy molecule set to the candidate receptor conformation to generate a set of binding postures and a set of scores and rankings for the candidate receptor conformation;

[0055] S4: For each candidate receptor conformation, calculate the differences in burial depth distribution and key contact rearrangement times in the binding posture set, and combine them with the scoring and ranking set to form a receptor conformation discrimination feature set;

[0056] S5: Identify the target receptor conformational combinations complementary to misclassified molecules based on the receptor conformational discrimination feature set, and determine the screening order according to the coverage order of the pocket state by the target receptor conformational combinations.

[0057] S6: Connect the candidate compound library to the target receptor conformation combination in the order of screening, and output the candidate compound screening results.

[0058] S1. Obtain the target protein molecular dynamics simulation trajectory and extract the receptor conformation set. Based on the contact relationships of binding pocket residues, establish a receptor conformational contact map, including:

[0059] The trajectory frames in the molecular dynamics simulation trajectory of the target protein are aligned with coordinates, and trajectory frames missing the coordinates of binding pocket residues are removed. The receptor conformation set is extracted according to the sampling interval.

[0060] Specifically, the target protein molecular dynamics simulation trajectory is a sequence of trajectory frames saved in chronological order of the simulation. Each trajectory frame contains atom name, residue name, residue number, chain identifier, and atomic three-dimensional coordinates. When the target protein molecular dynamics simulation trajectory uses nanometers as the unit of length, all atomic three-dimensional coordinates are multiplied by 10 to convert to angstroms; when the target protein molecular dynamics simulation trajectory already uses angstroms as the unit of length, the original atomic three-dimensional coordinates are retained. Atomic distances are uniformly calculated using Euclidean distance, which is obtained by taking the square root of the sum of the squares of the differences between the three-dimensional coordinates of two atoms along the three coordinate axes. Binding pocket residues are identified before trajectory frame processing. When the initial conformation of the target protein contains a reference ligand, the minimum distance between all heavy atoms of each amino acid residue and all heavy atoms of the reference ligand is calculated, and amino acid residues whose minimum distance is not greater than the pocket definition distance threshold are written into the binding pocket residue identifier table. Heavy atoms are non-hydrogen atoms. The pocket-bound distance threshold is determined by the minimum distance distribution of amino acid residues around the reference ligand. The minimum distance is counted by frequency analysis over a fixed bin width, which is determined according to the distance resolution requirements. For example, the bin width can be 0.1 Å. The center of the distance interval with the lowest frequency between the near-contact frequency peak and the far-non-contact frequency peak is then determined as the pocket-bound distance threshold. When there is no bimodal distribution of minimum distance, the pocket-bound distance threshold is determined by the maximum span of the reference ligand and the contact buffer distance. The maximum span of the reference ligand is the maximum distance between any two heavy atoms of the reference ligand. The contact buffer distance is used to cover the thermal fluctuations between the boundary residues of the binding pocket and the reference ligand. For example, the contact buffer distance can be 2.0 Å. The pocket-bound distance threshold can be, for example, 4.5 Å to 6.0 Å. When the initial conformation of the target protein does not have a reference ligand, the center of the active site is determined based on the geometric center of the heavy atoms of known functional residues, or the center of the cavity is determined based on the geometric center of the atoms on the inner surface of the cavity. Then, amino acid residues within a predetermined distance range around the center of the active site or the center of the cavity are written into a binding pocket residue identification table. The predetermined distance range is determined based on the cavity radius or the known active site scale; for example, the predetermined distance range can be 6.0 Å to 10.0 Å. The binding pocket residue identification table stores the residue name, residue number, and chain identifier for each binding pocket residue.

[0061] The coordinate alignment process begins with the determination of the coordinate alignment reference trajectory frame. From the target protein molecular dynamics simulation trajectory, the trajectory frame that is saved first in the sequence and simultaneously satisfies the following conditions: all binding pocket residues contain heavy atom coordinates, and the protein backbone atoms used for coordinate alignment are complete. Complete protein backbone atoms mean that each amino acid residue involved in the coordinate alignment simultaneously contains nitrogen, α-carbon, and carbonyl carbon atom coordinates. A one-to-one correspondence between the trajectory frame to be processed and the coordinate alignment reference trajectory frame is established according to atom name, residue name, residue number, and chain identifier. If the one-to-one correspondence of protein backbone atoms that can be used for alignment is less than 3 or if the spatial positions of the protein backbone atoms that can be used for alignment are collinear, the trajectory frame to be processed will not be included in the coordinate alignment process. After determining the one-to-one correspondence between protein backbone atoms, the geometric centers of the protein backbone atom set in the trajectory frame to be processed and the coordinate alignment reference trajectory frame are calculated separately. These two geometric centers are then aligned through translation. A covariance matrix is ​​constructed based on the one-to-one correspondence of the protein backbone atoms, and this covariance matrix is ​​decomposed to obtain a rigid body rotation matrix. Finally, the rigid body rotation matrix and translation vector are simultaneously applied to the three-dimensional coordinates of all atoms in the trajectory frame to be processed, completing the trajectory frame coordinate alignment. Coordinate alignment only changes the overall spatial position and orientation of the trajectory frame to be processed, without altering the relative distances between atoms within the trajectory frame.

[0062] Integrity checks are performed after trajectory frame coordinate alignment is complete. The heavy atom coordinates of each binding pocket residue in each trajectory frame are checked according to the binding pocket residue identifier table. If any binding pocket residue is missing a heavy atom coordinate used for atomic distance calculation, has a null heavy atom coordinate, has a non-numerical heavy atom coordinate, or its residue number or chain identifier is inconsistent with the binding pocket residue identifier table, the current trajectory frame is recorded as a trajectory frame missing binding pocket residue coordinates. Trajectory frames missing binding pocket residue coordinates are removed from the sampling objects, while retained trajectory frames continue to maintain their original sequential positions in the target protein molecular dynamics simulation trajectory. The number of retained trajectory frames is used to determine the sampling interval. If the number of retained trajectory frames is insufficient to support receptor conformation set extraction, the target protein molecular dynamics simulation trajectory is extended, and coordinate alignment and integrity checks are re-executed. Once the number of retained trajectory frames meets the requirements for receptor conformation set extraction, the sampling interval is determined.

[0063] The sampling interval is determined based on the redundancy of conformational changes of binding pocket residues in the retained trajectory frames. A minimum atomic distance matrix for binding pocket residues is established for each retained trajectory frame. The matrix elements of this matrix are the minimum distances among all heavy atom pairs of any two different binding pocket residues within the same retained trajectory frame. Then, using one trajectory frame as the minimum candidate sampling interval, a sequence of candidate sampling intervals is generated in an integer incrementing manner. The upper limit of the candidate sampling interval is determined based on statistical stability. The current candidate sampling interval is retained when the number of comparable trajectory frame pairs under each candidate sampling interval is not less than the lower limit of the number of trajectory frame pairs. The lower limit of the number of trajectory frame pairs is determined to be 10% to 20% of the total number of retained trajectory frames; for example, it can be 10 pairs. For each candidate sampling interval, the trajectory frame pairs corresponding to the interval are selected according to the original sequential positions of the retained trajectory frames. The absolute difference between the corresponding matrix elements of the minimum atomic distance matrices of the two binding pocket residues in the same trajectory frame pair is calculated. The average of the absolute differences of all matrix elements and all trajectory frame pairs is then calculated to obtain the average matrix difference corresponding to the current candidate sampling interval. The increment between the average differences of two adjacent matrices is calculated in ascending order of the candidate sampling intervals. The sampling stability threshold is used to identify when the average matrix difference enters a phase of gradual change. This threshold is determined based on the increment of the average matrix difference corresponding to the gradual change at the tail end of the candidate sampling interval sequence. This gradual change is defined as the last 10% to 30% of the candidate sampling interval sequence; for example, the last 20% could be considered. If the length of the candidate sampling interval sequence is insufficient to form a gradual change, the sampling stability threshold is the median of all increments of the average matrix difference. When the increments of the average matrix difference corresponding to the current candidate sampling interval and subsequent consecutive verification intervals are both no greater than the sampling stability threshold, the current candidate sampling interval is determined as the sampling interval. Continuous verification intervals are used to exclude random fluctuations; the length of a continuous verification interval can be, for example, two candidate sampling intervals. If no candidate sampling interval meets the conditions, the candidate sampling interval with the smallest increment of the average matrix difference is determined as the sampling interval.

[0064] Trajectory frames are sequentially extracted from the retained trajectory frames according to the sampling interval, and the three-dimensional coordinates of the target protein atoms and the binding pocket residue identification table corresponding to each extracted trajectory frame are written into the receptor conformation set. The order of the receptor conformation source trajectory frames is recorded according to the original order of the extracted trajectory frames in the target protein molecular dynamics simulation trajectory. Each receptor conformation in the receptor conformation set includes the three-dimensional coordinates of the target protein atoms, the binding pocket residue identification table, and the order of the receptor conformation source trajectory frames. The receptor conformation set is stored according to the order of the receptor conformation source trajectory frames.

[0065] The contact relationships of binding pocket residues for each receptor conformation are obtained based on the atomic distance relationships between binding pocket residues.

[0066] Specifically, a fixed list of binding pocket residue pairs is generated based on the binding pocket residue identifier table. This list is formed by pairwise combinations of different binding pocket residues. Each binding pocket residue pair in the list is arranged in the order of chain identifier, residue number, and residue name, and remains unchanged across all receptor conformations. Each receptor conformation in the receptor conformation set is traversed, and the minimum heavy atom distance is calculated for each binding pocket residue pair in the list. The minimum heavy atom distance is the minimum of the distances between all heavy atom pairs of two different binding pocket residues. The contact determination distance threshold is determined by the distribution of minimum heavy atom distances corresponding to all receptor conformations and all binding pocket residue pairs. The minimum heavy atom distance frequency is counted using a fixed bin width, determined according to distance resolution requirements. For example, a bin width of 0.1 Å is used. The center of the first distance interval between the short-distance and long-distance frequency peaks, where the frequency is simultaneously lower than both the previous and subsequent distance intervals, is defined as the contact determination distance threshold. If no distance interval exists where the minimum heavy atom distance distribution is simultaneously lower than both the previous and subsequent distance intervals, the contact determination distance threshold is determined by combining the sum of the van der Waals radii of the corresponding binding pocket residue heavy atoms with a distance buffer. The distance buffer covers fluctuations in minimum heavy atom distance caused by thermal motion and is determined by the standard deviation of all minimum heavy atom distance distributions on the short-distance side. For example, a distance buffer can be between 0.5 Å and 1.0 Å, and the contact determination distance threshold can be between 4.0 Å and 5.0 Å. When the minimum heavy atom distance is not greater than the contact determination distance threshold, the current binding pocket residue pair is considered to have contact in the current receptor conformation; when the minimum heavy atom distance is greater than the contact determination distance threshold, the current binding pocket residue pair is considered to have no contact in the current receptor conformation. The presence and absence of contacts for all binding pocket residue pairs in each receptor conformation together constitute the contact relationship of the binding pocket residues corresponding to the current receptor conformation.

[0067] By associating the contact relationships of the pocket residues according to the order of the receptor conformation origin trajectory frames, a receptor conformation contact map is constructed.

[0068] Specifically, the receptor conformation set is sorted according to the order of the receptor conformation source trajectory frames. Then, the contact relationships of the binding pocket residues corresponding to each receptor conformation are converted into contact state sequences according to the fixed order of the list of binding pocket residue pairs. In the contact state sequence, the presence of a contact is marked as 1, and the absence of a contact is marked as 0. All contact state sequences are written into a matrix row by row according to the order of the receptor conformation source trajectory frames, and the order of the list of binding pocket residue pairs is used as the column order, thereby establishing a receptor conformation contact map. In the receptor conformation contact map, each row corresponds to one receptor conformation source trajectory frame order, and each column corresponds to one binding pocket residue pair. The receptor conformation contact map simultaneously stores the receptor conformation source trajectory frame order and the contact relationships of the binding pocket residues.

[0069] S2. Identify pocket opening-closing transition nodes according to the receptor conformational contact map, and select candidate receptor conformations from both sides of the pocket opening-closing transition nodes, including:

[0070] Read the receptor conformation contact map in the order of the receptor conformation source trajectory frame, compare the contact relationship of the binding pocket residues corresponding to adjacent receptor conformations, and obtain the number of newly added contacts and the number of lost contacts.

[0071] Specifically, all row and column records in the receptor conformational contact map are read. The list of binding pocket residue pairs remains constant across all receptor conformations; therefore, the same column in the receptor conformational contact map represents the contact state changes of the same binding pocket residue pair across different rows. Adjacent rows are grouped into adjacent receptor conformational comparison pairs from front to back according to the receptor conformational origin trajectory frame order. The first row in each pair corresponds to the previous receptor conformation, and the second row corresponds to the next receptor conformation. The previous and next receptor conformations are consecutively arranged in the receptor conformational origin trajectory frame order.

[0072] For each adjacent receptor conformation comparison pair, the contact state between the preceding and following receptor conformations is compared column by column. If the value in the column corresponding to the preceding receptor conformation is 0 and the value in the column corresponding to the following receptor conformation is 1, the corresponding binding pocket residue pair is recorded as newly added contact. If the value in the column corresponding to the preceding receptor conformation is 1 and the value in the column corresponding to the following receptor conformation is 0, the corresponding binding pocket residue pair is recorded as lost contact. Binding pocket residue pairs with the same value in both columns are not included in the number of newly added and lost contacts. After comparing all columns of the current adjacent receptor conformation comparison pair, the number of newly added and lost contacts are counted separately. The number of newly added contacts is the number of binding pocket residue pairs recorded as newly added contacts in the current adjacent receptor conformation comparison pair, and the number of lost contacts is the number of binding pocket residue pairs recorded as lost contacts in the current adjacent receptor conformation comparison pair.

[0073] The source trajectory frame where the contact change direction reverses is determined based on the number of new contacts and the number of contacts disappearing, and the source trajectory frame where the contact change direction reverses is used as the pocket opening and closing transition node.

[0074] Specifically, to eliminate single-frame contact flipping caused by thermal fluctuations, orientation persistence verification is performed on binding pocket residue pairs marked as contact additions or loss in each adjacent receptor conformation comparison pair. The input data for orientation persistence verification consists of the state change records in the current adjacent receptor conformation comparison pair and the contact map row records from subsequent consecutive reads. The orientation persistence verification length is determined according to the resolution of the receptor conformation source trajectory frame order; for example, the orientation persistence verification length can be one or two adjacent receptor conformation comparison pairs. If a contact addition occurs in the current binding pocket residue pair within the current adjacent receptor conformation comparison pair, and the contact remains present within the orientation persistence verification length, the contact addition record for the current binding pocket residue pair is retained. If the current binding pocket residue pair reverts to a contact loss state within the orientation persistence verification length, the contact addition record for the current binding pocket residue pair is recorded as an instantaneous fluctuation and discarded. If a contact loss occurs in the current binding pocket residue pair within the current adjacent receptor conformational comparison pair, and the contact loss continues within the directional persistence verification length, the contact loss record for the current binding pocket residue pair is retained. If the current binding pocket residue pair recovers to a contact state within the directional persistence verification length, the contact loss record for the current binding pocket residue pair is recorded as an instantaneous fluctuation and discarded. After completing the directional persistence verification, the number of new contacts and the number of contact loss corresponding to each adjacent receptor conformational comparison pair are counted again.

[0075] When determining the direction of contact change based on the number of new contacts and the number of contacts lost, the contact change difference is calculated for each adjacent receptor conformation comparison pair. The contact change difference is the number of new contacts minus the number of contacts lost. When the contact change difference is greater than 0, the current adjacent receptor conformation comparison pair is recorded as the direction of increased contact. When the contact change difference is less than 0, the current adjacent receptor conformation comparison pair is recorded as the direction of decreased contact. When the contact change difference is equal to 0, the current adjacent receptor conformation comparison pair is recorded as the direction of contact equilibrium.

[0076] After obtaining the contact increase direction, contact decrease direction, and contact equilibrium direction, a direction sequence is established sequentially along the receptor conformation origin trajectory frames. Consecutive identical directions in the direction sequence are denoted as co-directional segments, while contact equilibrium directions form separate equilibrium segments. To identify when the contact change direction reverses, cases are identified where one co-directional segment, one equilibrium segment, and another co-directional segment are connected end-to-end, or where two co-directional segments are directly connected. A reversal of the contact change direction is defined as the preceding co-directional segment being a contact increase direction and the following co-directional segment being a contact decrease direction. Similarly, a reversal of the contact change direction is defined as the preceding co-directional segment being a contact decrease direction and the following co-directional segment being a contact increase direction. To avoid isolated direction changes consisting of only one adjacent receptor conformational comparison pair directly causing a reversal of the contact change direction, the number of adjacent receptor conformational comparison pairs contained in each of the preceding and following co-directional segments must not be less than a direction continuity threshold. The directional persistence threshold is determined based on the total number of receptor conformations. For example, the directional persistence threshold can be taken as the number of adjacent receptor conformation comparison pairs corresponding to 5% of the total number of receptor conformations. When the directional persistence threshold is less than 1, it is taken as 1.

[0077] When determining the source trajectory frame where the contact change direction reverses, for each interval where the contact change direction reverses, the number of new contacts and the number of contacts lost for all adjacent receptor conformation comparison pairs within the interval are read. If there is no equilibrium segment within the interval where the contact change direction reverses, the sequence of the source trajectory frames for the next receptor conformation of the last adjacent receptor conformation comparison pair in the previous same-direction segment is determined as the candidate node source trajectory frame sequence. Simultaneously, the sequence of the source trajectory frames for the previous receptor conformation of the first adjacent receptor conformation comparison pair in the next same-direction segment is determined as the candidate node source trajectory frame sequence. If the two candidate node source trajectory frame sequences are the same, the corresponding receptor conformation source trajectory frame sequence is determined as the pocket opening / closing transition node. If the two candidate node source trajectory frame sequences are different, the sum of the number of new contacts and the number of contacts lost for the adjacent receptor conformation comparison pairs corresponding to the two candidate node source trajectory frame sequences is compared, and the candidate node source trajectory frame sequence with the larger sum is determined as the pocket opening / closing transition node. When there is a balance segment in the interval where the contact change direction reverses, the sequence of all receptor conformation source trajectory frames covered by the balance segment is taken as the candidate node source trajectory frame sequence. The sum of the number of new contacts and the number of contacts lost for the adjacent receptor conformation comparison pairs on both sides of each candidate node source trajectory frame sequence is calculated. The candidate node source trajectory frame sequence with the largest sum is then determined as the pocket opening and closing transition node. When there are multiple candidate node source trajectory frame sequences with the same sum, the candidate node source trajectory frame sequence with the earlier receptor conformation source trajectory frame sequence is determined as the pocket opening and closing transition node.

[0078] Receptor conformations that maintain continuous contact relationships between binding pocket residues are selected on both sides of the pocket opening-closing transition node to form candidate receptor conformations.

[0079] Specifically, receptor conformations with continuously maintained binding pocket residue contact relationships are selected from both the area before and after the pocket opening / closing transition node. The area before the pocket opening / closing transition node starts with the receptor conformation immediately adjacent to the transition node and appearing earlier in the trajectory frame order. The area after the transition node starts with the receptor conformation immediately adjacent to the transition node and appearing later in the trajectory frame order. Using the preceding and following receptor conformations as centers, the process expands frame by frame away from the pocket opening / closing transition node, and the consistency of binding pocket residue contact relationships between adjacent receptor conformations is compared. The consistency comparison is calculated by counting the number of identical contact states between two adjacent receptor conformations across the entire list of binding pocket residue pairs, and then dividing this number by the total number of binding pocket residue pairs to obtain the contact consistency ratio. The contact consistency ratio ranges from 0 to 1; a higher ratio indicates a more stable binding pocket residue contact relationship between adjacent receptor conformations.

[0080] The continuity retention threshold corresponding to the contact consistency ratio is determined based on the contact consistency ratio distribution of all adjacent receptor conformations on both sides of the pocket opening / closing transition node. The contact consistency ratios of all adjacent receptor conformations before and after the pocket opening / closing transition node are statistically analyzed separately. These contact consistency ratios in both directions are then merged and sorted, and the values ​​at the 50% to 75% positions are used as candidate values ​​for the continuity retention threshold. Among the candidate values, the one that minimizes the fluctuation of the contact consistency ratio within both the preceding and following continuity retention intervals is determined as the continuity retention threshold. If a candidate value cannot simultaneously satisfy the requirement that both the preceding and following continuity retention intervals contain at least one receptor conformation, the median of all contact consistency ratios is determined as the continuity retention threshold. When the contact consistency ratio is not less than the continuity retention threshold, the corresponding receptor conformation remains within the continuity retention interval; when the contact consistency ratio is less than the continuity retention threshold, the corresponding receptor conformation stops expanding in the current direction. To avoid premature termination of the continuous hold interval due to a single abnormal adjacent receptor conformation, one isolated fluctuation is allowed on each side that is below the continuous hold threshold but subsequently rises above it again. Isolated fluctuations are not counted as termination conditions for the continuous hold interval. When the continuous hold threshold is lowered twice consecutively, the continuous hold interval terminates at the position where the continuous hold threshold is lowered for the first time.

[0081] After obtaining the continuous holding intervals before and after the pocket opening / closing transition node, candidate receptor conformations are selected from each continuous holding interval. During selection, the contact relationships of binding pocket residues corresponding to all receptor conformations within the same continuous holding interval are used as input data. The number of present contacts and the number of missing contacts are counted column by column. When the number of present contacts is greater than the number of missing contacts, the current column is recorded as the representative contact state on the same side. When the number of present contacts is less than the number of missing contacts, the current column is recorded as the representative contact state on the same side. When the number of present contacts equals the number of missing contacts, the contact state corresponding to the receptor conformation closer to the pocket opening / closing transition node in the current column is recorded as the representative contact state on the same side. The number of difference columns between each receptor conformation within the same continuous holding interval and the representative contact state on the same side is compared one by one. The receptor conformation with the smallest number of difference columns is determined as the current candidate receptor conformation. When multiple receptor conformations have the same number of difference columns, the receptor conformation whose source trajectory frame order is closer to the pocket opening / closing transition node is determined as the current candidate receptor conformation. The candidate receptor conformations formed by combining all candidate receptor conformations before and after the pocket opening-closing transition nodes together constitute the candidate receptor conformation. The candidate receptor conformation preserves the order of the receptor conformation origin trajectory frames, the contact relationships of binding pocket residues, and the correspondence between pocket opening-closing transition nodes.

[0082] S3. The baseline active molecule set and the decoy molecule set are docked to the candidate receptor conformation to generate a set of binding postures and a scoring and ranking set of the candidate receptor conformation, including:

[0083] The docking search region is determined by the contact relationship of the binding pocket residues corresponding to the candidate receptor conformation, and the conformation and charge state of the benchmark active molecule set and the decoy molecule set are standardized.

[0084] Specifically, when determining the docking search region based on the contact relationships of binding pocket residues corresponding to the candidate receptor conformation, the following information is retrieved from the candidate receptor conformation: the three-dimensional coordinates of target protein atoms, the binding pocket residue identification table, the binding pocket residue contact relationships, the order of the receptor conformation source trajectory frames, and the correspondence between pocket opening and closing transition nodes. In the binding pocket residue contact relationships, contacting binding pocket residue pairs reflect close residue combinations within the binding pocket of the candidate receptor conformation, while contactless binding pocket residue pairs reflect separated residue combinations within the binding pocket of the candidate receptor conformation. The docking search region is determined based on the spatial range of the binding pockets in the candidate receptor conformation, without altering the atomic three-dimensional coordinates of the candidate receptor conformation or the recorded binding pocket residue contact relationships. To ensure consistency between the docking search region and the pocket state of the candidate receptor conformation, the three-dimensional coordinates of all heavy atoms of binding pocket residues are extracted according to the binding pocket residue identification table, and then contacting binding pocket residue pairs are divided into contact connectivity groups according to the binding pocket residue contact relationships. The method for dividing contact connectivity groups is to merge any two pairs of binding pocket residues that share a common binding pocket residue and are both considered to be in contact into the same contact connectivity group, until no new pairs of binding pocket residues can be merged. Each contact connectivity group corresponds to a set of local pocket wall atoms.

[0085] For each contact connectivity group, the minimum and maximum coordinate values ​​of the local pocket inner wall atom set in the three coordinate axes are statistically analyzed to obtain the local enclosing interval. The same statistical analysis is repeated for all binding pocket residue heavy atoms to obtain the full pocket enclosing interval. The initial boundary of the docking search region is determined by superimposing the local enclosing interval and the full pocket enclosing interval. The superposition method is as follows: in each coordinate axis direction, the minimum coordinate value in all local enclosing intervals is compared with the minimum coordinate value in the full pocket enclosing interval, and the smaller value is taken as the initial lower boundary; the maximum coordinate value in all local enclosing intervals is compared with the maximum coordinate value in the full pocket enclosing interval, and the larger value is taken as the initial upper boundary. To avoid the initial boundary only covering the surface of the binding pocket residues without covering the ligand-accessible space, the initial lower and upper boundaries are extended. The boundary extension distance is determined based on the distance distribution from all binding pocket residue heavy atoms to the geometric center of the full pocket. The geometric center of the full pocket is the average coordinate of the three-dimensional coordinates of all binding pocket residue heavy atoms in the three coordinate axes. The distances from all heavy atoms of the binding pocket residues to the geometric center of the entire pocket are sorted from smallest to largest. Distances located at the 80% to 95% positions are selected as candidate boundary expansion distances. The proportion of the boundary covered by the expanded heavy atoms of all binding pocket residues is then used as a screening criterion to select the minimum boundary expansion distance that results in a coverage ratio of at least 0.95. If no distance value exists at the 80% to 95% positions, the median of all distances is used as the boundary expansion distance. After determining the boundary expansion distance, the boundary expansion distance is subtracted from the initial lower boundary and added to the initial upper boundary along each of the three coordinate axes to form the docking search region corresponding to the candidate receptor conformation. The docking search region corresponding to the candidate receptor conformation is represented by the lower and upper boundaries along the three coordinate axes, and a corresponding record is established with the source trajectory frame sequence of the candidate receptor conformation.

[0086] Before docking, both the baseline active molecule set and the decoy molecule set undergo conformational and charge state normalization. The baseline active molecule set is a collection of molecules with the target active label, while the decoy molecule set is a collection of molecules without the target active label that are comparable to the baseline active molecule set in terms of molecular size, number of hydrogen bond donors, number of hydrogen bond acceptors, number of rotatable bonds, or hydrophobicity range. Each molecular record in both the baseline active molecule set and the decoy molecule set includes at least the molecular identifier, molecular origin category, atomic connection relationships, atom type, bond type, and initial three-dimensional coordinates or two-dimensional connection information. Conformation normalization first performs a structural integrity check on each molecular record, including atomic valence state checks, bond order consistency checks, aromaticity representation unification, and multi-component record splitting. In the atomic valence state check, molecular records with a number of atomic connections incompatible with the allowed atomic valence states are marked as anomalous molecular records. In the bond order consistency check, molecular records where the same pair of atoms is assigned conflicting bond orders are marked as anomalous molecular records. In the aromaticity representation unification, atoms and bonds within the aromatic ring are uniformly converted into a consistent connection representation. In multi-component record splitting, for molecular records containing a main molecule portion, ions, or solvent molecules, only the main molecule portion with the highest number of heavy atoms and continuous covalent linkages is retained. Abnormal molecular records that can be corrected through atomic linkages are retained after correction; abnormal molecular records that cannot be corrected through atomic linkages are removed from the baseline active molecule set or decoy molecule set, while retaining the removal marker and molecule origin category.

[0087] Hydrogen atoms are added to the retained molecular records to generate three-dimensional conformations. The addition of hydrogen atoms is constrained by atomic valence equilibrium; a corresponding number of hydrogen atoms are added to each heavy atom with unsatisfied valence bonds. During three-dimensional conformation generation, rotatable single bonds, ring structures, and planar conjugated structures are identified, and torsion angles are assigned to each rotatable single bond. The torsion angle sampling step size is determined based on the number of rotatable single bonds; for example, 30 degrees can be used when the number of rotatable single bonds is no more than 4, and 60 degrees can be used when the number of rotatable single bonds is greater than 4. After generating candidate three-dimensional conformations according to the torsion angle sampling step size, collision checks are performed on each candidate three-dimensional conformation. During collision checks, if the atomic distance between any two non-bonded atoms is less than the van der Waals radius of the corresponding atom minus the collision tolerance, the current candidate three-dimensional conformation is recorded as a collision conformation and discarded. The collision tolerance is determined based on the atomic radius distribution of all retained molecular records, and can be, for example, between 0.2 Å and 0.5 Å. The retained candidate 3D conformations undergo geometric optimization. Geometric optimization reduces intramolecular tension by iteratively adjusting bond lengths, bond angles, and torsion angles. Geometric optimization terminates when the total energy change between two consecutive iterations does not exceed a predetermined convergence threshold. The predetermined convergence threshold is determined based on the total energy change distribution of the initial conformations of all retained molecules, and can be, for example, 1% to 5% of the median of all total energy changes. After geometric optimization, several 3D conformations with the lowest total energy and no redundancy are retained for each molecule record. 3D conformation redundancy is determined based on the deviation of all heavy atom coordinates between different 3D conformations of the same molecule record. If the deviation of heavy atom coordinates is lower than the conformation redundancy threshold, only the 3D conformation with the lower total energy is retained. The conformation redundancy threshold can be determined, for example, based on the values ​​of the top 20% of the deviation distribution of all 3D conformations.

[0088] Charge state normalization is performed after conformational normalization. Acidic, basic, and amphoteric functional groups are identified in each molecular record, and optional charge states are generated based on the target environment's acidity / base conditions. These target environment acidity / base conditions are determined based on experimental conditions related to the target protein or commonly used physiological conditions; for example, pH 6.5 to 8.0. Protonated and deprotonated states are generated for each ionizable functional group, and all ionizable functional group states are combined into optional charge states at the molecular level. When multiple optional charge states exist, the distribution ratio of each optional charge state under the target environment's acidity / base conditions is calculated, and optional charge states with a distribution ratio not lower than the state retention threshold are retained. The state retention threshold is determined based on state coverage requirements and the computational complexity of subsequent docking; for example, it can be 0.1 to 0.2. After charge state normalization, each retained molecular record forms one or more normalized molecular records. Each normalized molecular record stores the molecular identifier, molecular origin category, normalized 3D conformation, normalized atom type, and normalized charge state.

[0089] The standardized baseline active molecule set and decoy molecule set are docked to each candidate receptor conformation, respectively, preserving the binding posture that satisfies the docking search region constraints;

[0090] Specifically, when the standardized baseline active molecule set and decoy molecule set are docked to each candidate receptor conformation, the docking search region corresponding to the candidate receptor conformation is used as the spatial boundary, and the three-dimensional coordinates of the target protein atoms of the candidate receptor conformation are used as the receptor spatial constraints. Receptor preparation is performed on the candidate receptor conformation before docking. Receptor preparation includes adding polar hydrogen atoms, standardizing atom types, correcting the charge state of ionizable side chains in the binding pocket residues, and labeling immobile atoms. Immobile atoms are atoms other than side chain atoms near the binding pocket residues; their coordinates remain unchanged during docking. For each standardized molecule record, the standardized three-dimensional conformation is placed into the docking search region one by one, and pose sampling is performed along the translational, rotational, and rotatable single-bond degrees of freedom. During pose sampling, if any heavy atom in the normalized molecular record exceeds the boundary of the docking search region, the current pose is marked as invalid. If the atomic distance between any heavy atom in the normalized molecular record and any heavy atom in the candidate acceptor conformation is less than the collision determination distance, the current pose is marked as invalid. If no effective contact is formed between any heavy atom in the normalized molecular record and any heavy atom in the binding pocket residue, the current pose is marked as invalid. The collision determination distance is determined based on the van der Waals radii of the two atoms involved in the determination and the collision tolerance. Effective contact is determined based on the minimum atomic distance between the heavy atom in the normalized molecular record and the heavy atom in the binding pocket residue. Effective contact is defined as when the minimum atomic distance is not greater than the effective contact distance threshold, which can be, for example, 4.0 Å to 5.0 Å. Invalid poses are removed from subsequent scoring, and the retained poses are considered binding pose candidates.

[0091] Docking scoring is performed on each candidate binding posture. The docking score is determined by a combination of spatial matching score, intermolecular interaction score, and conformational deviation score. The spatial matching score is determined based on the degree of space filling of the candidate binding posture within the docking search region. The intermolecular interaction score is determined based on hydrogen bonding, hydrophobic contact, electrostatic attraction, and electrostatic repulsion between the candidate binding posture and the candidate acceptor conformation. The conformational deviation score is determined based on the degree of torsional deviation of the candidate binding posture relative to the normalized 3D conformation. To adapt the docking scoring to the differentiation requirements of the baseline active molecule set and the decoy molecule set, pre-test molecular records are extracted from both sets for pre-docking. The separation degree of the baseline active molecule set and the decoy molecule set under different weight combinations is then compared. Finally, the weight combination that results in the baseline active molecule set ranking higher than the decoy molecule set is selected as the official docking scoring weight. The number of pre-test molecular records can be, for example, 10% to 20% of the number of each set. For each normalized molecular record, several non-redundant candidate binding postures with optimal docking scores are retained within each candidate acceptor conformation. The redundancy determination of binding posture candidates is based on the deviation of all heavy atom coordinates between different binding posture candidates recorded by the same standardized molecule. When the deviation of all heavy atom coordinates is lower than the posture redundancy threshold, only the binding posture candidate with the better docking score is retained. The posture redundancy threshold can be determined, for example, based on the values ​​of the top 20% of the deviations of all binding posture candidates' heavy atom coordinates.

[0092] Record the molecular origin category, binding posture, and docking score corresponding to each candidate receptor conformation, and generate a set of binding postures and a set of scores for candidate receptor conformations according to the docking scores.

[0093] Specifically, the candidate receptor conformation origin trajectory frame sequence is used as the primary index, and the molecular identifier and normalized charge state in the normalized molecular record are used as the secondary index. Each record at least stores the candidate receptor conformation origin trajectory frame sequence, molecular identifier, molecular origin category, normalized charge state, three-dimensional coordinates of the binding posture atom, a list of contact residues corresponding to the binding posture, and a docking score. The list of contact residues corresponding to the binding posture is calculated based on the atomic distance between the three-dimensional coordinates of the binding posture atom and the heavy atoms of the binding pocket residues. After docking all normalized molecular records with all candidate receptor conformations, all retained records are first summarized according to the candidate receptor conformation origin trajectory frame sequence, and then arranged in descending order of docking score within each candidate receptor conformation. Two records with the same docking score are arranged in descending order of the number of binding pocket residues in the contact residue list corresponding to the binding posture; two records with the same number of binding pocket residues are arranged in ascending order of the total energy corresponding to the normalized three-dimensional conformation. After sorting, all retained binding posture records within the same candidate receptor conformation form the binding posture set of the candidate receptor conformation, and the sorting results corresponding to all retained records within the same candidate receptor conformation form the scoring and sorting set of the candidate receptor conformation. Both the binding posture set and the scoring and sorting set of the candidate receptor conformation are associated with the source trajectory frame order, molecular source category, molecular identifier, and normalized charge state of the candidate receptor conformation.

[0094] S4. For each candidate receptor conformation, calculate the differences in burial depth distribution and key contact rearrangement frequency in the binding posture set, and combine them with the scoring and ranking set to form a receptor conformation discrimination feature set, including:

[0095] The binding posture set of each candidate receptor conformation is read according to the molecular origin category. The burial depth is formed based on the distance from the heavy atom of the molecule in the binding posture to the boundary of the docking search region, and the difference in the burial depth distribution between the baseline active molecule set and the decoy molecule set is obtained.

[0096] Specifically, when reading the binding posture set and scoring ranking set for each candidate receptor conformation according to the molecular origin category, the binding posture set, scoring ranking set, docking search region boundary, binding pocket residue identifier table, and binding pocket residue contact relationship corresponding to each candidate receptor conformation are read sequentially according to the trajectory frame order of the candidate receptor conformation origin. Abnormal binding posture records are removed from the binding posture set corresponding to the current candidate receptor conformation. The binding posture set corresponding to the current candidate receptor conformation is split into a baseline active molecule set binding posture subset and a decoy molecule set binding posture subset according to the molecular origin category, while maintaining the order of the remaining records in the scoring ranking set corresponding to the current candidate receptor conformation.

[0097] When determining the burial depth based on the distance from the molecular heavy atom to the boundary of the docking search region in the binding posture, for each binding posture record corresponding to the current candidate acceptor conformation, the three-dimensional coordinates of the molecular heavy atom and the boundary of the docking search region corresponding to the current candidate acceptor conformation are read. The boundary of the docking search region is jointly represented by the lower boundary along the first coordinate axis, the upper boundary along the first coordinate axis, the lower boundary along the second coordinate axis, the upper boundary along the second coordinate axis, the lower boundary along the third coordinate axis, and the upper boundary along the third coordinate axis. For the same molecular heavy atom, the distances from the molecular heavy atom to the lower boundary along the first coordinate axis, the upper boundary along the first coordinate axis, the lower boundary along the second coordinate axis, the upper boundary along the second coordinate axis, the lower boundary along the third coordinate axis, and the upper boundary along the third coordinate axis are calculated respectively, and the minimum value among the six distances is determined as the boundary distance of the current molecular heavy atom. When the current molecular heavy atom is located inside the boundary of the docking search region, the boundary distance is non-negative; when the current molecular heavy atom is located outside the boundary of the docking search region, the boundary distance is negative. If any binding posture record contains a molecular heavy atom with a negative boundary distance, the current binding posture record is recorded as an out-of-bounds binding posture record. If no out-of-bounds binding posture record is found among all retained binding posture records, the burial depth calculation begins; if an out-of-bounds binding posture record exists, it is simultaneously removed from the binding posture set and scoring ranking set corresponding to the current candidate acceptor conformation.

[0098] For each non-boundary binding posture record, the boundary distances corresponding to all molecular heavy atoms are sorted from smallest to largest. To reduce the influence of flexible atoms at the molecular ends on the overall burial depth characterization, the boundary distances between the 25% and 75% positions after sorting are used as the central boundary distance sequence, and the average value of the central boundary distance sequence is determined as the burial depth of the current binding posture record. When the central boundary distance sequence is insufficient to form an average value, the average value of the boundary distances corresponding to all molecular heavy atoms is determined as the burial depth of the current binding posture record. The unit of burial depth is consistent with the unit of atomic three-dimensional coordinates. Since the atomic three-dimensional coordinates have been unified to angstroms in the previous steps, the unit of burial depth is angstroms. After calculating the burial depth of all non-boundary binding posture records, burial depth distributions are established for the binding posture subsets of the baseline active molecule set and the binding posture subsets of the decoy molecule set, respectively. The burial depth distribution includes at least the median burial depth, the interquartile range of burial depth, and the high burial proportion. The threshold for determining the high burial proportion is determined based on the overall burial depth distribution of all non-boundary binding posture records within the current candidate receptor conformation. The embedment depths of all non-boundary binding postures are sorted from smallest to largest, and the embedment depth at the 75th percentile after sorting is defined as the high embedment threshold. The number of binding posture records with embedment depths not less than the high embedment threshold in the current molecular source category is divided by the total number of all non-boundary binding posture records in the current molecular source category to obtain the high embedment ratio corresponding to the current molecular source category. Following the direction of subtracting the corresponding value of the decoy molecule set from the corresponding value of the baseline active molecule set, the median embedment depth difference, the interquartile range embedment depth difference, and the high embedment ratio difference are calculated respectively, which together constitute the embedment depth distribution difference corresponding to the current candidate receptor conformation.

[0099] An orientational contact relationship is formed based on the atomic distance between the heavy atoms of the molecule and the residues in the binding pocket;

[0100] Specifically, the three-dimensional coordinates of the target protein atoms, the binding pocket residue identification table, and the contact relationships of the binding pocket residues corresponding to the current candidate receptor conformation are read. Each binding pocket residue in the binding pocket residue identification table is uniquely identified by its residue name, residue number, and chain identifier. For each non-cross-boundary binding posture record corresponding to the current candidate receptor conformation, the minimum atomic distance between the molecular heavy atoms and all heavy atoms of each binding pocket residue is calculated. The posture contact determination distance threshold is determined based on the overall distribution of minimum atomic distances corresponding to all non-cross-boundary binding posture records within the current candidate receptor conformation. The frequency of all minimum atomic distances is counted according to a bin width of 0.1 Å, and the first distance interval between the short-distance frequency peak and the long-distance frequency peak that is simultaneously lower than the frequency of the adjacent preceding interval and the adjacent following interval is found in the minimum atomic distance frequency distribution. The center of the current distance interval is then determined as the posture contact determination distance threshold. When there is no distance interval in the minimum atomic distance frequency distribution, the sum of the van der Waals radius of the binding pocket residue heavy atoms and the van der Waals radius of the molecular heavy atoms, plus the distance buffer, is determined as the posture contact determination distance threshold. The distance buffer is determined based on the standard deviation of the distribution of all minimum atomic distances on the shorter distance side, and can be, for example, between 0.5 Å and 1.0 Å. When the minimum atomic distance between a molecular heavy atom and all heavy atoms of the current binding pocket residue is not greater than the posture contact determination distance threshold, the current binding pocket residue is recorded as having posture contact in the current binding posture record; when the minimum atomic distance between a molecular heavy atom and all heavy atoms of the current binding pocket residue is greater than the posture contact determination distance threshold, the current binding pocket residue is recorded as having posture contact absent in the current binding posture record. The posture contact presence and absence states corresponding to all binding pocket residues in the current binding posture record together form the posture contact relationship corresponding to the current binding posture record.

[0101] By comparing the postureal contact relationships with the contact relationships of binding pocket residues, the differences in the number of key contact rearrangements formed by contact addition, contact disappearance, and contact replacement were statistically analyzed.

[0102] Specifically, based on the list of binding pocket residue pairs, the posture contact relationship corresponding to each binding posture record is converted into a posture pairing contact state sequence. A given binding pocket residue pair in the list corresponds to two binding pocket residues. If both binding pocket residues corresponding to the current binding pocket residue pair are recorded as posture contacts in the current binding posture record, the current binding pocket residue pair is recorded as a posture pairing contact in the posture pairing contact state sequence; if either of the two binding pocket residues corresponding to the current binding pocket residue pair is recorded as a posture contact absence, the current binding pocket residue pair is recorded as a posture pairing contact absence in the posture pairing contact state sequence. The posture pairing contact state sequence is compared column-by-column with the binding pocket residue contact relationship corresponding to the current candidate receptor conformation, in the order of the binding pocket residue pair list. If the current column is recorded as a contact absence in the binding pocket residue contact relationship and as a posture pairing contact presence in the posture pairing contact state sequence, the corresponding binding pocket residue pair in the current column is recorded as a contact addition. When the current column is marked as having contact in the binding pocket residue contact relationship and as having missing contact in the posture pairing contact state sequence, the binding pocket residue pair corresponding to the current column is marked as having lost contact.

[0103] Contact substitution is determined according to the residue adjacency substitution rule. Each binding pocket residue pair corresponding to contact disappearance in each binding posture record is checked one by one. For the binding pocket residue pair corresponding to the current contact disappearance, the first and second binding pocket residues in the current binding pocket residue pair are read. Then, in the contact addition set of the same binding posture record, a contact addition corresponding to a binding pocket residue pair that shares one identical binding pocket residue with the first or second binding pocket residue is searched. After finding a contact addition corresponding to a contact addition pair that shares one identical binding pocket residue, it is then checked whether the two binding pocket residues that were not shared before and after sharing are recorded as contacts in the binding pocket residue contact relationship corresponding to the current candidate receptor conformation. When the two binding pocket residues that were not shared before and after sharing are recorded as contacts in the binding pocket residue contact relationship corresponding to the current candidate receptor conformation, the current contact disappearance and the current contact addition are recorded together as one contact substitution. After the current contact disappearance and the current contact addition are recorded as contact substitutions, they are not counted again in the contact disappearance and contact addition counts. After completing the matching of all contact disappearances and contact additions, the remaining number of new contacts, the remaining number of contact disappearances, and the number of contact replacements are counted for the current bonding posture record. The remaining number of new contacts, the remaining number of contact disappearances, and the number of contact replacements together constitute the key contact rearrangement count for the current bonding posture record, where the key contact rearrangement count is the sum of the remaining number of new contacts, the remaining number of contact disappearances, and the number of contact replacements.

[0104] A sorting cutoff is performed on the scoring and ranking set corresponding to the current candidate receptor conformation. The sorting cutoff position is determined based on the distribution of docking score differences between adjacent ranking positions within the current candidate receptor conformation. The docking score differences between two adjacent records are read from best to worst according to the scoring and ranking set. The position where three consecutive docking score differences are lower than the median of all adjacent docking score differences is determined as the starting position of the ranking platform. The ranking position preceding the starting position of the ranking platform is then determined as the sorting cutoff position. If no starting position of the ranking platform exists, the ranking positions corresponding to the first 30% of the scoring and ranking sets are determined as the sorting cutoff positions. Binding posture records whose ranking positions are not inferior to the sorting cutoff positions are retained in the binding posture subsets of the baseline active molecule set and the binding posture subset of the decoy molecule set, respectively. The distribution of critical contact rearrangement times corresponding to the retained records is then statistically analyzed. The critical contact rearrangement time distribution includes at least the median critical contact rearrangement time, the interquartile range of critical contact rearrangement times, and the high rearrangement ratio. The threshold for determining the high rearrangement ratio is determined based on the overall distribution of critical contact rearrangement times for all retained binding posture records after sorting cutoff. After sorting and truncating, all binding posture records retaining the key contact rearrangement counts are sorted from smallest to largest, and the value at the 75th percentile after sorting is determined as the high rearrangement threshold. The number of binding posture records in the current molecule source category with key contact rearrangement counts not less than the high rearrangement threshold is divided by the number of binding posture records in the current molecule source category that retain all binding posture records after sorting and truncating, yielding the high rearrangement ratio corresponding to the current molecule source category. Following the direction of subtracting the corresponding value of the decoy molecule set from the corresponding value of the baseline active molecule set, the median difference of key contact rearrangement counts, the interquartile range difference of key contact rearrangement counts, and the high rearrangement ratio difference are calculated respectively, which together constitute the key contact rearrangement count difference corresponding to the current candidate receptor conformation.

[0105] By correlating the differences in burial depth distribution, the differences in the number of key contact rearrangements, and the scoring and ranking sets, a set of receptor conformation discrimination features is formed.

[0106] Specifically, using the order of candidate receptor conformation source trajectory frames as the index unit, the differences in burial depth distribution, key contact rearrangement frequency, and scoring and ranking sets are simultaneously read within the current candidate receptor conformation, and the scoring and ranking sets are converted into ranking distinction indicators. The ranking distinction indicators include at least the proportion of the baseline active molecule set in the top 10% of the ranking positions, the first occurrence position of the baseline active molecule set, the first occurrence position of the decoy molecule set, and the number of intersections of molecule source categories. The proportion of the baseline active molecule set in the top 10% of the ranking positions is the ratio obtained by dividing the number of records in the top 10% of the scoring and ranking sets corresponding to the current candidate receptor conformation whose molecule source category belongs to the baseline active molecule set by the total number of records in the top 10% of the ranking positions. The first occurrence position of the baseline active molecule set is the ranking position corresponding to the first record in the scoring and ranking set corresponding to the current candidate receptor conformation whose molecule source category belongs to the baseline active molecule set. The first occurrence position of the decoy molecule set is the ranking position corresponding to the first record in the scoring and ranking set corresponding to the current candidate receptor conformation whose molecule source category belongs to the decoy molecule set. The molecular origin category crossover count is the number of times the molecular origin category of two adjacent records changes when the scoring and ranking set corresponding to the current candidate receptor conformation is read from best to worst according to the ranking position. The differences in median burial depth, interquartile range of burial depth, high burial ratio, median difference in key contact rearrangement count, interquartile range of key contact rearrangement count, high rearrangement ratio, the proportion of the baseline active molecule set in the top 10% ranking positions, the first occurrence position of the baseline active molecule set, the first occurrence position of the decoy molecule set, and the molecular origin category crossover count are written into the receptor conformation discrimination feature record corresponding to the current candidate receptor conformation in a fixed field order. All receptor conformation discrimination feature records corresponding to all candidate receptor conformations together form the receptor conformation discrimination feature set. Each receptor conformation discrimination feature record in the receptor conformation discrimination feature set maintains a one-to-one correspondence with the candidate receptor conformation origin trajectory frame order.

[0107] S5. Identify target receptor conformational combinations complementary to misclassified molecules based on the receptor conformational discrimination feature set, and determine the screening order according to the coverage order of the pocket states by the target receptor conformational combinations, including:

[0108] Based on the scoring and ranking set and the molecular origin category, molecules ranked lower than the baseline active molecule set and molecules ranked higher than the baseline active molecule set are labeled to form a misclassified molecule set for each candidate receptor conformation.

[0109] Specifically, based on the order of the source trajectory frames of the candidate receptor conformation, the scoring and ranking set, receptor conformation discrimination feature set record, pocket opening and closing transition node correspondence, and source trajectory frame order information on both sides of the pocket opening and closing transition node are read one by one for each candidate receptor conformation. Each ranking record in the scoring and ranking set includes at least the candidate receptor conformation source trajectory frame order, molecular identifier, molecular source category, normalized charge state, ranking position, and docking score. The molecular source category only includes the baseline active molecule set and the decoy molecule set. An integrity check is performed on the scoring and ranking set corresponding to the current candidate receptor conformation. Ranking records with missing molecular identifiers, missing molecular source categories, missing normalized charge states, missing ranking positions, missing docking scores, non-numerical ranking positions, non-numerical docking scores, or multiple ranking positions for the same molecular identifier under the same normalized charge state are recorded as abnormal ranking records. Abnormal ranking records are removed from the scoring and ranking set corresponding to the current candidate receptor conformation. After removing abnormal sorting records, the remaining sorting records are rearranged from best to worst according to their sorting position, and then split into a baseline active molecule set sorting subset and a decoy molecule set sorting subset according to the molecule source category.

[0110] A sorting inversion count is established within the current candidate receptor conformation. For each sorting record in the subset of the baseline active molecule set, the number of decoy molecule set sorting records preceding the current sorting record is counted to obtain the forward decoy count corresponding to the current baseline active molecule set sorting record. For each sorting record in the subset of the decoy molecule set, the number of baseline active molecule set sorting records following the current sorting record is counted to obtain the backward activity count corresponding to the current decoy molecule set sorting record. When the current sorting record and adjacent sorting records have the same docking score, the sorting records with the same docking score are treated as a whole segment, and forward decoy counting or backward activity counting begins at the first sorting record with a different docking score after the end of the segment, to avoid misclassification marker shift caused by the same docking score. After completing the forward decoy counting and backward activity counting, all forward decoy counts and all backward activity counts within the current candidate receptor conformation are sorted from smallest to largest, and the value at the 75th percentile after sorting is determined as the misclassification counting threshold. The sorted records of the baseline active molecule set with forward decoy counts not less than the misclassification count threshold are denoted as the baseline active molecule set misclassification records. The sorted records of the decoy molecule set with backward activity counts not less than the misclassification count threshold are denoted as the decoy molecule set misclassification records. When the number of non-zero records of forward decoy counts or backward activity counts within the current candidate receptor conformation is insufficient to determine 75% of the position values, the median of all non-zero counts is determined as the misclassification count threshold; when none of the non-zero counts exist, the misclassification count threshold is recorded as 0.

[0111] After determining the misclassification count threshold, the misclassified molecule labels are constrained based on the distribution of sorting positions. The cumulative proportions of the baseline active molecule set and the decoy molecule set are calculated for all sorting records corresponding to the current candidate receptor conformation. The cumulative proportion of the baseline active molecule set is the number of records whose molecular origin category belongs to the baseline active molecule set at the current sorting position and all preceding sorting records, divided by the total number of records at the current sorting position and all preceding sorting records. The cumulative proportion of the decoy molecule set is the number of records whose molecular origin category belongs to the decoy molecule set at the current sorting position and all preceding sorting records, divided by the total number of records at the current sorting position and all preceding sorting records. The sorting position where the cumulative proportions of the baseline active molecule set and the cumulative proportions of the decoy molecule set first reverse in size is determined as the sorting boundary position. When multiple adjacent sorting positions simultaneously satisfy the size reversal, the sorting position preceding the adjacent sorting position with the largest docking score difference is determined as the sorting boundary position. When the sorting position of a misclassified molecule record in the baseline active molecule set is after the sorting boundary position and the forward decoy count is not less than the misclassification count threshold, the corresponding molecule is recorded as a misclassified molecule in the baseline active molecule set. When a misclassified molecule in the decoy molecule set is located before the sorting boundary and its backward activity count is not less than the misclassification count threshold, the corresponding molecule identifier is recorded as a misclassified molecule in the decoy molecule set. The misclassified molecules in the baseline active molecule set and the misclassified molecules in the decoy molecule set together constitute the misclassified molecule set corresponding to the current candidate receptor conformation. When there is no labeled molecule identifier for the current candidate receptor conformation, the misclassified molecule set corresponding to the current candidate receptor conformation is recorded as an empty set, while maintaining a one-to-one correspondence between the source trajectory frame order of the current candidate receptor conformation and the empty set.

[0112] Based on the receptor conformation discrimination feature set, the overlap relationship of misclassified molecule sets, the direction of difference in burial depth distribution, and the direction of difference in the number of key contact rearrangements are compared, and candidate receptor conformations that are complementary to the misclassified molecule sets are selected as the target receptor conformation combination.

[0113] Specifically, the sets of misclassified molecules and the set of receptor conformation discrimination features corresponding to all candidate receptor conformations are read in the order of the trajectory frames from which the candidate receptor conformations originate. For any two candidate receptor conformations, the number of intersections, unions, intersections, and unions of misclassified molecules in the baseline active molecule set, respectively, are calculated. When the union of misclassified molecules in the baseline active molecule set is greater than 0, the misclassification overlap ratio of the baseline active molecule set is obtained by dividing the number of intersections by the number of unions; when the union of misclassified molecules in the baseline active molecule set is equal to 0, the misclassification overlap ratio of the baseline active molecule set is recorded as 0. When the union of misclassified molecules in the decoy molecule set is greater than 0, the misclassification overlap ratio of the decoy molecule set is obtained by dividing the number of intersections by the number of unions; when the union of misclassified molecules in the decoy molecule set is equal to 0, the misclassification overlap ratio of the decoy molecule set is recorded as 0. The weights for the overlap ratios of the baseline active molecule set and the decoy molecule set are determined based on the proportions of misclassified molecules in the baseline active molecule set and the decoy molecule set, respectively. These weighted sums are then used to obtain the total overlap ratio of the misclassified molecule sets. Subtracting the total overlap ratio from 1 yields the complementarity ratio of the misclassified molecule sets. A higher complementarity ratio indicates greater complementarity between the misclassified molecule sets corresponding to the two candidate receptor conformations.

[0114] The direction of differences in burial depth distribution and the direction of differences in key contact rearrangement times are determined by the directional consistency of the corresponding differences recorded in the receptor conformation discrimination feature set. For each candidate receptor conformation, the median difference in burial depth, the interquartile range difference in burial depth, and the high burial proportion difference are read. Then, the median of the absolute values ​​of the corresponding differences for all candidate receptor conformations is calculated to obtain the median difference buffer threshold, the interquartile range difference buffer threshold, and the high burial proportion difference buffer threshold. The current difference whose absolute value is less than 20% of the corresponding buffer threshold is recorded as 0. The above 20% is the directional buffer proportion, which can be, for example, 10% to 30%. After directional buffering, if at least two of the three burial depth differences are greater than 0, the burial depth distribution difference direction corresponding to the current candidate receptor conformation is recorded as positive; if at least two of the three burial depth differences are less than 0, the burial depth distribution difference direction corresponding to the current candidate receptor conformation is recorded as negative; otherwise, the burial depth distribution difference direction corresponding to the current candidate receptor conformation is recorded as the equilibrium direction. The difference direction of key contact rearrangement times is determined in the same way. Directional buffering is performed on the median difference, interquartile range difference, and high rearrangement ratio difference of key contact rearrangement times. After directional buffering, if at least two of the three key contact rearrangement times differences are greater than 0, the key contact rearrangement times difference direction corresponding to the current candidate receptor conformation is recorded as positive; if at least two of the three key contact rearrangement times differences are less than 0, the key contact rearrangement times difference direction corresponding to the current candidate receptor conformation is recorded as negative; otherwise, the key contact rearrangement times difference direction corresponding to the current candidate receptor conformation is recorded as the equilibrium direction.

[0115] After obtaining the complementarity ratio of the misclassified molecule sets, the direction of difference in burial depth distribution, and the direction of difference in the number of key contact rearrangements, candidate receptor conformations complementary to the misclassified molecule sets are selected as target receptor conformation combinations. Pairwise combinations are performed on all candidate receptor conformations, and a combination screening value is calculated for each candidate receptor conformation combination. The combination screening value is jointly determined by the complementarity ratio of the misclassified molecule sets, the matching results of the direction of difference in burial depth distribution, and the matching results of the direction of difference in the number of key contact rearrangements. When the burial depth distribution direction of the two candidate receptor conformations in the current candidate receptor conformation combination is opposite to the direction of difference in burial depth distribution, the matching result of the direction of difference in burial depth distribution is marked as satisfied; when the direction of difference in the number of key contact rearrangements of the two candidate receptor conformations is opposite to the direction of difference in burial depth distribution, the matching result of the direction of difference in the number of key contact rearrangements is marked as satisfied. When both the matching results of the direction of difference in burial depth distribution and the matching results of the direction of difference in the number of key contact rearrangements are satisfied, the current candidate receptor conformation combination is retained as a first-level candidate combination. If only one of the following conditions is met—matching the burial depth distribution direction and the key contact rearrangement frequency direction—and the complementary ratio of the misclassified molecule set is not lower than the complementary threshold, the current candidate receptor conformation combination is retained as a secondary combination candidate. The complementary threshold is determined based on the complementary ratio distribution of the misclassified molecule sets obtained from pairwise comparisons of all candidate receptor conformations; for example, the complementary threshold can be the median of the complementary ratios of all misclassified molecule sets. When primary combination candidates exist, the complementary ratios of the misclassified molecule sets of all primary combination candidates are compared, and the primary combination candidate with the largest complementary ratio is determined as the target receptor conformation combination. When multiple primary combination candidates have the same complementary ratio of their misclassified molecule sets, the sum of the proportions of the baseline active molecule sets in the top 10% of the ranking positions corresponding to each candidate receptor conformation in the current primary combination candidate is compared, and the primary combination candidate with the larger sum of the proportions of the baseline active molecule sets in the top 10% of the ranking positions is determined as the target receptor conformation combination. If the sum of the proportions of the baseline active molecules in the top 10% of the ranking positions is still the same, compare the sum of the number of cross-references of the molecular source categories for each candidate receptor conformation in the current primary candidate combination, and determine the primary candidate combination with the smaller sum of the number of cross-references of the molecular source categories as the target receptor conformation combination. If there are no primary candidate combinations but there are secondary candidate combinations, determine the target receptor conformation combination from the secondary candidate combinations according to the same rules. If there are no primary or secondary candidate combinations, increase the number of candidate receptor conformations in descending order of the complementarity ratio of the misclassified molecule sets to form a candidate receptor conformation combination containing 3 or more candidate receptor conformations. It is required that any two candidate receptor conformations in the current candidate conformation combination satisfy at least one difference direction matching result. Then, compare the average complementarity ratio of the misclassified molecule sets, the sum of the proportions of the baseline active molecules in the top 10% of the ranking positions, and the sum of the number of cross-references of the molecular source categories in turn to finally determine the target receptor conformation combination.

[0116] The screening order is determined according to the sequence of source trajectory frames on both sides of the pocket opening-closing transition node corresponding to the target receptor conformation combination.

[0117] Specifically, the process involves reading the pocket opening / closing transition node correspondence, the order of source trajectory frames preceding and following each candidate receptor conformation within the target receptor conformation combination. Each candidate receptor conformation originates from either the area preceding or following a pocket opening / closing transition node. To ensure the screening order reflects the continuous coverage of pocket state changes, candidate receptor conformations within the target receptor conformation combination are grouped by pocket opening / closing transition nodes, resulting in node groups. Candidate receptor conformations preceding and following the same pocket opening / closing transition node are grouped into the same node group. Different node groups are sorted in ascending order based on the position of the pocket opening / closing transition node within the receptor conformation source trajectory frame sequence. Within each node group, the screening order is determined according to the source trajectory frame sequence. When a node group contains both candidate receptor conformations before and after a pocket-opening / closing transition node, the order of the source trajectory frames before the pocket-opening / closing transition node is read, followed by the order of the source trajectory frames after the pocket-opening / closing transition node. The candidate receptor conformations within the node group are then arranged in ascending order according to the source trajectory frame order. When a node group contains only one candidate receptor conformation, the current candidate receptor conformation is used as the unique selection position within the current node group. After sorting all node groups internally, the selection positions within all node groups are concatenated according to the order between node groups to form the selection order corresponding to the target receptor conformation combination. The selection order must at least preserve the source trajectory frame order of each candidate receptor conformation in the target receptor conformation combination, the pocket-opening / closing transition node correspondence, the position within the node group, and the overall sorting position.

[0118] S6. Connect the candidate compound library to the target receptor conformation combinations in the order of screening, and output the candidate compound screening results, including:

[0119] The candidate compound library is standardized in terms of conformation and charge state to form a standardized candidate compound library;

[0120] Specifically, all candidate compound records in the candidate compound library are read. Each candidate compound record includes at least a candidate compound identifier, atomic connection relationships, atom types, bond types, and initial three-dimensional coordinates or two-dimensional connection information. A structural integrity check is performed on all candidate compound records, including atomic valence state checks, bond order consistency checks, aromaticity representation unification, and multi-component record splitting. In the atomic valence state check, candidate compound records with atomic connection numbers incompatible with allowed atomic valence states are marked as anomalous candidate compound records. In the bond order consistency check, candidate compound records where the same pair of atoms is assigned conflicting bond orders are marked as anomalous candidate compound records. In the aromaticity representation unification, atoms and bonds within the aromatic ring are uniformly converted into a consistent connection representation. In the multi-component record splitting, for candidate compound records containing a main molecule portion, ion pairs, or solvent molecules, only the main molecule portion with the most heavy atoms and continuous covalent connections is retained. If an anomalous candidate compound record can be corrected through atomic connection relationships, it is retained after correction; if an anomalous candidate compound record cannot be corrected through atomic connection relationships, it is removed from the candidate compound library, and the correspondence between the removal marker and the candidate compound identifier is retained.

[0121] Hydrogen atoms are added to the records of retained candidate compounds to generate three-dimensional conformations. The addition of hydrogen atoms is constrained by atomic valence equilibrium; a corresponding number of hydrogen atoms are added to each heavy atom with an unsatisfied valence bond. During three-dimensional conformation generation, rotatable single bonds, ring structures, and planar conjugated structures are identified, and torsion angles are assigned to each rotatable single bond. The torsion angle sampling step size is determined based on the number of rotatable single bonds. When the number of rotatable single bonds is no more than 4, the torsion angle sampling step size can be, for example, 30 degrees; when the number of rotatable single bonds is greater than 4, the torsion angle sampling step size can be, for example, 60 degrees. After generating candidate three-dimensional conformations according to the torsion angle sampling step size, collision checks are performed on each candidate three-dimensional conformation. During collision checks, if the atomic distance between any two non-bonded atoms is less than the van der Waals radius of the corresponding atom minus the collision tolerance, the current candidate three-dimensional conformation is recorded as a collision conformation and discarded. The collision tolerance is determined based on the atomic radius distribution of all retained candidate compound records; the collision tolerance can be, for example, 0.2 Å to 0.5 Å. The retained candidate 3D conformations undergo geometric optimization. Geometric optimization reduces intramolecular tension by iteratively adjusting bond lengths, bond angles, and torsion angles. Geometric optimization terminates when the total energy change between two consecutive iterations does not exceed a predetermined convergence threshold. The predetermined convergence threshold is determined based on the total energy change distribution of the initial conformations of all retained candidate compounds, and can be, for example, 1% to 5% of the median of all total energy changes. After geometric optimization, several non-redundant 3D conformations with the lowest total energy are retained for each candidate compound. 3D conformation redundancy is determined based on the deviation of all heavy atom coordinates between different 3D conformations of the same candidate compound. If the deviation of all heavy atom coordinates is lower than the conformation redundancy threshold, only the 3D conformation with the lower total energy is retained. The conformation redundancy threshold is determined based on the distribution of all 3D conformation deviations, and can be, for example, the value at the top 20% of the distribution.

[0122] Charge state normalization is performed after conformational normalization. Acidic, basic, and amphoteric functional groups are identified in each candidate compound record, and optional charge states are generated based on the target environment's acid-base conditions. These target environment acid-base conditions are consistent with those used in the conformational and charge state normalization of the baseline active molecule set and decoy molecule set; for example, pH 6.5 to 8.0. Protonated and deprotonated states are generated for each ionizable functional group, and all ionizable functional group states are combined to form optional charge states at the candidate compound level. When multiple optional charge states exist, the distribution ratio of each optional charge state under the target environment's acid-base conditions is calculated, and optional charge states with a distribution ratio not lower than the state retention threshold are retained. This state retention threshold is consistent with the state retention threshold used in the conformational and charge state normalization of the baseline active molecule set and decoy molecule set; for example, 0.1 to 0.2. After charge state normalization, each retained candidate compound record forms one or more normalized candidate compound records. Each normalized candidate compound record stores the candidate compound identifier, normalized 3D conformation, normalized atom types, and normalized charge states. All normalized candidate compound records together form a normalized candidate compound library.

[0123] The standardized candidate compound library is sequentially docked to the target receptor conformation in the target receptor conformation combination according to the screening order, and the binding posture that satisfies the docking search region constraint corresponding to the target receptor conformation is retained.

[0124] Specifically, the process reads the target receptor conformation combination, screening order, the 3D coordinates of the target protein atoms corresponding to the target receptor conformation, and the docking search region corresponding to the target receptor conformation. The screening order at least saves the source trajectory frame order of the target receptor conformation, the correspondence of pocket opening and closing transition nodes, the position within the node group, and the overall sorting position. Each target receptor conformation is read sequentially from front to back according to the overall sorting position, and all standardized candidate compound records in the standardized candidate compound library are docked to the current target receptor conformation one by one. Before docking, receptor preparation is performed on the current target receptor conformation. Receptor preparation includes adding polar hydrogen atoms, unifying the atom type, correcting the charge state of ionizable side chains in the binding pocket residues, and marking immobile atoms. Immobile atoms are atoms other than side chain atoms near the binding pocket residues, and the coordinates of immobile atoms remain unchanged during the docking process of the current target receptor conformation. For each standardized candidate compound record, the standardized 3D conformation is placed into the docking search region corresponding to the current target receptor conformation, and the pose is sampled along the translational degrees of freedom, rotational degrees of freedom, and rotatable single bond degrees of freedom. During posture sampling, if any heavy atom in the standardized candidate compound record exceeds the boundary of the docking search region corresponding to the current target receptor conformation, the current posture is recorded as invalid. If the atomic distance between any heavy atom in the standardized candidate compound record and any heavy atom in the current target receptor conformation is less than the collision determination distance, the current posture is recorded as invalid. If no effective contact is formed between any heavy atom in the standardized candidate compound record and any heavy atom in the binding pocket residues, the current posture is recorded as invalid. The collision determination distance is determined based on the van der Waals radii of the two atoms involved in the determination and the collision tolerance. Effective contact is determined based on the minimum atomic distance between the heavy atom in the standardized candidate compound record and the heavy atom in the binding pocket residues; if the minimum atomic distance is not greater than the effective contact distance threshold, it is recorded as effective contact. The effective contact distance threshold is consistent with the effective contact distance threshold used when docking the baseline active molecule set and the decoy molecule set to the candidate receptor conformation; for example, the effective contact distance threshold can be 4.0 Å to 5.0 Å. Invalid postures are removed from the scoring process, and the retained postures are used as candidate binding postures under the current target receptor conformation.

[0125] Docking scores are applied to each candidate binding posture. The docking score is determined by a combination of spatial matching score, intermolecular interaction score, and conformational deviation score. The spatial matching score is determined based on the degree of space filling in the docking search region corresponding to the current target receptor conformation. The intermolecular interaction score is determined based on hydrogen bonding, hydrophobic contact, electrostatic attraction, and electrostatic repulsion between the candidate binding posture and the current target receptor conformation. The conformational deviation score is determined based on the degree of torsional deviation of the candidate binding posture relative to the normalized 3D conformation. To ensure comparability of docking scores between different target receptor conformations, the calculation rules for spatial matching score, intermolecular interaction score, and conformational deviation score are consistent with those used when docking the baseline active molecule set and decoy molecule set to candidate receptor conformations. The docking score weights are also consistent with those used when docking the baseline active molecule set and decoy molecule set to candidate receptor conformations, and the weights are not readjusted. For each normalized candidate compound, several non-redundant binding posture candidates with optimal docking scores are recorded in each target receptor conformation. The redundancy determination of binding posture candidates is based on the deviation of all heavy atom coordinates between different binding posture candidates recorded by the same standardized candidate compound. When the deviation of all heavy atom coordinates is lower than the posture redundancy threshold, only the binding posture candidate with the better docking score is retained. The posture redundancy threshold is consistent with the posture redundancy threshold used when docking the baseline active molecule set and the decoy molecule set to the candidate receptor conformation. After completing the docking of all standardized candidate compound records for the current target receptor conformation, all retained binding posture candidates corresponding to the current target receptor conformation are written into the docking result record of the current target receptor conformation, and the next target receptor conformation is read in the screening order until all target receptor conformations in the target receptor conformation combination are docked.

[0126] Record candidate compounds, target receptor conformations, binding postures, and docking scores. Based on the docking scores and binding posture retention of the same candidate compounds in the target receptor conformational combinations, form the candidate compound screening results.

[0127] Specifically, the candidate compound identifier is used as the primary index, the target acceptor conformation source trajectory frame order as the secondary index, and the normalized charge state as the tertiary index. Each record stores at least the candidate compound identifier, the target acceptor conformation source trajectory frame order, the pocket opening / closing transition node correspondence, the normalized charge state, the three-dimensional coordinates of the binding posture atoms, the contact residue list corresponding to the binding posture, the docking score, and the overall sorting position. The contact residue list corresponding to the binding posture is calculated based on the atomic distance between the three-dimensional coordinates of the binding posture atoms and the heavy atoms of the binding pocket residues. To avoid duplicate records for the same candidate compound under the same target acceptor conformation, deduplication is performed on all records corresponding to the same candidate compound identifier, the same target acceptor conformation source trajectory frame order, and the same normalized charge state. During deduplication, records with better docking scores are retained first; among records with the same docking scores, records with more binding pocket residues in the contact residue list corresponding to the binding posture are retained first; among records with the same number of binding pocket residues, records with better conformation deviation scores corresponding to the three-dimensional coordinates of the binding posture atoms are retained first. After deduplication of all records, a candidate compound docking record set is formed.

[0128] When forming candidate compound screening results based on the docking scores and binding posture retention of the same candidate compounds in target receptor conformation combinations, the docking record sets of candidate compounds are first merged according to the candidate compound identifier to obtain the target receptor conformation retention record set for each candidate compound. The number of records in the target receptor conformation retention record set is recorded as the binding posture retention count of the current candidate compound. A binding posture retention count greater than 0 indicates that the current candidate compound retains a binding posture that satisfies the docking search region constraint corresponding to the current target receptor conformation in at least one target receptor conformation. All docking scores in the target receptor conformation retention record set corresponding to the current candidate compound are sorted from best to worst to obtain the docking score ranking within the target receptor conformation combination of the current candidate compound. The first docking score in the docking score ranking within the target receptor conformation combination is recorded as the first docking score. When the number of records in the docking score ranking within the target receptor conformation combination is not less than 2, the average of the first 2 docking scores is recorded as the double docking score; when the number of records in the docking score ranking within the target receptor conformation combination is only 1, the first docking score is also recorded as the double docking score. The binding posture retention rate is obtained by dividing the number of times the current candidate compound retains its binding posture by the total number of target receptor conformations in the target receptor conformation combination. To balance local optimal binding ability and cross-target receptor conformation retention ability, all candidate compounds are ranked from best to worst based on their first docking score. Then, among candidate compounds with the same first docking score, they are ranked from highest to lowest based on their binding posture retention rate. Candidate compounds with the same binding posture retention rate are further ranked from best to worst based on their two-position docking score. Candidate compounds with the same first docking score, binding posture retention rate, and two-position docking score are ranked from first to last according to the overall ranking position corresponding to the earliest record in the target receptor conformation source trajectory frame sequence. After ranking, the candidate compound identifier, first docking score, two-position docking score, number of binding posture retentions, binding posture retention rate, target receptor conformation source trajectory frame sequence list, and best binding posture record for each candidate compound are written into the candidate compound screening results.

[0129] Example 2

[0130] The difference between Embodiment 2 and Embodiment 1 is that this embodiment introduces a drug virtual screening system based on molecular docking and kinetic simulation.

[0131] Figure 2 A schematic diagram of a drug virtual screening system based on molecular docking and kinetic simulation is provided. The drug virtual screening system based on molecular docking and kinetic simulation includes:

[0132] Receptor mapping module: Obtains the molecular dynamics simulation trajectory of the target protein and extracts the receptor conformation set, and establishes a receptor conformation contact map based on the contact relationship of binding pocket residues;

[0133] The transition node identification module identifies pocket-opening and closing transition nodes according to the receptor conformation contact map and selects candidate receptor conformations from both sides of the pocket-opening and closing transition nodes.

[0134] Molecular docking evaluation module: docks the baseline active molecule set and the decoy molecule set to the candidate receptor conformation, generating a set of binding postures and a set of scores and rankings for the candidate receptor conformation;

[0135] Conformation discrimination analysis module: For each candidate receptor conformation, calculate the differences in burial depth distribution and key contact rearrangement times in the binding posture set, and combine them with the scoring and ranking set to form a receptor conformation discrimination feature set;

[0136] Combination sorting determination module: Based on the receptor conformation discrimination feature set, identify the target receptor conformation combination complementary to the misclassified molecules, and determine the screening order according to the coverage order of the pocket state of the target receptor conformation combination.

[0137] Virtual screening output module: Connects the candidate compound library to the target receptor conformation combination in the screening order and outputs the candidate compound screening results.

[0138] The above embodiments are only used to illustrate the technical methods of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical methods of the present invention without departing from the spirit and scope of the technical methods of the present invention.

Claims

1. A drug virtual screening strategy based on molecular docking and dynamics simulation, characterized in that, Includes the following steps: S1: Obtain the molecular dynamics simulation trajectory of the target protein and extract the receptor conformation set, and establish a receptor conformation contact map based on the contact relationship of binding pocket residues; S2: Identify pocket opening-closing transition nodes according to the receptor conformational contact map, and select candidate receptor conformations from both sides of the pocket opening-closing transition nodes; S3: Dock the baseline active molecule set and the decoy molecule set to the candidate receptor conformation to generate a set of binding postures and a set of scores and rankings for the candidate receptor conformation; S4: For each candidate receptor conformation, calculate the differences in burial depth distribution and key contact rearrangement times in the binding posture set, and combine them with the scoring and ranking set to form a receptor conformation discrimination feature set; S5: Identify the target receptor conformational combinations complementary to misclassified molecules based on the receptor conformational discrimination feature set, and determine the screening order according to the coverage order of the pocket state by the target receptor conformational combinations. S6: Connect the candidate compound library to the target receptor conformation combination in the order of screening, and output the candidate compound screening results. 2.The drug virtual screening strategy based on molecular docking and dynamics simulation according to claim 1, wherein, S1, specifically: The trajectory frames in the molecular dynamics simulation trajectory of the target protein are aligned with coordinates, and trajectory frames missing the coordinates of binding pocket residues are removed. The receptor conformation set is extracted according to the sampling interval. The contact relationships of binding pocket residues for each receptor conformation are obtained based on the atomic distance relationships between binding pocket residues. By associating the contact relationships of the pocket residues according to the order of the receptor conformation origin trajectory frames, a receptor conformation contact map is constructed.

3. The drug virtual screening strategy based on molecular docking and dynamics simulation according to claim 2, characterized in that, S2, specifically: Read the receptor conformation contact map in the order of the receptor conformation source trajectory frame, compare the contact relationship of the binding pocket residues corresponding to adjacent receptor conformations, and obtain the number of newly added contacts and the number of lost contacts. The source trajectory frame where the contact change direction reverses is determined based on the number of new contacts and the number of contacts disappearing, and the source trajectory frame where the contact change direction reverses is used as the pocket opening and closing transition node. Receptor conformations that maintain continuous contact relationships between binding pocket residues are selected on both sides of the pocket opening-closing transition node to form candidate receptor conformations.

4. The drug virtual screening strategy based on molecular docking and dynamics simulation according to claim 3, characterized in that, S3, specifically: The docking search region is determined by the contact relationship of the binding pocket residues corresponding to the candidate receptor conformation, and the conformation and charge state of the benchmark active molecule set and the decoy molecule set are standardized. The standardized baseline active molecule set and decoy molecule set are docked to each candidate receptor conformation, respectively, preserving the binding posture that satisfies the docking search region constraints; Record the molecular origin category, binding posture, and docking score corresponding to each candidate receptor conformation, and generate a set of binding postures and a set of scores for candidate receptor conformations according to the docking scores.

5. The drug virtual screening strategy based on molecular docking and dynamics simulation according to claim 4, characterized in that, S4, specifically: The binding posture set of each candidate receptor conformation is read according to the molecular origin category. The burial depth is formed based on the distance from the heavy atom of the molecule in the binding posture to the boundary of the docking search region, and the difference in burial depth distribution between the baseline active molecule set and the decoy molecule set is obtained. An orientational contact relationship is formed based on the atomic distance between the heavy atoms of the molecule and the residues in the binding pocket; By comparing the postureal contact relationships with the contact relationships of binding pocket residues, the differences in the number of key contact rearrangements formed by contact addition, contact disappearance, and contact replacement were statistically analyzed. By correlating the differences in burial depth distribution, the differences in the number of key contact rearrangements, and the scoring and ranking sets, a set of receptor conformation discrimination features is formed.

6. The drug virtual screening strategy based on molecular docking and dynamics simulation according to claim 5, characterized in that, S5, specifically: Based on the scoring and ranking set and the molecular origin category, molecules ranked lower than the baseline active molecule set and molecules ranked higher than the baseline active molecule set are labeled to form a misclassified molecule set for each candidate receptor conformation. Based on the receptor conformation discrimination feature set, the overlap relationship of misclassified molecule sets, the direction of difference in burial depth distribution, and the direction of difference in the number of key contact rearrangements are compared, and candidate receptor conformations that are complementary to the misclassified molecule sets are selected as the target receptor conformation combination. The screening order is determined according to the sequence of source trajectory frames on both sides of the pocket opening and closing transition node corresponding to the target receptor conformation combination.

7. The drug virtual screening strategy based on molecular docking and dynamics simulation according to claim 6, characterized in that, S6, specifically: The candidate compound library is standardized in terms of conformation and charge state to form a standardized candidate compound library; The standardized candidate compound library is sequentially docked to the target receptor conformation in the target receptor conformation combination according to the screening order, and the binding posture that satisfies the docking search region constraint corresponding to the target receptor conformation is retained. Record candidate compounds, target receptor conformations, binding postures, and docking scores. Based on the docking scores and binding posture retention of the same candidate compounds in the target receptor conformation combinations, form the candidate compound screening results.

8. A system for virtual screening of drugs based on molecular docking and dynamics simulation, for implementing a strategy for virtual screening of drugs based on molecular docking and dynamics simulation according to any one of claims 1 to 7, characterized in that, include: Receptor mapping module: Obtains the molecular dynamics simulation trajectory of the target protein and extracts the receptor conformation set, and establishes a receptor conformation contact map based on the contact relationship of binding pocket residues; The transition node identification module identifies pocket-opening and closing transition nodes according to the receptor conformation contact map and selects candidate receptor conformations from both sides of the pocket-opening and closing transition nodes. Molecular docking evaluation module: docks the baseline active molecule set and the decoy molecule set to the candidate receptor conformation, generating a set of binding postures and a set of scores and rankings for the candidate receptor conformation; Conformation discrimination analysis module: For each candidate receptor conformation, calculate the differences in burial depth distribution and key contact rearrangement times in the binding posture set, and combine them with the scoring and ranking set to form a receptor conformation discrimination feature set; Combination sorting determination module: Based on the receptor conformation discrimination feature set, identify the target receptor conformation combination complementary to the misclassified molecules, and determine the screening order according to the coverage order of the pocket state of the target receptor conformation combination. Virtual screening output module: Connects the candidate compound library to the target receptor conformation combination in the screening order and outputs the candidate compound screening results.

Images