A molecular generation method based on pharmacophore perception prior and multi-target flow matching

Abstract

The application discloses a molecule generation method based on pharmacophore perception prior and multi-target flow matching, and relates to the technical field of artificial intelligence and computer-aided drug design. The method comprises the following steps: determining the number of ligand atoms and initializing a noise state; constructing a fusion guide velocity field, embedding multiple attribute constraints such as pharmacophore specificity, drug-likeness and synthetic accessibility into an unconditional vector field in a weighted projection manner, and eliminating guide direction contradictions through a conflict perception mechanism; driving a common differential equation solver from noise integration to a terminal state by using the velocity field, performing discrete decoding on the terminal state to construct a molecule graph; screening effective molecules through chemical valence and connectivity verification, and triggering resampling until the standard is met or the upper limit of retry is reached. Multi-target collaborative optimization is realized, the generated molecules have high target binding capacity, good drug-likeness and synthetic feasibility, and the structural effectiveness and spatial rationality are significantly better than those of existing methods, and the method is suitable for target-oriented rational drug design.

Description

Technical Field

[0001] This invention relates to the fields of artificial intelligence and computer-aided drug design technology, and in particular to a molecular generation method based on pharmacophore-aware priors and multi-objective flow matching. Background Technology

[0002] Structure-based drug design (SBDD) aims to generate small ligand molecules that bind with high affinity and selectivity to target proteins based on their three-dimensional spatial structure. Deep generative models, particularly diffusion and flow-matching models, have made significant progress in this field. These models can learn chemical principles from data to generate novel molecules.

[0003] However, existing deep generative methods still face many challenges in practical applications. First, most methods only implicitly learn protein-ligand interaction patterns, lacking explicit guidance on key physicochemical interactions such as hydrogen bonding and hydrophobic interactions, resulting in insufficient chemical complementarity between the generated molecules and the target. Second, drug discovery is a typical multi-objective optimization problem, requiring a simultaneous trade-off between multiple key attributes such as binding affinity, drug-likeness, synthetic accessibility, lipid-water partition coefficient, and topological polar surface area. Existing methods often focus on optimizing a single objective (such as binding affinity) or simply weighting and merging multiple objectives, failing to achieve flexible and decoupled control of attributes and struggling to handle potential conflicts between different attribute optimization directions. Finally, molecular generation requires ensuring the rationality of both its two-dimensional topological structure and three-dimensional geometric conformation, but existing methods often treat these two aspects separately, potentially leading to the generation of invalid or non-physical molecular structures.

[0004] Therefore, there is an urgent need for an efficient molecular generation method that can integrate complementary prior knowledge of targets, synergistically optimize multiple drug properties, and ensure the rationality of molecular structure. Summary of the Invention

[0005] The purpose of this invention is to provide a molecular generation method based on pharmacophore-aware prior and multi-target flow matching, aiming to solve or improve at least one of the above-mentioned technical problems.

[0006] To achieve the above objectives, the present invention provides the following solution: A molecular generation method based on pharmacophore-aware prior and multi-objective flow matching includes: Obtain the atomic coordinates and atomic features of the target protein pocket, and use a pre-trained physicochemical perception prior network to predict the multi-channel pharmacophore distribution prior corresponding to the target protein pocket; The pharmacophore distribution prior of the multi-channel pharmacophore is aggregated and fused with the global features of the protein to obtain the pharmacophore condition vector. The pharmacophore condition vector is then used as the input of the generation condition to couple the SE3-isovariant generation network, so that the atomic coordinates, atomic type and chemical bond type of the ligand molecule are updated together in the same generation process. At each time step of the flow matching sampling, the difference between the conditional vector field and the unconditional vector field of multiple drug-likeness attributes is calculated by a multi-objective classifier-free guidance mechanism. Conflict-aware projection resolution is performed on the guidance direction of attributes with negative correlations, and the resolved attribute guidance directions are weighted and fused to obtain the fused guidance velocity field. The fusion-guided velocity field drives the ordinary differential equation solver to integrate from the noisy state to the final state. The continuous variables of the final state are discretized and decoded to construct a molecular graph. After chemical valence state and connectivity verification and screening, the effective target molecules that meet the multi-objective optimization requirements are output. If the verification fails, a resampling mechanism is triggered until an effective molecule is obtained or the preset retry limit is reached.

[0007] Further training of the physicochemical perception prior network includes: For the protein pocket of the first The atom in the first The pharmacodynamic field intensity on each channel is calculated using a Gaussian kernel function, and the expression is: In the formula, For the real scene; For the first A set of reference ligand atom indices for pharmacophore-like groups; For the corresponding to the first Preset Gaussian kernel width for each channel; For the pocket number A three-dimensional Cartesian coordinate vector of an atom; For reference ligands, the first A three-dimensional Cartesian coordinate vector of an atom; An SE3-isovariable graph neural network was constructed as a physicochemical perception prior network, with the atomic coordinates and atomic features of the protein pocket as input; Train the SE3-isovariable graph neural network to minimize the mean square error between the pharmacodynamic field predicted by the SE3-isovariable graph neural network and the actual field; Through training, the SE3-isovariant graph neural network learns to predict the spatial distribution of multichannel pharmacophores directly from the pocket structure.

