Molecular docking affinity prediction method and system based on transfer learning, and medium

Abstract

The application discloses a molecular docking affinity prediction method and system based on transfer learning, and a medium, relates to the technical field of bioinformatics and artificial intelligence, and inputs a predicted protein-ligand combination into a prediction model to output a binding affinity index; the training process of the prediction model is as follows: a pre-training data set and a fine-tuning data set are constructed; the prediction model is constructed, including a structure encoder and an affinity predictor connected in sequence; in the pre-training stage, the prediction model is trained end to end by using the pre-training data set and pre-training labels; in the fine-tuning stage, the fine-tuning data set is used as input, the structure encoder parameters are fixed, and the affinity predictor parameters are fine-tuned with the prediction of the binding affinity index as the target; the method improves the prediction accuracy, generalization ability and multi-index comprehensive evaluation level of the prediction model in a specific target field.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics and artificial intelligence, and in particular to a method, system, and medium for predicting molecular docking affinity based on transfer learning. Background Technology

[0002] The interaction networks between proteins and small molecules in biological systems are extremely complex; human cells contain approximately 2 × 10⁻⁶ such networks. 4 These proteins can theoretically form up to 10 8 Potential interactions of magnitude. However, traditional prediction tasks rely on datasets only 10... 3 Up to 10 4 The limited number of records, coupled with their small size and incompleteness, severely restricts the generalization ability and prediction accuracy of models. Currently, the field also faces core bottlenecks such as difficulty in addressing the "activity cliff" phenomenon and the reliance on a single prediction indicator. Furthermore, traditional methods face a dilemma when dealing with molecular interactions: while traditional software based on physical force fields has a robust theoretical framework, it is inherently flawed due to approximations in fundamental physics; while machine learning models, due to a lack of training data, struggle to fully capture the complex characteristics of molecular interactions.

[0003] Existing solutions can be broadly categorized into three types. The first type is traditional molecular docking software, such as AutoDock Vina. These solutions primarily rely on physical force fields for theoretical support, but their computational accuracy is limited by the fundamental physical approximations of interatomic interactions, making it difficult to completely eliminate inherent errors. The second type consists of conventional machine learning prediction models, which predict affinity by directly training on known data. However, due to insufficient experimental data, these models exhibit weak generalization ability and limited accuracy when facing specific target domains. The third type is traditional sequence alignment schemes based on algorithms such as BLAST. These schemes typically rely on direct matching of amino acid letters to find homologous proteins. Their limitation lies in the tendency to miss distant homologous phenomena, meaning they cannot effectively identify key samples with large sequence differences but consistent physicochemical property distribution frequencies and similar binding pockets.

[0004] The main drawbacks of existing technologies lie in the extremely limited size and incompleteness of the datasets, making it difficult for machine learning models to fully capture the complex characteristics of molecular interactions. This limits their generalization ability and prediction accuracy, especially when facing the activity cliff phenomenon and the simultaneous prediction of multi-dimensional affinity indices. Furthermore, traditional docking software is limited by approximate models of physical force fields, resulting in inherent errors. Traditional sequence alignment algorithms, based on letter matching, are prone to missing distant homologous sequences with consistent physicochemical property distribution patterns. Summary of the Invention

[0005] Based on the technical problems existing in the background technology, this invention proposes a molecular docking affinity prediction method, system and medium based on transfer learning, which improves the prediction accuracy, generalization ability and multi-index comprehensive evaluation level of the prediction model in specific target fields.

[0006] The molecular docking affinity prediction method proposed in this invention, based on transfer learning, inputs the protein-ligand combination to be predicted into the prediction model and outputs a binding affinity index. The training process of the prediction model is as follows: Protein samples are selected from a pre-defined general set, multiple ligands are matched for each protein sample and molecular docking is performed to construct a pre-training dataset, and energy term parameters and docking scores are calculated as pre-training labels. The target protein sequence is sequentially encoded and subjected to quantum Fourier transform. The extracted quantum frequency domain fingerprint is compared with the samples in a preset general set to calculate the fingerprint fidelity. Homologous samples are then selected to construct a fine-tuning dataset. Construct a predictive model, which includes a structural encoder and an affinity predictor connected in sequence; Pre-training phase: The prediction model is trained end-to-end using the pre-training dataset and pre-training labels; Fine-tuning phase: Using the fine-tuning dataset as input, with the structure encoder parameters fixed, the affinity predictor parameters are fine-tuned with the goal of predicting the binding affinity index.

[0007] Furthermore, the construction of the pre-trained dataset includes: Protein structure identifiers were screened from the general PDBbind database to create a protein sample library; For each protein, multiple ligands were matched from the ChEMBL database for molecular docking. The energy parameters, including van der Waals attraction, van der Waals repulsion, solvation energy, and hydrogen bond energy, as well as the docking score, were calculated using molecular docking software and used as pre-training labels.

[0008] Furthermore, the sequential encoding of the target protein sequence specifically involves: The entire amino acid chain of the target protein sequence is considered as a continuous biological signal flow; Each residue of the target protein sequence is encoded using 8-dimensional physicochemical properties, and an adaptive resampling technique is used to uniformly map the encoded feature vector along the sequence length dimension. At each location point, a resampled physicochemical signal matrix is ​​obtained. This physicochemical signal matrix carries the global physicochemical information of the target protein sequence. The 8-dimensional physicochemical properties include hydrophobicity, molecular weight, isoelectric point, number of hydrogen bond donors, number of hydrogen bond acceptors, polarity, side chain flexibility, and aromaticity. This represents the number of qubits used in the quantum Fourier transform.

[0009] Furthermore, the quantum Fourier transform includes: Based on the dimensions of the physicochemical signal matrix Dynamic allocation One quantum bit resource; The global physicochemical information carried in the physicochemical signal matrix is ​​loaded into a quantum superposition state using amplitude encoding; The quantum superposition state is transformed from the spatial domain to the frequency domain by quantum Fourier transform, and the quantum frequency domain fingerprint of the target protein sequence is extracted.