[0008] Furthermore, the mean square error between the pharmacophore field predicted by the SE3-isovariant graph neural network and the actual field is expressed as: In the formula, Mean square error; This represents the total number of pharmacophore channels. The number of atoms in the pocket; For the predicted efficacy of the medicine; This is a real-world scenario.

[0009] Furthermore, the multi-channel pharmacophore distribution prior is aggregated and fused with the global protein features to obtain the pharmacophore condition vector. This pharmacophore condition vector is then used as the input to generate the SE3-isovariant generative network, which includes: The number of pocket-ligand complexes of all proteins in the training set is counted, and a conditional probability distribution of pocket atom number versus ligand atom number is constructed and fitted to a truncated normal distribution conditioned on the pocket atom number. Based on the protein pocket volume and the number of atoms in the pocket, the mode is extracted from the truncated normal distribution as the number of atoms in the ligand molecule to be generated. ; Global protein features were obtained by global pooling of pocket atom features using the SE3-equivariant encoder. The multi-channel pharmacophore distribution prior is aggregated and fused with the global features of the protein to obtain a fixed-dimensional pharmacophore condition vector, and the pharmacophore condition vector is used as the input of the generation condition to couple the SE3-equivariant generative network. Coupled SE3-equivariant generative network in time molecular state Protein pocket information, pharmacophore condition vector, and time. As input, output a joint prediction of the final atomic coordinates, atomic type, and chemical bond type.

[0010] Furthermore, the atomic coordinates, atom types, and chemical bond types of the ligand molecule are updated collaboratively during the same generation process, including: Within the SE3-equivariant generative network, the first The layer update process simultaneously handles atomic features, atomic coordinates, and chemical bond features, and injects pharmacophore condition vectors into the updates of atomic and chemical bond features, expressed as: In the formula, and Atoms In the Layer and first Hidden feature vectors of the layer; and Atoms In the Layer and first The three-dimensional Cartesian coordinate vector of the layer; For atoms and The chemical bond between them is in the first The feature vector of the layer; , , These are the learnable multilayer perceptrons corresponding to atomic feature updates, coordinate updates, and chemical bond feature updates, respectively. For atoms The set of neighboring atoms; This is the pharmacophore condition vector.

[0011] Furthermore, the coupled SE3-equivariant generative network is trained by minimizing the total loss function, including: The total loss is a weighted sum of coordinate matching loss, atom type loss, chemical bond type loss, and geometric consistency loss, expressed as: In the formula, For coordinate matching loss; Atom type loss weighted by pharmacophore; For geometrically perceived loss of chemical bond types; For geometric consistency loss based on van der Waals radius; , and These are the preset weight hyperparameters for atom type loss, chemical bond type loss, and geometric consistency loss, respectively. pharmacophore-weighted atom type loss The expression is: In the formula, To measure the sampling time step Initial atom type distribution and Real Atom Type at Moments The joint expectation operator; Here, K represents the weighting coefficient corresponding to the k-th pharmacophore channel; K is the number of pharmacophore channels. The Kullback-Leibler divergence; For the forward process of flow matching, from the initial state Evolution to time The conditional probability distribution of the actual atom type at time; To generate the network in the 1st Under the condition of individual pharmacophore channels Predicted probability distribution of atomic type at time; For time steps The actual atom type label of the time ligand molecule; The pharmacophore condition vector; Geometrically Perceived Loss of Chemical Bond Type The expression is: In the formula, For time steps Initial bond order distribution and true key level The joint expectations; Category weights for different chemical bond types; This is a distance-based weighting function; For atoms and The actual Euclidean distance between them; Atoms in their initial state Interval True one-hot labels resembling chemical bonds; Atoms predicted by the network Interval Probability of chemical bonds; Let Euclidean distance be the variable between any two atoms; For atoms and The actual Euclidean distance between them; The Gaussian variance of the distance weighting function; and For atoms and 3D coordinate vector; Geometric consistency loss based on van der Waals radius The expression is: In the formula, These are the weighting coefficients for the loss term; and Atoms and atoms The van der Waals radius; This is the preset tolerance coefficient; and These are the final state atoms predicted by the network. and The coordinate vector.

[0012] Furthermore, at each time step of the flow matching sampling, a multi-objective classifier-free guidance mechanism is used to calculate the difference between the conditional vector field and the unconditional vector field of multiple drug-likeness attributes. Conflict-aware projection resolution is performed on the guidance directions of attributes with negative correlations, and the resolved attribute guidance directions are weighted and fused to obtain the fused guidance velocity field, including: The target attributes are encoded using the Gaussian radial basis function (RBF). Each target attribute is transformed and normalized according to the optimization direction. After the interval, it is mapped to a high-dimensional feature; After encoding all target attributes separately, they are mapped and concatenated into a unified attribute condition vector through a shared projection multilayer perceptron. During the training phase, Bernoulli mask variables are independently sampled for all target attributes in the attribute condition vector. Construct training conditions; During the inference generation phase, for each target attribute, the difference between the corresponding conditional guidance vector field and the unconditional vector field is calculated to obtain the independent guidance direction of the current target attribute. When the inner product of two attribute guiding vectors is negative, it indicates that there is an optimization direction conflict between them. Remove the component that is negatively correlated with the other target attribute guiding vector from the current target attribute guiding vector to obtain the projected guiding vector. At each time step of the flow matching sampling, the guiding vectors of each target attribute are weighted by a preset intensity coefficient and then superimposed onto the unconditional vector field to form the final fused guiding velocity field.

[0013] Further target attributes include: binding affinity (VinaScore), drug-likeness (QED), synthetic accessibility (SA), lipid-water partition coefficient (LogP), and topological polar surface area (TPSA).

[0014] Furthermore, the training conditions are expressed as follows: In the formula, The masked conditional vector that is actually input into the model during the training phase; For protein pocket context condition vector, it represents the fusion representation of pharmacophore condition vector and protein global features; and These are Bernoulli mask random variables for the VinaScore and TPSA attributes, respectively; the other target attributes are similarly represented. and These are the encoding vectors for the VinaScore and TPSA attributes, respectively; the other target attributes are encoded similarly. and These are the zero vector placeholders for the corresponding attributes, with dimensions consistent with the attribute encoding vector.

[0015] Furthermore, the fusion-guided velocity field drives the ordinary differential equation solver to integrate from the noisy state to the final state. The continuous variables of the final state are discretized and decoded to construct a molecular graph. After screening by chemical valence state and connectivity checks, the effective target molecules that meet the multi-objective optimization requirements are output. If the checks fail, a resampling mechanism is triggered until an effective molecule is obtained or a preset retry limit is reached, including: Based on the number of ligand atoms from the standard normal distribution The initial noise state during sampling is expressed as: In the formula, This is the initial atomic coordinate noise matrix; Let be the uniform probability tensor of the initial atom type. Number of atom type categories; To initialize the uniform probability tensor of the key type, Let be the number of key type categories, and satisfy . Symmetric constraints; It is a zero vector. It is the identity matrix; Using the initial noise state as initial values, a numerical ODE solver is employed to solve the ordinary differential equation. from Points to ; When the points reach When, the continuous final state is obtained. ;right Along the atomic type dimension, Perform the argmax operation along the bond type dimension to convert it into a discrete sequence of atom type labels and a chemical bond adjacency matrix; Preserve as a three-dimensional Cartesian coordinate matrix; Molecular graphs are constructed based on discrete atom types and chemical bond adjacency matrices. If the molecular graph satisfies the basic chemical valence state constraints and is a single connected component, the current molecular graph is used as the target molecule for the final output; otherwise, the current sample is discarded and resampled until a valid molecule is obtained or the maximum number of retries is reached.

[0016] According to specific embodiments provided by the present invention, the present invention discloses the following technical effects: This invention discloses a molecular generation method based on pharmacophore-aware priors and multi-objective flow matching. The method constructs a fused guiding velocity field, uniformly encoding pharmacophore priors and multiple druggability attribute constraints into the flow matching generation process, effectively solving the attribute conflict problem in multi-objective molecular optimization. A conflict-aware projection mechanism is employed to resolve conflicting guiding directions, ensuring that optimization objectives are synergistic rather than antagonistic. ODE integration maintains multi-objective guidance throughout the process, and combined with discretization decoding and chemical valence state verification, it significantly improves the binding affinity, drug-likeness, and synthetic accessibility of the generated molecules, while ensuring structural validity and spatial rationality. This provides an efficient and controllable molecular generation scheme for target-guided rational drug design. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 This is a schematic flowchart of the method of the present invention; Figure 2 This is a schematic diagram of the structure of the multi-channel pharmacodynamic field construction and physicochemical perception prior network in this embodiment; Figure 3 This is a schematic diagram of the pharmacophore-guided geometry-aware conditional flow matching framework in this embodiment; Figure 4 This is a schematic diagram of the multi-target classifier-less guidance mechanism in this embodiment. Detailed Implementation

[0019] 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 skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] The purpose of this invention is to provide a molecular generation method based on pharmacophore-aware prior and multi-target flow matching, aiming to solve or improve at least one of the above-mentioned technical problems.

[0021] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0022] Definition: Basic chemical valence state constraints refer to the set of atomic-level bonding rules used to verify the chemical rationality of generated molecules when decoding continuous relaxed states into discrete molecular diagrams. Its core definition is: the total number of actual bonds formed by each non-hydrogen atom in the molecular diagram (i.e., the sum of the bond orders of all chemical bonds connected to that atom) must not exceed the maximum allowed number of covalent bonds for that atom type in its current charge state. The specific parameter values ​​for the basic chemical valence state constraints in the embodiments of this application are not subjectively set, but are implemented based on common knowledge and standard tool libraries in the field of cheminformatics.