[0010] Furthermore, before inputting the pre-trained dataset or fine-tuned dataset into the prediction model, a fusion graph needs to be constructed based on the samples. The specific construction process is as follows; A protein pocket diagram is constructed using the Cα atoms of protein binding pocket residues as residue nodes and the spatial distance between residues less than a preset threshold as edges. Small molecule graphs are constructed using ligand atoms as nodes and the chemical bonds between ligand atoms as edges. The protein pocket diagram is merged with the small molecule diagram, and cross-molecular connection edges are established between protein residue nodes and ligand atom nodes based on spatial proximity to obtain the fusion diagram.

[0011] Furthermore, the binding affinity indices in the fine-tuning stage include one or more of the dissociation constant, inhibition constant, and half-maximal inhibitory concentration, and the fine-tuning employs multi-task learning to achieve simultaneous prediction of multiple indices.

[0012] Furthermore, the construction of the fine-tuning dataset also includes: when the target protein has a clear domain division, the fine-tuning dataset is directly constructed using a domain-specific database.

[0013] Furthermore, the energy term parameters and docking fractions are calculated by AutoDock Vina.

[0014] A computer system includes a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the method described above.

[0015] A computer-readable storage medium storing a plurality of computer programs for being invoked by a processor and executing the method described above.

[0016] The advantages of the molecular docking affinity prediction method, system and medium based on transfer learning provided by this invention are as follows: By constructing a two-stage technical system of general knowledge foundation and domain characteristic calibration through transfer learning, it aims to solve the problem of high data dependence by using large-scale pre-training, and introduce quantum Fourier transform to achieve efficient homology retrieval and domain fine-tuning, thereby systematically improving the prediction accuracy, generalization ability and multi-index comprehensive evaluation level of the model in a specific target domain. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the structure of the present invention; Figure 2 A comparison chart of predicted RMSE under different training strategies. Detailed Implementation

[0018] The technical solution of the present invention will now be described in detail through specific embodiments. Many specific details are set forth in the following description to provide a thorough understanding of the invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0019] like Figure 1 and Figure 2 As shown, the molecular docking affinity prediction method based on transfer learning proposed in this invention inputs the protein-ligand combination to be predicted into the prediction model and outputs the binding affinity index. The training process of the prediction model is as follows: Protein samples are selected from a pre-defined general set, multiple ligands are matched for each protein sample and molecular docking is performed to construct a pre-training dataset, and energy term parameters and docking scores are calculated as pre-training labels. The target protein sequence is sequentially encoded and subjected to quantum Fourier transform. The extracted quantum frequency domain fingerprint is compared with the samples in a preset general set to calculate the fingerprint fidelity. Homologous samples are then selected to construct a fine-tuning dataset. Construct a predictive model, which includes a structural encoder and an affinity predictor connected in sequence; Pre-training phase: The prediction model is trained end-to-end using the pre-training dataset and pre-training labels; Fine-tuning phase: Using the fine-tuning dataset as input, with the structure encoder parameters fixed, the affinity predictor parameters are fine-tuned with the goal of predicting the binding affinity index.

[0020] The purpose of this embodiment is to construct a two-stage technical system of general knowledge foundation and domain characteristic calibration through transfer learning. It aims to solve the problem of high data dependence by using large-scale pre-training and introduce quantum Fourier transform to achieve efficient homology retrieval and domain fine-tuning, thereby systematically improving the prediction accuracy, generalization ability and multi-index comprehensive evaluation level of the prediction model in a specific target domain.

[0021] The aforementioned technical aspects are not isolated improvements, but rather an interconnected organic whole: the construction of a fine-tuning dataset enhanced by quantum Fourier transform ensures the learning quality during the fine-tuning stage; large-scale pre-training ensures that the structure encoder learns reliable general physical laws; the multi-level feature representation of the fusion graph enriches the feature vectors extracted by the structure encoder; the combination of these three aspects ultimately enables the prediction model to achieve prediction accuracy and generalization ability in small-sample, specific target domains that traditional solutions cannot match.

[0022] I. Data Preparation: Constructing the pre-trained dataset and the fine-tuning dataset.

[0023] 1.1. Construct a pre-trained dataset; N (e.g., N=1000) protein structure identifiers (PDB IDs) are selected from large datasets (e.g., PDBbind, CASF, etc.) through random sampling to establish the core protein sample library; Subsequently, based on these N proteins, a large-scale data augmentation operation was carried out, matching each selected protein with M (e.g., M=1000) different ligand molecules (randomly selected from the ChEMBL database), and conducting systematic molecular docking experiments, finally constructing a dataset containing M×N (e.g., M×N=1 million) sets of protein-ligand pairing combinations. Finally, using AutoDock Vina's molecular docking software (for example, setting a cube search box covering the entire pocket with the geometric center of the docking pocket as the origin, and expanding each dimension by approximately...), The margin, docking accuracy, and number of energy assessments are controlled by the exhaustiveness parameter of AutoDock Vina (default value is 8). Quantitative analysis of energy and docking effect is performed on each protein-ligand combination, accurately calculating key energy parameters such as van der Waals attraction, van der Waals repulsion, solvation energy, and hydrogen bond energy, while outputting the corresponding docking score. A pre-training dataset is constructed based on the protein-ligand pairing dataset, energy parameters, and docking score. The energy parameters and docking score serve as label data in the subsequent prediction model pre-training process.

[0024] PDBbind is an authoritative database of biomolecular complexes created by the Shanghai Institute of Organic Chemistry, Chinese Academy of Sciences. ChEMBL is a manually curated database of bioactive molecules with drug-grade properties. It integrates chemical, biological activity, and genomic data to help translate genomic information into effective new drugs.

[0025] 1.2. Constructing the fine-tuning dataset; This embodiment employs quantum Fourier transform (QFT) for efficient homology retrieval on the PDBbind general set, thereby constructing a fine-tuning dataset. The core of this embodiment lies in leveraging the quantum parallelism and efficient frequency domain transformation capabilities of QFT to address the computationally intensive and inefficient nature of traditional homology retrieval on the PDBbind general set. The overall approach involves preprocessing the PDBbind general set and target sequences, then using QFT for feature extraction and matching, ultimately achieving rapid screening of homology sequences. A similarity threshold is set to construct a fine-tuning dataset of Q (e.g., Q = 2000-5000) entries.