[0023] like Figure 1 As shown, this invention provides a molecular generation method based on pharmacophore-aware prior and multi-target flow matching, comprising: In this embodiment, a phased training strategy is adopted to sequentially construct a physicochemical perception prior network, a coupled SE(3)-equal variation generative network, and a multi-objective classifier-free guided capability. The training objectives at each stage are independent, and the data requirements are progressive, collectively forming a complete technical solution, including: like Figure 2 As shown, step one involves obtaining the atomic coordinates and atomic features of the target protein pocket, and using a pre-trained physicochemical perception prior network to predict the multi-channel pharmacophore distribution prior corresponding to the target protein pocket, including: The training process for the physicochemical perception prior network includes: The definition includes at least three basic channels: hydrophobic, hydrogen bond donor, and hydrogen bond acceptor. In one embodiment, channels for aromatic interactions, positive ionization centers, negative ionization centers, and metal coordination can be further expanded.

[0024] For the protein pocket of the first The atom in the first The pharmacodynamic field intensity on each channel is calculated using a Gaussian kernel function, and the expression is: In the formula, For the real scene; For the first A set of reference ligand atom indices for pharmacophore-like groups is used to provide monitoring signals for the corresponding channels, for example, when When used as a drainage channel, Includes the indices of all carbon and sulfur atoms in the reference ligand; For the corresponding to the first The preset Gaussian kernel width for each channel is used to control the decay rate of the pharmacophore's radius of action; different channels can be set with different kernel widths. The value is used to reflect the differences in their physical range of action; For the pocket number A three-dimensional Cartesian coordinate vector of an atom; For reference ligands, the first A three-dimensional Cartesian coordinate vector of atoms, where the reference ligand is a known active ligand co-crystallized with the current protein pocket in the training set; SE3-isovariable graph neural network was constructed as a physicochemical perception prior network, with the atomic coordinates and atomic features of the protein pocket as input; among them, the atomic features include embedding vectors of atomic type, residue type, partial charge, hydrophobicity identifier, etc. Train the SE3-isomorphic graph neural network to minimize the mean square error (MSE) between the pharmacodynamic field predicted by the SE3-isomorphic graph neural network and the actual field. The expression is as follows: In the formula, Mean square error; This represents the total number of pharmacophore channels. The number of atoms in the pocket; For the predicted efficacy of the medicine; This is a real-world scenario.

[0025] Through training, the SE3-isovariant graph neural network learns to predict the spatial distribution of multichannel pharmacophores directly from the pocket structure, without relying on a reference ligand.

[0026] In one embodiment, the multi-channel pharmacophore distribution prior is represented as a fixed-dimensional representation through attention pooling, ensemble pooling, or point cloud Transformer aggregation to adapt to cases where different protein pockets have different numbers of atoms, and to avoid dimensional inconsistencies caused by directly flattening the variable-length pharmacophore field.

[0027] like Figure 3 As shown, in step two, the multi-channel pharmacophore distribution prior is aggregated and fused with the global protein features to obtain the pharmacophore condition vector. This pharmacophore condition vector is then used as the input to the SE3-isovariant generation network, enabling the atomic coordinates, atomic types, and chemical bond types of the ligand molecule to be updated collaboratively in the same generation process. This includes: In this embodiment, the number of atoms of the ligand molecule to be generated It is not a fixed value; it is dynamically determined based on the current protein pocket.

[0028] The number of pocket-ligand complexes of all proteins in the training set is counted, and a conditional probability distribution of pocket atom number versus ligand atom number is constructed and fitted to a truncated normal distribution conditioned on the pocket atom number. Based on the protein pocket volume and the number of atoms in the pocket, the mode is extracted from the truncated normal distribution as the number of atoms in the ligand molecule to be generated. ; Global protein features were obtained by global pooling of pocket atom features using the SE3-equivariant encoder. The multi-channel pharmacophore distribution prior is aggregated and fused with the global features of the protein to obtain a fixed-dimensional pharmacophore condition vector, and the pharmacophore condition vector is used as the input of the generation condition to couple the SE3-equivariant generative network. Coupled SE3-equivariant generative network in time molecular state Protein pocket information, pharmacophore condition vector, and time. As input, the output is a joint prediction of the final atomic coordinates, atomic type, and chemical bond type; among which, the molecular state... This includes the current atomic coordinates, atom type, and chemical bond type; protein pocket information includes the original atomic coordinates and atomic features. Within the SE3-equivariant generative network, the first The layer update process simultaneously handles atomic features, atomic coordinates, and chemical bond features, and injects pharmacophore condition vectors into the updates of atomic and chemical bond features, expressed as: In the formula, and Atoms In the Layer and first Hidden feature vectors of the layer; and Atoms In the Layer and first The three-dimensional Cartesian coordinate vector of the layer; For atoms and The chemical bond between them is in the first The feature vector of the layer; , , These are the learnable multilayer perceptrons (MLPs) corresponding to atomic feature updates, coordinate updates, and chemical bond feature updates, respectively. For atoms The set of neighboring atoms is truncated by distance in this embodiment; The pharmacophore condition vector; The coupled SE3-isovariant generative network is trained by minimizing the total loss function. The total loss is a weighted sum of coordinate matching loss, atom type loss, chemical bond type loss, and geometric consistency loss, expressed as: In the formula, For coordinate matching loss; Atom type loss weighted by pharmacophore; For geometrically perceived loss of chemical bond types; For geometric consistency loss based on van der Waals radius; , and These are preset weight hyperparameters for atom type loss, chemical bond type loss, and geometric consistency loss, used to balance the gradients during multi-task training. In this embodiment... , , ; pharmacophore-weighted atom type loss The expression is: In the formula, To measure the sampling time step Initial atom type distribution and Real Atom Type at Moments The joint expectation operator; is the weighting coefficient corresponding to the kth pharmacophore channel, reflecting the importance of this channel in predicting the atom type; K is the number of pharmacophore channels, and K=3 when only three channels, namely hydrophobic, hydrogen bond donor, and hydrogen bond acceptor, are used. The Kullback-Leibler divergence; For the forward process of flow matching, from the initial state Evolution to time The conditional probability distribution of the actual atom type at time; To generate the network in the first Under the condition of individual pharmacophore channels, Predicted probability distribution of atom type at time step; For time steps The actual atom type label of the time ligand molecule; The pharmacophore condition vector; Geometrically Perceived Loss of Chemical Bond Type The expression is: In the formula, For time step Initial bond order distribution and true key level The negative sign before the joint expectation transforms maximizing likelihood into minimizing loss; Category weights are assigned to different chemical bond types, specifically including no bond, single bond, double bond, and triple bond; This is a distance-based weighting function; For atoms and The actual Euclidean distance between them; Initial state ( (lower atom) Interval True one-hot labels resembling chemical bonds; Atoms predicted by the network Interval Probability of chemical bonds; Let Euclidean distance be the variable between any two atoms; For atoms and The actual Euclidean distance between them; The Gaussian variance of the distance weighting function controls the tolerance range of bond length. and For atoms and 3D coordinate vector; Geometric consistency loss based on van der Waals radius The expression is: In the formula, These are the weighting coefficients for the loss term; and Atoms and atoms The van der Waals radius; A preset tolerance coefficient is used to allow slight overlap between atoms to simulate quantum effects, with a value ranging from 0.8 to 0.9. and These are the final states predicted by the network ( )atom and The coordinate vector.

[0029] In this embodiment, the AdamW optimizer is used, with a learning rate of 5×10⁻⁴, a batch size of 4, a gradient accumulation step of 8, and a training epoch of 50.

[0030] like Figure 4 As shown, in step three, at each time step of the flow matching sampling, a multi-objective classifier-free guidance mechanism is used to calculate the difference between the conditional vector field and the unconditional vector field of multiple drug-likeness attributes. Conflict-aware projection resolution is performed on the guidance directions of attributes with negative correlations, and the resolved attribute guidance directions are weighted and fused to obtain the fused guidance velocity field, including: The five target attributes of binding affinity (VinaScore), drug-likeness (QED), synthetic accessibility (SA), lipid-water partition coefficient (LogP), and topological polar surface area (TPSA) were encoded using Gaussian radial basis functions (RBF). Each target attribute is transformed and normalized according to the optimization direction. After the interval is defined, it is mapped to a high-dimensional feature, expressed as: In the formula, For target attributes After the first scalar response values ​​after Gaussian radial basis kernel mapping; To normalize to The drug-like properties of the interval; For the first The center point of a Gaussian kernel, multiple Evenly distributed in Within the interval, to achieve discretized and fine-grained encoding of continuous attributes; The preset smoothing coefficient controls the degree of overlap between adjacent Gaussian kernels; The total number of Gaussian radial basis kernels used for each target attribute; Among them, the lower the VinaScore and SA values, the better; the higher the QED value, the better; and LogP and TPSA are normalized according to the preset target interval.

[0031] After encoding all target attributes separately, they are mapped and concatenated into a unified attribute condition vector through a shared projection multilayer perceptron, expressed as: In the formula, To create a unified attribute condition vector; A multilayer perceptron projection module with shared parameters is used to align heterogeneous attribute encodings to the latent space dimension of the generative network. These are the feature vectors of VinaScore, QED, SA, LogP, and TPSA attributes encoded using RBF, respectively. , , , , These are the normalized VinaScore, QED, SA, LogP, and TPSA attribute values, respectively. During the training phase, random masks are applied independently to each attribute component of the attribute condition vector, enabling the model to learn the generation capability under arbitrary attribute subset conditions.

[0032] Independently sample Bernoulli mask variables for all target attributes in the attribute condition vector. The training conditions are constructed as follows: In the formula, The masked conditional vector that is actually input into the model during the training phase; For protein pocket context condition vector, it represents the fusion representation of pharmacophore condition vector and protein global features; and The Bernoulli mask random variables for the VinaScore and TPSA attributes are respectively (0 indicates retention, 1 indicates discard), and the other target attributes are similarly represented; and These are the encoding vectors for the VinaScore and TPSA attributes, respectively; the other target attributes are encoded similarly. and These are the zero vector placeholders for the corresponding attributes, with the same dimensions as the attribute encoding vector; During the inference generation phase, for each target attribute Calculate the difference between the corresponding conditional guidance vector field and the unconditional vector field to obtain the independent guidance direction of the current target attribute, expressed as: In the formula, Enable only attributes The difference between the conditional vector field and the unconditional vector field of the full mask represents the attribute. Independent optimization direction; To generate a network in a given molecular state ,time and the velocity field output under the given conditions; For containing only attributes Conditional concatenation vector; A conditional concatenation vector in which all attributes are masked (with all zeros as placeholders); When the inner product of two attribute guiding vectors is negative, it indicates a conflict in their optimization directions. The component negatively correlated with the other target attribute guiding vector is removed from the current target attribute guiding vector to obtain the projected guiding vector, expressed as: In the formula, Attributes after conflict-aware projection resolution Guiding vector; This is a preset numerical stability constant to prevent division by zero errors; At each time step of the flow matching sampling, the target attribute guidance vectors, after conflict resolution processing, are weighted by a preset intensity coefficient and then superimposed onto the unconditional vector field to form the final fused guidance velocity field, expressed as: In the formula, The fusion-guided velocity field used at the current time step to drive ODE / SDE sampling; The basic generated velocity field under full-attribute mask conditions; The guidance strength coefficient for each target attribute can be dynamically adjusted according to different optimization scenarios.

[0033] Step 4: Using a fused guided velocity field to drive the ordinary differential equation solver, the process integrates from the noisy state to the final state. The continuous variables of the final state are then discretized and decoded to construct a molecular graph. After verification and screening based on chemical valence state and connectivity, the effective target molecules that meet the multi-objective optimization requirements are output. If the verification fails, a resampling mechanism is triggered until an effective molecule is obtained or a preset retry limit is reached, including: Among them, the ordinary differential equation solver, in order to describe the noise distribution and recover the data distribution, generates the data along the direction of integration. Execution in a positive direction.

[0034] Based on the number of ligand atoms from the standard normal distribution The initial noise state during sampling is expressed as: In the formula, This is the initial atomic coordinate noise matrix; Let be the uniform probability tensor of the initial atom type. Number of atom type categories; To initialize the uniform probability tensor of the key type, The number of key type categories (including no-key categories) and satisfying the following conditions: Symmetric constraints; It is a zero vector. It is the identity matrix; Using the initial noise state as initial values, a numerical ODE solver is employed to solve the ordinary differential equation. from Points to ; In this embodiment, the Dormand-Prince adaptive step size method is used, and the relative tolerance is set to 1. Absolute tolerance is .

[0035] In each integral substep, the fusion-guided velocity field at the current moment is calculated. This ensures that multi-target guidance is applied throughout the entire generation trajectory.

[0036] When the points reach When, the continuous final state is obtained. ;right Along the atomic type dimension, Perform the argmax operation along the bond type dimension to convert it into a discrete sequence of atom type labels and a chemical bond adjacency matrix; Preserve as a three-dimensional Cartesian coordinate matrix; A molecular graph is constructed based on discrete atom types and chemical bond adjacency matrices; the largest connected subgraph is extracted to eliminate fragmented results, and then the valence state rules and formal charges of each atom are checked according to the built-in valence state table of RDKit. If the molecular graph satisfies the basic chemical valence state constraints and is a single connected component, the current molecular graph is used as the target molecule for the final output; otherwise, the current sample is discarded and resampled until a valid molecule is obtained or the maximum number of retries is reached. In this embodiment, the maximum number of retries is 50.

[0037] Through the above sampling and decoding process, the guiding velocity field is fused. The encoded pharmacophore specificity and multiple drug-likeness optimization directions are fully mapped into discrete molecular structures, thereby obtaining candidate molecules that simultaneously meet multiple optimization requirements such as target binding affinity and drug-likeness, synthetic accessibility, lipid-water partition coefficient, and topological polar surface area.

[0038] To verify the technical effect of the method of the present invention, 100 untrained protein pockets were selected in the CrossDocked2020 test set. 100 candidate molecules were generated for each pocket using the method of this application. Binding ability, drug-likeness, synthetic accessibility, physicochemical properties, structural validity, spatial rationality, and interaction recovery ability were statistically analyzed according to a unified molecular post-processing and evaluation process. The experimental results are shown in Table 1.

[0039] Table 1

[0040] Experimental results show that the average Top-10 Vina Score per pocket of molecules generated by the method of this application is -7.8 kcal / mol, the average QED is 0.58, and the average SA Score is 3.4; the proportion of generated molecules with LogP in the range of 1.0 to 3.5 is 70%, and the TPSA is in the range of 40 to 90 Å. 2 The proportion of the interval was 70%; the validity of the generated molecule was 88%, the clash ratio was 12%, and the H-bond recovery was 55%.

[0041] The generated conformations were scored using AutoDock Vina based on affinity, and the average Vina Score of the Top-10 molecules in each pocket was recorded. QED, SA, LogP, and TPSA were calculated using RDKit, where a lower SA score indicates easier synthesis, and the proportion of molecules falling into the preset target range was counted for LogP and TPSA.

[0042] Structural rationality was evaluated using Validity and Clash ratio, while interaction recovery was evaluated using H-bond recovery. Validity represents the proportion of molecules that passed the tests for chemical valence, connectivity, and RDKit sanitization; Clash ratio represents the proportion of generated conformations that underwent severe spatial collisions with the protein pocket; and H-bond recovery represents the proportion of generated molecules that recovered key hydrogen bond interactions with the reference ligand.

[0043] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0044] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A molecular generation method based on pharmacophore-aware prior and multi-target flow matching, characterized in that, include: Obtain the atomic coordinates and atomic features of the target protein pocket, and use a pre-trained physicochemical perception prior network to predict the multi-channel pharmacophore distribution prior corresponding to the target protein pocket; The pharmacophore distribution prior of the multi-channel pharmacophore is aggregated and fused with the global features of the protein to obtain the pharmacophore condition vector. The pharmacophore condition vector is then used as the input of the generation condition to couple the SE3-isovariant generation network, so that the atomic coordinates, atomic type and chemical bond type of the ligand molecule are updated together in the same generation process. At each time step of the flow matching sampling, the difference between the conditional vector field and the unconditional vector field of multiple drug-likeness attributes is calculated by a multi-objective classifier-free guidance mechanism. Conflict-aware projection resolution is performed on the guidance direction of attributes with negative correlations, and the resolved attribute guidance directions are weighted and fused to obtain the fused guidance velocity field. The solution of ordinary differential equations driven by the fusion guided velocity field integrates from the noisy state to the final state. The continuous variables of the final state are discretized and decoded to construct a molecular graph. After being screened by chemical valence state and connectivity verification, the effective target molecules that meet the requirements of multi-objective optimization are output. If the verification fails, a resampling mechanism is triggered until a valid molecule is obtained or the preset retry limit is reached.