[0026] Before performing quantum computing, biological sequence information (referring to the amino acid letter sequence extracted from the PDBbind universal set or the target protein) needs to be transformed into a mathematical expression that can be processed by quantum algorithms. Specifically, the target protein sequence to be retrieved is encoded, and its features are converted into numerical vectors. The known protein sequences in the PDBbind universal set undergo the same preprocessing (i.e., encoding) to ensure that the data format is consistent with the target sequence. Through quantum Fourier transform, the structural or compositional features of the protein sequence are transformed from the time domain (spatial domain) to the frequency domain. Utilizing the superposition state of qubits, the model can process multiple sequence features simultaneously, significantly accelerating feature extraction.

[0027] In this embodiment, the frequency domain fingerprint of the protein is constructed using quantum Fourier transform (QFT), and the calculation process is as follows: The first step is amino acid encoding: each amino acid is represented by an 8-dimensional physicochemical property vector (hydrophobicity, molecular weight, isoelectric point, polarity, hydrogen bond donor ability, polarity, side chain flexibility, and aromaticity).

[0028] The second step is length normalization: the protein sequence signal is resampled using linear interpolation. = 256 dimensions ( (To be determined dynamically based on actual circumstances).

[0029] The third step is frequency domain transformation: Quantum Fourier transform is performed on each of the eight channels to preserve complete complex information (amplitude + phase) and obtain a quantum frequency domain fingerprint in complex form.

[0030] Step 4, similarity calculation: The QFT fingerprint similarity between two proteins is defined as the square of the modulus of the complex inner product, normalized to [0, 1], and used as the similarity score for homology retrieval.

[0031] Specifically: (a1) The first and second steps above are described in further detail below.

[0032] Universal coding of global physicochemical characteristics of proteins; The entire amino acid chain of the target protein sequence is considered as a continuous biological signal stream. Utilizing eight core physicochemical properties (as shown in Table 1) – hydrophobicity, molecular weight, isoelectric point, number of hydrogen bond donors, number of hydrogen bond acceptors, polarity, side chain flexibility, and aromaticity – each residue of the target protein sequence is encoded into an 8-dimensional feature vector.

[0033] Table 1: Eight core physicochemical properties

[0034] Since the number of residues L varies among different proteins, each protein first independently determines the number of qubits based on its own length. (Pick ,Right now (The least power of 2, not less than L), and then each residue is encoded into an 8-dimensional physicochemical property vector to form an L×8 matrix. To achieve uniform processing of proteins of different lengths, an adaptive resampling technique is used (e.g., one-dimensional linear interpolation is performed on each physicochemical property channel along the sequence length direction, and the encoded 8-dimensional physicochemical property vector is uniformly resampled to...). (at each location point), after resampling The physicochemical signal matrix carries the global physicochemical information of the target protein sequence, thereby preserving the whole-chain evolutionary fingerprint while achieving universal coverage of various types of proteins (such as kinases, transmembrane proteins, etc.). Ultimately, each protein is represented by its own defined... Each qubit independently enters the quantum circuit to perform a quantum Fourier transform, achieving adaptive, dimension-fixed, and unified cross-length processing. Among these... This represents the number of qubits used in the quantum Fourier transform.

[0035] Among them, global physicochemical information refers to the multi-channel numerical signal stream with fixed dimensions that is generated by encoding the full-length amino acid sequence of the target protein or PDBbind universal centralized sample through the above-mentioned 8-dimensional physicochemical properties and adaptive resampling. This signal stream carries the distribution law of hydrophobicity, charge, flexibility and other properties on the entire protein chain in the form of numerical sequence, and serves as the input data for amplitude encoding operation.

[0036] It should be noted that the construction of the fine-tuning dataset also includes: when the target protein has a clear domain division, a domain-specific database is used to construct the fine-tuning dataset.

[0037] (a2) The third step above will be described in further detail below.

[0038] Frequency domain conversion based on dynamic qubits and quantum Fourier transform; Based on the dimensions of the resampled physicochemical signal matrix Dynamic allocation Each quantum bit resource (e.g., 12 quantum bits can carry 4096-dimensional features) is used, and amplitude encoding is employed to load the global physicochemical information carried in the physicochemical signal matrix into a quantum superposition state. Subsequently, a standard quantum Fourier transform is performed... (that is, only about) It can be processed by one quantum bit. Frequency domain transformation can be performed with the time complexity of 1D data. Due to the position-independent nature of the quantum Fourier transform, it can automatically capture frequency patterns hidden in long chains. For example, protein sequences have a natural periodic structure: α-helices have 3.6 residues per turn, and β-sheets have approximately 2 residues per cycle. These structural constraints are reflected in the periodic fluctuations of the physicochemical properties of amino acids (hydrophobicity, charge, etc.). Homologous proteins (such as hemoglobin from different species), despite significant sequence differences, share the same tertiary folding structure, resulting in consistent spatial frequency distributions of their physicochemical properties, thus enabling frequency domain fingerprint matching. That is, even if the relative positions of homologous fragments in the sequence change, the quantum frequency domain fingerprint of the target protein sequence can be accurately extracted, thereby solving the problem of traditional algorithms struggling to match distant homologous sequences.

[0039] The mechanism of the quantum Fourier transform is to transform a sequence signal into the frequency domain. In the actual computation of this embodiment, a quantum simulator is used. Its operating environment and resource requirements are determined based on the scale, for example, using more than 64GB of memory and a CPU: Intel(R) Xeon(R) Platinum 8336c CPU @ 2.30GHz.

[0040] Table 2: Physicochemical properties of major amino acids

[0041] This embodiment employs a standard quantum Fourier transform circuit. This circuit takes a physicochemical signal matrix encoded with amplitude and loaded into a quantum superposition state as input, and performs a mathematical mapping through quantum Fourier transform. This digital mapping transforms the sequence signal carrying global physicochemical information from the spatial domain (sequence position) to the frequency domain (frequency components). The complex information output after the transformation directly serves as the quantum frequency domain fingerprint of the protein, completely preserving amplitude and phase information. This standard quantum Fourier transform circuit is a well-known architecture in the field of quantum computing; therefore, this embodiment will not elaborate on the circuit details, but focuses on utilizing its frequency domain transformation capability to extract fingerprint features for subsequent homology retrieval.

[0042] III. Autonomous homology retrieval and construction of a Q-level fine-tuning dataset; Efficient autonomous retrieval is achieved by comparing the quantum fingerprint fidelity of target protein sequences with samples in the PDBbind universal dataset. Specifically, the system automatically sorts the entire database based on frequency domain similarity scores, identifying samples with highly correlated global physicochemical properties without manual intervention. By setting a rigorous similarity threshold, Q high-quality data points (e.g., 2000 to 5000) are accurately extracted from the universal dataset. This self-constructed fine-tuning dataset not only possesses strong target specificity but also significantly enhances the performance of transfer learning models in the domain characteristic calibration stage due to its coverage of global homology information.

[0043] The fourth step above is described in further detail below: the quantum fingerprint fidelity is used as the similarity score, and the similarity calculation steps are as shown in (b1) to (b3).

[0044] (b1) The quantum superposition state carrying global physicochemical information is transformed from the spatial domain to the frequency domain by quantum Fourier transform, and a complex quantum frequency domain fingerprint vector (i.e. the complex vector output by quantum Fourier transform itself) is obtained.

[0045] The target protein and the PDBbind universal centralized samples were processed through the above process to obtain their respective fingerprint vectors (target sequence fingerprint and sample sequence fingerprint).

[0046] (b2) Calculate the vector dot product; The similarity of the frequency domain fingerprints of two proteins is measured using the quantum state inner product (i.e., multiplying the conjugate of the previous complex vector with the subsequent vector element-wise and summing the results). The closer the squared modulus of this inner product is to 1, the more consistent the frequency distribution of the global physicochemical properties of the two proteins are.

[0047] (b3) Calculate the square of the modulus (similarity score); The similarity score is the square of the absolute value of the inner product of the quantum states.

[0048] Traditional sequence alignment (such as BLAST) is based on direct matching of letters (A, W, G, etc.). However, in biology, there is the phenomenon of distant homology: two proteins may have very different letter sequences, but due to the consistent frequency distribution of their physicochemical properties, they share similar 3D structures and binding pockets. In this embodiment, the quantum fingerprint fidelity (similarity score) calculated by quantum Fourier transform reflects the resonance of physicochemical laws. Even if the amino acids are different, as long as the hydrophobic-hydrophilic rhythm is consistent, the quantum fingerprint fidelity will be high. If the quantum fingerprint fidelity is greater than the similarity threshold, it can retrieve key data that was missed by traditional algorithms from the PDBbind general set.

[0049] This embodiment employs a two-stage transfer learning strategy: large-scale pre-training followed by specific fine-tuning. The pre-training stage utilizes 1 million (optional) sets of large-scale molecular docking data (1000 (optional) proteins × 1000 (optional) ligands) to enable the prediction model to learn the general physical laws governing intermolecular interactions. The fine-tuning stage requires only 2000 to 5000 (optional) high-quality homology samples to achieve high-precision predictions for specific target domains. This contrasts sharply with traditional machine learning models that directly rely on scarce experimental affinity data (10... 3 Up to 10 4 This embodiment significantly reduces the reliance on high-quality experimental data while maintaining extremely high prediction accuracy in small-sample, specific target domains.

[0050] This embodiment applies Quantum Fourier Transform (QFT) to protein homology retrieval. Unlike the traditional BLAST algorithm, which relies on direct matching of amino acid letters, this embodiment transforms protein sequences into a continuous stream of physicochemical signals (hydrophobic waveforms, charge waveforms, etc.) through 8-dimensional physicochemical property encoding. After quantum Fourier transform to the frequency domain, quantum frequency domain fingerprints are extracted. Even if two protein sequences have significantly different amino acid letters, as long as their physicochemical property distribution patterns (such as the periodicity of the hydrophobic core and the charge distribution pattern) are consistent, they can be identified as homologous samples through high-frequency domain similarity scores. This mechanism can effectively retrieve distant homologous samples missed by traditional algorithms from the PDBbind general set, significantly broadening the coverage of the fine-tuning dataset.

[0051] Before inputting the pre-trained dataset or fine-tuned dataset into the prediction model, this embodiment needs to construct a fusion graph based on the samples. The specific construction process is as follows (c1) to (c3): (c1) Construct a protein pocket diagram with Cα atoms of protein binding pocket residues as residue nodes and the spatial distance between residues less than a preset threshold as edges; (c2) Construct a small molecule graph using ligand atoms as nodes and chemical bonds between ligand atoms as edges; Example: (c2-1) Determination of the characteristics of ligand atom nodes.

[0052] To enable the structure encoder to fully perceive the chemical nature of the ligand molecule, this embodiment constructs a 15-dimensional initial feature vector for each ligand atom. This feature vector covers four levels: atom type, electronic structure, topological environment, and physicochemical properties, as specifically defined in Table 3.

[0053] Table 3: Ligand Atom Node Feature Definitions (15 Dimensions)

[0054] The selection of the 15 features in Table 3 is based on the following: Atom type, atomic mass, and Gasteiger charge determine the fundamental strength of intermolecular van der Waals forces and electrostatic interactions (corresponding to van der Waals attraction / repulsion in the pre-training labels); hybridization type, bond valence, and implicit / explicit valence reflect the saturation and spatial orientation of atomic bonding, directly affecting the three-dimensional conformational stability of ligand molecules; in-ring atoms and aromaticity markers are used to identify potential donors for π-π stacking interactions (corresponding to aromatic interactions in Table 1); the total number of hydrogen atoms and hydrogen bond donor / acceptor capacity (implicit in Gasteiger charge and atom type) directly determine the accuracy of hydrogen bond energy calculation. By normalizing all 15 features to the [0,1] interval, it can be ensured that features of different dimensions maintain gradient stability during backpropagation of the neural network, avoiding slow convergence caused by differences in feature scale.

[0055] (c2-2) Determination of the chemical bond edge characteristics of ligands.

[0056] After constructing the ligand atom nodes, this embodiment further defines 4D edge features for the chemical bonds between atoms to characterize the physicochemical properties of the bonds, as shown in Table 4.

[0057] Table 4: Definition of Ligand Chemical Bond Edge Characteristics (4-dimensional)

[0058] In Table 4, bond-level encoding is obtained by using the `GetBondType()` method of RDKit to retrieve the bond type (SINGLE, DOUBLE, TRIPLE, AROMATIC), mapped sequentially, and normalized by dividing by 3.0. Conjugation is determined using the `GetIsConjugated()` method, which identifies a multiple bond connected by a single bond to another multiple bond as conjugated. In-ring encoding is determined by the `IsInRing()` method to determine if a chemical bond belongs to a ring system. Bond length encoding is obtained by using the `MolTransforms.getBondLength()` method to retrieve the Euclidean distance (in units) between two bonding atoms in the 3D conformation. Then divide by 5.0 for normalization. All of the above RDKit interfaces are well-known tool functions in the field of cheminformatics. Those skilled in the art can fully reproduce the extraction process of ligand chemical bond edge features based on the above description.

[0059] The selection of the 4D edge features in Table 4 takes into account the following factors: bond order encoding distinguishes single, double, triple, and aromatic bonds. Different bond orders correspond to different bond energies and rotational barriers, directly determining the distribution of flexible regions of ligand molecules, and thus affecting their induced fit to the protein binding pocket (corresponding to side chain flexibility in Table 1); conjugation markers are used to identify delocalization. Electronic systems, which often possess rigid planar structures, can significantly enhance the interaction with aromatic residues (such as phenylalanine and tyrosine). Stacking effects; intra-ring markers are used to distinguish between ring-based and open-chain structures. Ring-based structures have lower conformational entropy and lose less free energy during bonding, which is beneficial for enhancing affinity; bond length encoding provides spatial geometric constraints through measured distances in the 3D conformation, assisting the graph attention network in correctly perceiving the spatial proximity relationships between atoms. All the above edge features are also normalized and together with the node features, constitute the complete input representation of the ligand small molecule graph.

[0060] After assigning values ​​to the ligand atom node features (Table 3) and chemical bond edge features (Table 4), the complete ligand small molecule diagram data structure is obtained, and then the merging operation (c3) is performed with the aforementioned protein pocket diagram.

[0061] (c3) The protein pocket diagram is merged with the small molecule diagram, and cross-molecular connection edges are established between protein residue nodes and ligand atom nodes according to spatial proximity to obtain the fusion diagram.

[0062] Specifically, the Euclidean distance between the Cα atom of the protein pocket residue and the ligand atom is calculated. If it is less than a preset threshold (e.g., ...), the result is considered. Then, a transmolecular connection edge is established.

[0063] The pre-training dataset or fine-tuning dataset, which has been transformed into a fusion graph, is partitioned. Specifically, the pre-training dataset is divided into a pre-training set and a pre-test set according to a certain ratio (e.g., 7:3 or 8:2, etc.), and the fine-tuning dataset is divided into a fine-tuning set and a fine-tuning test set. The prediction model is trained using the pre-training set and tested using the pre-test set; the trained prediction model is fine-tuned using the fine-tuning set and tested using the fine-tuning test set.

[0064] This embodiment merges the protein pocket map and the small molecule map into a unified fusion map. Compared to traditional methods that process proteins and ligands separately or employ mid-to-late stage fusion strategies, the fusion map enables the structure encoder to model intermolecular interactions simultaneously at the atom-to-atom, residue-to-residue, and atom-to-residue levels. Combined with deep feature extraction at the graph attention layer, this embodiment achieves a more fine-grained and accurate representation of interaction features than the physical force field approximation of traditional docking software.

[0065] II. Model Architecture and Training; 2.1. The prediction model consists of a structural encoder and an affinity predictor connected in sequence; The structure encoder takes the fused graph as input and performs deep processing through the graph attention layer. By capturing complex interaction information at the atomic and residue levels layer by layer, it finally generates a feature vector containing comprehensive structural and interaction features, providing feature support for subsequent prediction tasks.

[0066] The affinity predictor, consisting of fully connected layers, has the core function of receiving the feature vector output by the structural encoder. Through the stepwise mapping and feature transformation of the fully connected layers, it predicts the various energy terms (such as van der Waals attraction / repulsion, hydrogen bond energy, solvation energy, etc.) and the total docking score calculated by the molecular docking software of AutoDock Vina during the pre-training stage, based on the prediction results corresponding to different outputs in the training stage.

[0067] For example, the structural encoder uses a graph attention network (GAT) with 4 layers. Its input is the fused graph, and its output is a 512-dimensional feature vector. The affinity predictor uses a fully connected network (MLP) with 4 layers. Its input is the 512-dimensional feature vector output by the structural encoder. During the pre-training phase, its output is a 5-dimensional feature vector of [van der Waals attraction, van der Waals repulsion, solvation energy, hydrogen bond energy, docking fraction]. During the fine-tuning phase, its output is [dissociation constant]. Suppression constant Half-maximal inhibitory concentration The affinity predictor consists of four layers: FC1, FC2, FC3, and FC4. FC1 uses Linear+ReLU+Dropout with a 512-dimensional input and output; FC2 uses Linear+ReLU+Dropout with a 512-dimensional input and 256-dimensional output; FC3 uses Linear+ReLU+Dropout with a 256-dimensional input and 128-dimensional output; and FC4 uses Linear with a 128-dimensional input and either a 5-dimensional or 3-dimensional output. Here, Linear is a linear layer, ReLU is an activation layer, and Dropout is a regularization parameter that can be set to 0.2.

[0068] During the fine-tuning phase, the binding affinity indices of protein-small molecule interactions are accurately predicted. These indices include the dissociation constant, inhibition constant, and half-maximal inhibitory concentration (IC50).

[0069] In molecular docking and ligand-receptor interaction studies, affinity indices are used to measure the strength of molecular binding, and dissociation constants are also considered. This describes the equilibrium state of reversible binding between ligands and receptors; the smaller the value, the stronger the affinity. Furthermore, when half of the receptor binding sites are occupied... Equal to free ligand concentration; inhibition constant ( The half-maximal inhibitory concentration (MCC) reflects the affinity between the inhibitor and the target; the smaller the value, the stronger the inhibitory effect, which can be deduced through enzyme kinetics. The concentration at which an inhibitor achieves 50% inhibition is determined by experimental conditions such as substrate concentration and reaction system. It directly reflects the inhibitory activity of the inhibitor, and together the three factors provide the core basis for the quantitative analysis of intermolecular binding interactions.

[0070] This embodiment employs multi-task learning during the fine-tuning phase, enabling simultaneous output of the dissociation constant ( ), suppression constant ( ) and half-maximal inhibitory concentration (IC50) Three core binding affinity indicators. Compared to traditional methods that only predict a single docking fraction, this embodiment provides a more comprehensive quantitative analysis of intermolecular binding interactions.

[0071] 2.2 Training Strategy; Pre-training phase: Learning general structural features; The training objective of this stage is to predict the docking score and various energy terms of AutoDock Vin. During the training process, the prediction model architecture is fixed, and a large-scale docking PDBbind general set of 1 million pairs (configurable) is used to perform end-to-end training on the structure encoder and affinity predictor. This allows the prediction model to learn the general structural interaction rules between proteins and small molecule atoms. The key result is that the prediction model can reproduce the various energy terms output by traditional docking software with high accuracy.

[0072] During the pre-training phase, the loss function used is MSE: ,in, This is the batch sample size, for example, 50. For the first The model prediction vector for each sample (pre-trained in 5 dimensions / fine-tuned in 3 dimensions). For the first The true label vector of each sample It is a 2-norm. This represents the batch size.

[0073] Fine-tuning phase: Optimizing specific affinity prediction; This training phase is set to fix the pre-trained structural encoder parameters to retain the general structural knowledge of proteins and small molecules it has learned, and only the affinity predictor parameters are fine-tuned. The training objective is to accurately predict the three core binding affinity indices of protein-small molecule interactions (dissociation constant, inhibition constant, and half-inhibition concentration). The training method is multi-task learning to achieve simultaneous prediction of multiple indices.

[0074] III. Application of Predictive Models; After pre-training and fine-tuning, the prediction model can be directly used to predict protein-small molecule binding affinity.

[0075] (d1) Input data preprocessing: For the target protein to be predicted, the binding pocket region is first extracted using molecular simulation tools. The Cα atom of the protein binding pocket residue is used as the residue node, and the spatial distance between residues less than a preset threshold is used as the edge to construct the protein pocket graph. At the same time, the small molecule Mol2 files to be screened are processed using the RDKit tool (an open-source cheminformatics toolkit) to remove molecules containing invalid bonds or unreasonable configurations. Then, the small molecule graph is generated using small molecule atoms as nodes and the chemical bonds between ligand atoms as edges. The two types of graph data are then merged into a fusion graph, which is used as the input for the prediction model.

[0076] (d2) Predictive model inference: The preprocessed fusion graph is input into the trained predictive model. The structure encoder will automatically extract the key feature vectors of proteins and small molecules based on the general structural interaction rules learned in the pre-training stage, without the need to retrain the structure encoder. Then, the affinity predictor will output the three core affinity indicators of the protein-small molecule combination and the energy terms (van der Waals attraction / repulsion, hydrogen bond energy, solvation energy, etc.) based on the fine-tuned parameters.

[0077] (d3) Results Analysis: Candidate small molecules were screened based on the predicted affinity index values— , The smaller the value, The lower the value, the stronger the binding ability of the small molecule to the target protein and the better the inhibitory activity. Simultaneously, combining energy data allows for further analysis of the dominant factors in intermolecular interactions (such as hydrogen bonding or van der Waals forces), providing direction for molecular modification and optimization. This predictive model can be widely applied in scenarios such as virtual screening of drug molecules and optimization of lead compounds. Compared to traditional molecular docking software, it can improve prediction efficiency and the ability to simultaneously evaluate multiple indicators.

[0078] IV. Example 1 To verify the effectiveness of the large-scale pre-training + domain fine-tuning two-stage transfer learning strategy proposed in this embodiment, the following comparative experiment was designed: using the same prediction model architecture (structural encoder GAT + affinity predictor MLP), the final prediction accuracy of the two training strategies was compared on the same kinase target fine-tuning dataset.

[0079] (e1) Experimental setup Dataset: The kinase family was selected as the target domain, and a fine-tuned dataset containing 2,000 samples was constructed (covering major kinase subtypes such as CDK, EGFR, and VEGFR).

[0080] Comparison of options: Option A (transfer learning, i.e. this embodiment): First, pre-train on 1 million sets of general molecular docking data (PDBbind general set) for 150 rounds, then fine-tune on kinase fine-tuning dataset for 100 rounds after fixing the structure encoder parameters.

[0081] Option B (direct training, i.e. no transfer learning): without pre-training, train directly from scratch on the same kinase fine-tuning dataset for 250 epochs (the total number of epochs is the same as the total number of pre-training + fine-tuning epochs in Option A, to eliminate the impact of the difference in the number of training steps).

[0082] Evaluation metrics: Root mean square error (RMSE) is used as a quantitative indicator of prediction accuracy. Three affinity indices (dissociation constant, ...) are calculated separately. ), suppression constant ( ), half-maximal inhibitory concentration ( The predicted RMSE value on the independent test set.

[0083] (e2) Experimental Results and Analysis The verification results are as follows Figure 2 As shown in Table 5.

[0084] Table 5: Comparison of predicted RMSE under different training strategies (kinase dataset)

[0085] As shown in Table 5, scheme A, which employs a transfer learning strategy, significantly outperforms scheme B, which is trained directly from scratch, across all three affinity metrics. The average RMSE decreased from 1.95 to 1.74, a reduction of 11.2%. Further analysis is needed to understand the reasons for this advantage.

[0086] First, the effective transfer of prior knowledge of physical laws. After pre-training on 1 million sets of large-scale general docking data, the structural encoder has fully learned the general non-bonded interaction rules between protein and ligand atoms (such as van der Waals potential surfaces, hydrogen bond directionality, hydrophobic collapse driving, etc.). When transferring to kinase targets, the encoder does not need to learn the basic physical rules between atoms from scratch, but directly uses the general knowledge as a starting point, only fine-tuning to adapt to the residue distribution and charge characteristics specific to the kinase pocket. In contrast, the direct training scheme is limited by 2,000 small sample data, making it difficult to fit the complex physical mapping in the high-dimensional feature space from scratch, leading to overfitting and insufficient generalization ability on the test set.

[0087] Second, overfitting suppression in small-sample scenarios. The structural encoder in this embodiment contains approximately 3.8 million trainable parameters (4-layer GAT). If trained directly from scratch on a small-scale kinase dataset, the number of parameters (3.8M) is far greater than the sample size (2K), making overfitting highly likely. The transfer learning strategy, by freezing the encoder parameters and fine-tuning only the affinity predictor (approximately 0.43M), compresses the effective number of trainable parameters to about one-tenth of the original, significantly reducing the model's dependence on sample size. This mechanism enables high-accuracy prediction results even in specific target domains with scarce labeled data (such as insufficiently studied kinase mutants), validating the practical value of this approach in data-scarce scenarios.

[0088] Third, improved convergence speed and training efficiency. Experimental data shows that scheme A only requires about 30 epochs to reach convergence during the fine-tuning phase (total epochs 150 + 30 = 180), while scheme B requires about 180 epochs to approach convergence from scratch, and the final convergence value is much higher than that of scheme A. Transfer learning not only improves the final accuracy but also significantly accelerates the convergence process and reduces the computational resource consumption of model training.

[0089] (e3) Overall Conclusion Based on the comparative experiments above, the transfer learning strategy proposed in this embodiment reduces the average RMSE by 11.2% on the kinase dataset, demonstrating the significant superiority of the two-stage training paradigm of "large-scale pre-training + domain fine-tuning" in molecular docking affinity prediction tasks. This strategy successfully solves the technical bottleneck of traditional deep learning models' inaccurate predictions due to insufficient data in small-sample, specific target domains, and is the core technical support for this invention.

[0090] V. Verification of QFT fingerprint homology retrieval performance To verify the ability of the quantum Fourier transform (QFT) fingerprint in this embodiment to retrieve distant homologous proteins, representative proteins covering multiple SCOP fold types were selected as the verification dataset, as shown in Table 6.

[0091] Table 6: Composition of the validation dataset

[0092] Table 6 shows that UniProt (UniProt Knowledgebase, or UniProtKB for short) is currently the most comprehensive and information-rich database of protein sequences and functions internationally. Globin is a globin, Serine Protease is a serine protease, Cytochrome c is cytochrome c, Lysozyme is lysozyme, RNase is a ribonuclease, and Ubiquitin is ubiquitin.

[0093] The homology criteria for the validation datasets in Table 6 are as follows: protein homology was determined by SCOP 2.08 (Structural Classification of Proteins) through DALI structure alignment. Long-range homology equals the same SCOP fold with sequence identity <30%; non-homology equals different SCOP folds. Based on these criteria, a total of 66 pairs of analytical objects were constructed, including 10 pairs of long-range homology pairs and 56 pairs of non-homology pairs.

[0094] (f1) Validate the analysis of the criteria for selecting the dataset; Table 6 shows the proteins selected, covering the three most representative fold spaces in the SCOP database: class a (fully α-helical, such as the globin family and cytochrome c), class b (fully β-sheet, such as the trypsin-like serine proteases family), and class d (α+β mixed folds, such as lysozymes and nucleases). Each fold contains members from different species with varying sequence lengths and sequence similarities below 30% (e.g., human hemoglobin α-chain and soybean leghemoglobin, with only about 24% sequence similarity, but sharing the same globin fold), thus providing a challenging validation scenario for long-distance homology retrieval. Simultaneously, different folds (e.g., globin and trypsin) constitute a clearly defined set of non-homologous negative samples. This design ensures that the validation dataset comprehensively examines the retrieval robustness of the QFT fingerprinting method under different fold types, sequence lengths, and evolutionary distances.

[0095] (f2) Detection performance comparison and analysis; The QFT fingerprinting method of this embodiment was compared with the classic sequence alignment algorithm BLAST (which uses Smith-Waterman local alignment and BLOSUM62 scoring matrix, with the threshold set to E-value < 0.01 and including database size correction) on the same validation set. The results are shown in Table 7.

[0096] Table 7: Comparison of Remote Homologous Detection Performance

[0097] As shown in Table 7, the BLAST method detected 9 pairs of long-range homologs on this validation set (missing 1 pair); while the QFT fingerprint method in this embodiment detected 8 pairs (missing 2 pairs). Although the detection rate of QFT fingerprint is slightly lower than that of BLAST, it should be noted that: firstly, both methods maintained a zero false positive rate on 56 pairs of non-homologous samples, indicating that QFT fingerprint does not misclassify non-homologous proteins as homologous, and its retrieval reliability fully meets the requirements for constructing the fine-tuning dataset; secondly, the two pairs of long-range homologs missed by QFT fingerprint are proteins with large differences in sequence length (such as Human Ubiquitin and some d-type folded proteins) or with relatively special domain compositions, belonging to marginal cases with blurred boundaries in the current SCOP classification, and traditional algorithms also miss such cases. Overall, QFT fingerprint achieves a long-range homolog identification capability comparable to or even close to traditional sequence alignment methods with a detection rate of 8 pairs (missing 2 pairs), fully verifying its effectiveness as a homologous retrieval tool.

[0098] (f3) Comparative Analysis of Algorithm Characteristics and Efficiency The core computational characteristics of the two methods are further compared, as shown in Table 8.

[0099] Table 8: Comparison of characteristics between QFT fingerprinting and BLAST methods

[0100] Further analysis based on Table 8 reveals that the core bottleneck of the traditional BLAST algorithm lies in its dynamic programming alignment complexity of O(mn). When processing large-scale databases (such as the PDBbind general set containing tens of thousands of sequences), the computational cost of a single retrieval increases quadratically with the sequence length. In contrast, this embodiment encodes the protein sequence into a fixed-dimensional frequency domain fingerprint (fixed to 256 dimensions after resampling; in actual implementation, this can be extended to 256 dimensions based on available qubit resources). n The QFT fingerprint method simplifies the similarity calculation between any two sequences to a complex inner product between pre-computed fingerprints, reducing the complexity of a single comparison to O(1) and the complexity of a full database search to only O(M), completely decoupling it from sequence length. In large-scale virtual screening scenarios, this efficiency advantage will amplify exponentially with the increase in database size. More importantly, the QFT fingerprint method has a natural quantized transfer path: on a quantum computer, the quantum Fourier transform requires only O(log²N) quantum gate operations, and combined with quantum superposition states, it can realize parallel encoding and measurement of fingerprints across the entire database, which is fundamentally impossible for the serial architecture of the classical BLAST algorithm.

[0101] (f4) Overall Conclusion Based on the above verification results, although the QFT fingerprinting method in this embodiment has a slightly lower detection rate than BLAST, its advantages are reflected in three aspects: First, the zero false positive rate ensures the purity of the retrieval results and avoids contamination of the fine-tuning dataset by non-homologous samples; second, constant-time alignment significantly improves the retrieval efficiency of large-scale databases, making it suitable for screening pre-trained general sets containing massive amounts of samples; third, the algorithm architecture is naturally compatible with future quantum computing hardware and has good potential for technological evolution. Therefore, using QFT fingerprinting for autonomous homologous retrieval to construct a fine-tuning dataset can significantly improve efficiency while ensuring retrieval quality, and is an indispensable key technical step in this embodiment.

[0102] In this embodiment, BLAST is based on positional alignment of sequence letters, requiring amino acid residues to be consistent in position. When distant homologous proteins undergo prolonged divergent evolution, the amino acid sequence may become unrecognizable, rendering BLAST unable to identify them. In contrast, QFT encodes the entire protein sequence into a continuous biological signal stream (hydrophobic waveforms, charge waveforms, etc.) through 8-dimensional physicochemical properties. After Fourier transform to the frequency domain, it extracts the periodic distribution pattern of global physicochemical properties—such as the frequency of hydrophobic cores and the spacing of charge-rich regions. Even if two proteins have completely different amino acid sequences, as long as their physicochemical property distribution rhythms are consistent (the frequency of alternating hydrophobic and hydrophilic regions is the same), their frequency domain fingerprints will be highly similar. This physicochemical resonance recognition mechanism helps bypass letter-level differences and capture the physicochemical essence that determines the three-dimensional structure and binding pocket functional properties, thus retrieving distant homologous samples missed by BLAST. This improves the prediction accuracy, generalization ability, and multi-index comprehensive evaluation level of the prediction model in this application for specific target areas.

[0103] Based on the above description of the embodiments, those skilled in the art will understand that the molecular docking affinity prediction method, system, and medium based on transfer learning described in this embodiment can be implemented in pure software or deployed and run on general-purpose or dedicated computing hardware platforms. Based on this essence, the technical solution of this embodiment can be specifically implemented in the form of a software product containing program instructions. This software product can be stored on various non-volatile storage media or directly deployed as a local or cloud service. The program instructions are used to cause computer devices with processing capabilities—including but not limited to personal computers, server clusters, mobile terminals, or other network devices—to execute the steps described in this embodiment.

[0104] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A molecular docking affinity prediction method based on transfer learning, characterized in that, Input the protein-ligand combination to be predicted into the prediction model, and output the binding affinity index. The training process of the prediction model is as follows: Protein samples are selected from a pre-defined general set, multiple ligands are matched for each protein sample and molecular docking is performed to construct a pre-training dataset, and energy term parameters and docking scores are calculated as pre-training labels. The target protein sequence is sequentially encoded and subjected to quantum Fourier transform. The extracted quantum frequency domain fingerprint is compared with the samples in a preset general set to calculate the fingerprint fidelity. Homologous samples are then selected to construct a fine-tuning dataset. Construct a predictive model, which includes a structural encoder and an affinity predictor connected in sequence; Pre-training phase: The prediction model is trained end-to-end using the pre-training dataset and pre-training labels; Fine-tuning phase: Using the fine-tuning dataset as input, with the structure encoder parameters fixed, the affinity predictor parameters are fine-tuned with the goal of predicting the binding affinity index.

2. The method according to claim 1, characterized in that, The construction of the pre-trained dataset includes: Protein structure identifiers were screened from the general PDBbind database to create a protein sample library; For each protein, multiple ligands were matched from the ChEMBL database for molecular docking. The energy parameters, including van der Waals attraction, van der Waals repulsion, solvation energy, and hydrogen bond energy, as well as the docking score, were calculated using molecular docking software and used as pre-training labels.

3. The method according to claim 1, characterized in that, The specific steps for encoding the target protein sequence are as follows: The entire amino acid chain of the target protein sequence is considered as a continuous biological signal flow; Each residue of the target protein sequence is encoded using 8-dimensional physicochemical properties, and an adaptive resampling technique is used to uniformly map the encoded feature vector along the sequence length dimension. At each location point, a resampled physicochemical signal matrix is ​​obtained, which carries the global physicochemical information of the target protein sequence; The eight physicochemical properties include hydrophobicity, molecular weight, isoelectric point, number of hydrogen bond donors, number of hydrogen bond acceptors, polarity, side chain flexibility, and aromaticity. This represents the number of qubits used in the quantum Fourier transform.

4. The method according to claim 3, characterized in that, The quantum Fourier transform includes: Based on the dimensions of the physicochemical signal matrix Dynamic allocation One quantum bit resource; The global physicochemical information carried in the physicochemical signal matrix is ​​loaded into a quantum superposition state using amplitude encoding; The quantum superposition state is transformed from the spatial domain to the frequency domain by quantum Fourier transform, and the quantum frequency domain fingerprint of the target protein sequence is extracted.

5. The method according to claim 1, characterized in that, Before inputting the pre-trained or fine-tuned dataset into the prediction model, a fusion graph needs to be constructed based on the samples. The specific construction process is as follows: A protein pocket diagram is constructed using the Cα atoms of protein binding pocket residues as residue nodes and the spatial distance between residues less than a preset threshold as edges. Small molecule graphs are constructed using ligand atoms as nodes and the chemical bonds between ligand atoms as edges. The protein pocket diagram is merged with the small molecule diagram, and cross-molecular connection edges are established between protein residue nodes and ligand atom nodes based on spatial proximity to obtain the fusion diagram.

6. The method according to claim 1, characterized in that, The binding affinity indices in the fine-tuning phase include one or more of the dissociation constant, inhibition constant, and half-maximal inhibitory concentration. The fine-tuning employs multi-task learning to achieve simultaneous prediction of multiple indices.

7. The method according to claim 1, characterized in that, The construction of the fine-tuning dataset also includes: when the target protein has a clear domain division, the fine-tuning dataset is directly constructed using a domain-specific database.

8. The method according to claim 1, characterized in that, The energy parameters and docking fractions were calculated by AutoDockVina.

9. A computer system comprising a memory, a processor, and a computer program stored in the memory, characterized in that, The processor executes the computer program to implement the method according to any one of claims 1-8.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a plurality of computer programs, which are used to be invoked by a processor and to execute the method as described in any one of claims 1-8.

Images