2. The molecular generation method according to claim 1, characterized in that, The training of the physicochemical perception prior network includes: For the protein pocket of the first The atom in the first The pharmacodynamic field intensity on each channel is calculated using a Gaussian kernel function, and the expression is: In the formula, For the real scene; For the first A set of reference ligand atom indices for pharmacophore-like groups; For the corresponding to the first Preset Gaussian kernel width for each channel; For the pocket number A three-dimensional Cartesian coordinate vector of an atom; For reference ligands, the first A three-dimensional Cartesian coordinate vector of an atom; An SE3-isovariable graph neural network was constructed as a physicochemical perception prior network, with the atomic coordinates and atomic features of the protein pocket as input; Train the SE3-isovariable graph neural network to minimize the mean square error between the pharmacodynamic field predicted by the SE3-isovariable graph neural network and the actual field; Through training, the SE3-isovariant graph neural network learns to predict the spatial distribution of multichannel pharmacophores directly from the pocket structure.

3. The molecular generation method according to claim 2, characterized in that, The mean square error between the pharmacodynamic field predicted by the SE3-isovariant graphical neural network and the actual field is expressed as follows: In the formula, Mean squared error; This represents the total number of pharmacophore channels. Number of atoms in the pocket; For the predicted efficacy of the medicine; This is a real-world scenario.

4. The molecular generation method according to claim 1, characterized in that, The process of aggregating and fusing multi-channel pharmacophore distribution priors with global protein features to obtain a pharmacophore condition vector, and using this pharmacophore condition vector as the input to generate the SE3-isovariant generative network, includes: The number of pocket-ligand complexes of all proteins in the training set is counted, and a conditional probability distribution of pocket atom number versus ligand atom number is constructed and fitted to a truncated normal distribution conditioned on the pocket atom number. Based on the protein pocket volume and the number of atoms in the pocket, the mode is extracted from the truncated normal distribution as the number of atoms in the ligand molecule to be generated. ; Global protein features were obtained by global pooling of pocket atom features using the SE3-equivariant encoder. The multi-channel pharmacophore distribution prior is aggregated and fused with the global features of the protein to obtain a fixed-dimensional pharmacophore condition vector, and the pharmacophore condition vector is used as the input of the generation condition to couple the SE3-equivariant generative network. Coupled SE3-equivariant generative network in time molecular state Protein pocket information, pharmacophore condition vector, and time. As input, output a joint prediction of the final atomic coordinates, atomic type, and chemical bond type.

5. The molecular generation method according to claim 1, characterized in that, The atomic coordinates, atom types, and chemical bond types of the ligand molecule are updated collaboratively during the same generation process, including: Within the SE3-equivariant generative network, the first The layer update process simultaneously handles atomic features, atomic coordinates, and chemical bond features, and injects pharmacophore condition vectors into the updates of atomic and chemical bond features, expressed as: In the formula, and Atoms In the Layer and first Hidden feature vectors of the layer; and Atoms In the Layer and first The three-dimensional Cartesian coordinate vector of the layer; For atoms and The chemical bond between them is in the first The feature vector of the layer; , , These are the learnable multilayer perceptrons corresponding to atomic feature updates, coordinate updates, and chemical bond feature updates, respectively. For atoms The set of neighboring atoms; This is the pharmacophore condition vector.

6. The molecular generation method according to claim 1, characterized in that, The coupled SE3-equivariant generative network is trained by minimizing the total loss function, including: The total loss is a weighted sum of coordinate matching loss, atom type loss, chemical bond type loss, and geometric consistency loss, expressed as: In the formula, For coordinate matching loss; Atom type loss weighted by pharmacophore; For geometrically perceived loss of chemical bond types; For geometric consistency loss based on van der Waals radius; , and These are the preset weight hyperparameters for atom type loss, chemical bond type loss, and geometric consistency loss, respectively. pharmacophore-weighted atom type loss The expression is: In the formula, To measure the sampling time step Initial atom type distribution and Real Atom Type at Moments The joint expectation operator; Here, K represents the weighting coefficient corresponding to the k-th pharmacophore channel; K is the number of pharmacophore channels. The Kullback-Leibler divergence; For the forward process of flow matching, from the initial state Evolution to time The conditional probability distribution of the actual atom type at time; To generate the network in the first Under the condition of individual pharmacophore channels Predicted probability distribution of atom type at time step; For time steps The actual atom type label of the time ligand molecule; The pharmacophore condition vector; Geometrically Perceived Loss of Chemical Bond Type The expression is: In the formula, For time step Initial bond order distribution and true key level The joint expectations; Category weights for different chemical bond types; This is a distance-based weighting function; For atoms and The actual Euclidean distance between them; Atoms in their initial state Interval True one-hot labels resembling chemical bonds; Atoms predicted by the network Interval Probability of chemical bonds; Let Euclidean distance be the variable between any two atoms; For atoms and The actual Euclidean distance between them; The Gaussian variance of the distance weighting function; and For atoms and 3D coordinate vector; Geometric consistency loss based on van der Waals radius The expression is: In the formula, These are the weighting coefficients for the loss term; and Atoms and atoms The van der Waals radius; This is the preset tolerance coefficient; and These are the final state atoms predicted by the network. and The coordinate vector.

7. The molecular generation method according to claim 1, characterized in that, At each time step of the flow matching sampling, a multi-objective classifier-free guidance mechanism is used to calculate the difference between the conditional vector field and the unconditional vector field of multiple drug-likeness attributes. Conflict-aware projection resolution is performed on the guidance directions of attributes with negative correlations, and the resolved attribute guidance directions are weighted and fused to obtain a fused guidance velocity field, including: The target attributes are encoded using the Gaussian radial basis function (RBF). Each target attribute is transformed and normalized according to the optimization direction. After the interval, it is mapped to a high-dimensional feature; After encoding all target attributes separately, they are mapped and concatenated into a unified attribute condition vector through a shared projection multilayer perceptron. During the training phase, Bernoulli mask variables are independently sampled for all target attributes in the attribute condition vector. Construct training conditions; During the inference generation phase, for each target attribute, the difference between the corresponding conditional guidance vector field and the unconditional vector field is calculated to obtain the independent guidance direction of the current target attribute. When the inner product of two attribute guiding vectors is negative, it indicates that there is an optimization direction conflict between them. Remove the component that is negatively correlated with the other target attribute guiding vector from the current target attribute guiding vector to obtain the projected guiding vector. At each time step of the flow matching sampling, the guiding vectors of each target attribute are weighted by a preset intensity coefficient and then superimposed onto the unconditional vector field to form the final fused guiding velocity field.

8. The molecular generation method according to claim 7, characterized in that, The target attributes include: binding affinity (VinaScore), drug-likeness (QED), synthetic accessibility (SA), lipid-water partition coefficient (LogP), and topological polar surface area (TPSA).

9. The molecular generation method according to claim 8, characterized in that, The training condition is expressed as follows: In the formula, The masked conditional vector that is actually input into the model during the training phase; For protein pocket context condition vector, it represents the fusion representation of pharmacophore condition vector and protein global features; and These are Bernoulli mask random variables for the VinaScore and TPSA attributes, respectively; the other target attributes are similarly represented. and These are the encoding vectors for the VinaScore and TPSA attributes, respectively; the other target attributes are encoded similarly. and These are the zero vector placeholders for the corresponding attributes, with dimensions consistent with the attribute encoding vector.

10. The molecular generation method according to claim 1, characterized in that, The method utilizes a fusion-guided velocity field to drive an ordinary differential equation solver to integrate from the noisy state to the final state, performs discretization decoding on the continuous variables of the final state to construct a molecular graph, and outputs effective target molecules that meet the requirements of multi-objective optimization after chemical valence state and connectivity verification. If the verification fails, a resampling mechanism is triggered until a valid molecule is obtained or the preset retry limit is reached, including: Based on the number of ligand atoms from the standard normal distribution The initial noise state during sampling is expressed as: In the formula, This is the initial atomic coordinate noise matrix; Let be the uniform probability tensor of the initial atom type. Number of atom type categories; To initialize the uniform probability tensor of the key type, Let be the number of key type categories, and satisfy . Symmetric constraints; It is a zero vector. It is the identity matrix; Using the initial noise state as initial values, a numerical ODE solver is employed to solve the ordinary differential equation. from Points to ; When the points reach At that time, the continuous final state is obtained. ;right Along the atomic type dimension, Perform the argmax operation along the bond type dimension to convert it into a discrete sequence of atom type labels and a chemical bond adjacency matrix; Preserve as a three-dimensional Cartesian coordinate matrix; Molecular graphs are constructed based on discrete atom types and chemical bond adjacency matrices. If the molecular graph satisfies the basic chemical valence state constraints and is a single connected component, the current molecular graph is used as the target molecule for the final output; otherwise, the current sample is discarded and resampled until a valid molecule is obtained or the maximum number of retries is reached.

Images