Drug molecular structure three-dimensional modeling and optimization processing method

Abstract

The invention relates to the technical field of pharmaceutical chemistry and bioinformatics, in particular to a drug molecular structure three-dimensional modeling and optimization processing method, which comprises the following steps: constructing a three-dimensional virtual space containing a pathogen and a drug agent; setting a connection logic based on distance and energy, so that the drug agents are self-assembled to form molecular clusters when meeting conditions; dynamic simulation is operated, and molecular clusters are driven to form space surrounding on the surface of the pathogen; calculating the surface coverage rate of the cluster to the pathogen key sites and generating a fitness score; and carrying out updating such as hydrophobic modification and the like on the drug molecule topological structure based on the score by utilizing an evolutionary algorithm, and iteratively generating an optimal structure. According to the method, a multi-agent cooperation mechanism is introduced into a calculation model, and drug self-assembly and virus space encirclement are incorporated into screening indexes, so that the method can discover a novel drug structure for physically blocking virus infection by forming a large-volume cluster, and the problem of screening leakage in a traditional method is effectively solved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of pharmaceutical chemistry and bioinformatics, in particular to a drug molecule structure three-dimensional modeling and optimization processing method. BACKGROUND

[0002] In modern drug research and development, computer-aided drug design (CADD) has become a key technology to shorten the development cycle and reduce costs. The existing mainstream technology mainly includes structure-based drug design (SBDD), and its typical representatives are molecular docking (Molecular Docking) and molecular dynamics simulation (MD). The existing technical solutions usually have the following process: obtaining the static crystal structure of the virus target protein, using the scoring function to calculate the binding free energy of a single drug small molecule ligand and the target active pocket, and screening potential drugs, for example, AutoDock or Schrödinger software adopts such a "lock and key model" for screening.

[0003] However, the existing technology has the following obvious defects: Firstly, the swarm effect is ignored, and the existing method is mostly based on the "one-to-one" binding mode for prediction. However, in the high-concentration local drug delivery scene in the real physiological environment, drug molecules often self-assemble into dimers or clusters through intermolecular weak interactions (such as van der Waals force, Piling up), this "molecular cluster" can wrap the virus surface through a larger steric hindrance volume, producing an inhibitory effect that monomers do not have, and the existing technology cannot simulate this dynamic assembly process, resulting in the screening of drugs that have weak monomer binding but strong group blocking ability; Secondly, the static optimization path is single, and the traditional quantitative structure-activity relationship (QSAR) analysis often relies on manual experience to modify functional groups, and lacks a closed-loop mechanism that can automatically evolve molecular structures according to the "space surrounding ability", resulting in low efficiency of structure optimization.

[0004] Therefore, a drug molecule structure three-dimensional modeling and optimization processing method is proposed. SUMMARY

[0005] The purpose of the present application is to provide a drug molecule structure three-dimensional modeling and optimization processing method to solve the technical problems that the existing drug design method ignores the synergistic surrounding effect of drug molecules forming clusters in the real environment, and lacks an automatic structure evolution mechanism based on the steric hindrance blocking effect.

[0006] To solve the above technical problems, the purpose of the present application is to provide a drug molecule structure three-dimensional modeling and optimization processing method, comprising the following steps: S1, constructing a three-dimensional virtual space in a computing device, mapping a target pathogen as a target object, and mapping an initial candidate drug molecule as a drug agent object with autonomous motion attributes; S2, setting a drug molecule interconnection logic, the connection logic being used to monitor the relative states between the drug agent objects in real time, and establishing a virtual rigid constraint to lock a plurality of drug agent objects as a co-moving molecular cluster object when a preset combination condition is met; S3, running a dynamics simulation in the three-dimensional virtual space, controlling the position update of the drug agent objects or the molecular cluster object based on a physical potential field, and driving the formation of a spatial enclosure of the target object by the drug agent objects or the molecular cluster object; S4, calculating the surface coverage of the target object key site by the molecular cluster object, and generating a fitness score of the current drug molecule based on the surface coverage; S5, based on the fitness score, updating the topological structure of the drug molecule using an evolutionary algorithm, generating a new generation of drug agent objects and returning to step S2 until a preset convergence condition is met, and outputting an optimal drug molecule structure.

[0007] As a further improvement of the technical solution, in step S2, the preset combination condition specifically includes: extracting the charge distribution data of the contact surface, and calculating the electrostatic potential product of the contact area; When the Euclidean distance between two drug agent objects is less than the van der Waals contact threshold, and the electrostatic potential product is negative and meets the energy threshold requirement, it is determined that the preset combination condition is met.

[0008] As a further improvement of the technical solution, in step S2, the virtual rigid constraint is established, specifically including: defining the first drug agent object and the second drug agent object as a rigid combination in the physical engine, locking the relative position and relative angle of the two objects relative to the common center of mass, so that the force applied to the rigid combination in the subsequent simulation produces a whole translation or rotation motion.

[0009] As a further improvement of the technical solution, in step S2, the connection logic further includes a rupture judgment: In each simulation frame, the external shear stress borne by the virtual rigid constraint inside the molecular cluster object is calculated, and the corresponding equivalent deformation energy is calculated in combination with the geometric parameters of the constraint scope; When the equivalent deformation energy is greater than the binding energy threshold, the virtual rigid constraint is released, and the molecular cluster object is decomposed into independent drug agent objects.

[0010] As a further improvement of the technical solution, in step S3, the running dynamics simulation specifically includes: For each simulation frame, the force vector acting on each object is calculated, which at least includes an electrostatic attraction component pointing to the target object and a random disturbance component simulating thermal motion; The spatial coordinates and attitude of the object are updated according to the force vector using a numerical integration algorithm.

[0011] As a further improvement of the technical solution, in step S3, the position update of the drug agent object or the molecular cluster object also includes attitude adjustment logic: When it is detected that the distance between the drug agent object or the molecular cluster object and the target object is less than a preset sensing threshold, the electrostatic torque acting on it is calculated based on the electric dipole moment of the drug agent object or the molecular cluster object and the potential gradient of the surface of the target object; The torque is applied to the drug agent object or the molecular cluster object to perform a rotation transformation, so that its binding site is oriented towards the target object.

[0012] As a further improvement of the technical solution, in step S1, the initialization simulation space is also constructed to include: A plurality of environmental disturbance particles are generated in the three-dimensional virtual space, which have a preset mass attribute and an initial velocity vector; During the simulation process of step S3, when the coordinates of the environmental disturbance particles coincide with the coordinates of the drug agent or the molecular cluster object, the velocity vector after collision is calculated according to the law of conservation of momentum, thereby changing the motion trajectory of the drug agent or the molecular cluster object.

[0013] As a further improvement of the technical solution, in step S4, the generation logic of the fitness score is: The fitness score is calculated by weighted linear combination; The variables of the weighted linear function at least include the surface coverage, the total bond energy of all virtual rigid constraints inside the molecular cluster object, and the number of drug agents constituting the molecular cluster object.

[0014] As a further improvement of the technical solution, in step S5, the topology of the drug molecule is updated using an evolutionary algorithm, specifically including: The topology of the high fitness score molecule is converted into a molecular graph; Non-skeletal nodes in the molecular graph are identified, which are defined as side chain atoms or terminal atoms that do not belong to the molecular pharmacophore skeleton; Randomly select one of the non-backbone nodes, replace it with an atom or functional group fragment with a higher hydrophobicity parameter to enhance the intermolecular binding capacity, and generate the topological structure data of the offspring molecules.

[0015] Compared with the prior art, the present application has the following beneficial effects: 1. In the drug molecule structure three-dimensional modeling and optimization processing method, the multi-agent collaborative mechanism is introduced into the calculation model, the drug self-assembly and virus space surrounding are included in the screening index, so that the method can find new drug structures that physically block virus infection by forming large volume clusters, effectively solving the problem of missing screening in traditional methods.

[0016] 2. In the drug molecule structure three-dimensional modeling and optimization processing method, an automatic closed loop from drug efficacy evaluation (coverage) to structure modification (topology update) is established, which can automatically evolve more beneficial to group and combination of molecular structure according to the effect of collaborative surrounding, thereby improving the efficiency of drug optimization.

[0017] 3. In the drug molecule structure three-dimensional modeling and optimization processing method, by introducing environmental disturbance particles and momentum conservation calculation, the Brownian motion disturbance in the real body fluid environment is simulated, so that the prediction result has higher robustness and clinical reference value.

[0018] 4. In the drug molecule structure three-dimensional modeling and optimization processing method, a weighted linear function based on coverage, bond energy and quantity is used as fitness score, which not only ensures the pursuit of drug efficacy, but also considers the dosage cost and binding stability of the drug, providing multi-dimensional optimization direction. BRIEF DESCRIPTION OF DRAWINGS

[0019] Figure 1 The flowchart of the present application. DETAILED DESCRIPTION

[0020] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0021] Example 1 At present, the existing technical solutions usually have the following process: obtaining the static crystal structure of the viral target protein, calculating the binding free energy of a single drug small molecule ligand and the target active pocket by using a scoring function, and screening potential drugs, for example, the software such as AutoDock or Schrödinger adopts such a "lock and key model" for screening; However, the existing drug design method has the problems of ignoring the synergistic surrounding effect of drug molecules after forming clusters in the real environment, and lacking an automatic structure evolution mechanism based on the space steric hindrance blocking effect.

[0022] In view of this, as Figure 1 The purpose of the embodiment is to provide a drug molecule structure three-dimensional modeling and optimization processing method, taking the optimization of small molecule inhibitors of the spike protein of the new coronavirus as an application scenario, which comprises the following steps: Since an approximate virtual space similar to the real physiological environment is constructed, it provides a basis for subsequent simulation of the interaction between drug molecules and pathogens; mapping the entity molecules into virtual objects facilitates dynamic control and monitoring through computer algorithms, therefore, in step S1, a three-dimensional virtual space is constructed in a computing device, a target pathogen is mapped into a target object, and an initial candidate drug molecule is mapped into a drug intelligent agent object with autonomous motion attributes; The specific implementation manner is as follows: Firstly, a workstation is selected, which is equipped with an Intel Xeon W-3275 processor (28 cores and 56 threads), an NVIDIA A100 GPU (40GB of video memory), and an Ubuntu 22.04 LTS operating system; The virtual space construction software adopts Unity 3D 2022.3 version, and the space coordinate system is built by using the physical engine module (Cartesian coordinate system, unit: angstrom Å); Secondly, the target pathogen is selected as the spike protein of the new coronavirus (S protein), the three-dimensional structure data of which is obtained from the PDB database (number 6VSB), and the structure is pretreated (water molecules and ligands are removed, and the hydrogen bond network is optimized) by using PyMOL2.5 software; The pretreated structure is imported into Unity3D and mapped into a "target object"; the size ratio of the target object in the virtual space is set to 1:1000 (that is, 1 Å in reality corresponds to 1000 unit lengths in the virtual space), and the key sites (ACE2 binding sites, including 20 amino acid residues such as Lys417, Tyr453, His475, which are determined by the S protein-ACE2 binding structure reported in the literature [Nature, 2020]) are marked; The initial candidate drug molecule is a known small molecule inhibitor of the S protein (e.g., the SMILES string is “C1=CC=C(C=C1)CN2C(=O)C3=C(C(=O)N(C2=O)C)CCCC3”). RDKit version 2023.03 is used to generate a three-dimensional structure based on the SMILES string, and its atomic charges (using the AM1-BCC charge calculation method) and van der Waals radii (C atom 1.7Å, N atom 1.55Å, O atom 1.52Å) are calculated. This three-dimensional structure is then imported into a virtual space and mapped as a “drug agent object”, endowed with autonomous motion attributes: the initial velocity is set to 0.5Å / ps (based on the Brownian motion velocity of small molecules at physiological temperature 37℃), and the direction of motion is random. A C# script is written to implement the motion rule of “updating the direction of motion every 1ps, with a deviation range of ±15°”.

[0023] The three-dimensional virtual space constructed using the above-mentioned technical means can accurately reproduce the structure and key sites of the S protein. The autonomous motion properties of the drug agent conform to the motion law of actual small molecules, laying a reliable foundation for subsequent dynamic binding and envelopment simulation.

[0024] Because the simulation of specific binding between multiple drug molecules forms synergistic molecular clusters, it enhances the coverage of key sites of the target object; at the same time, rigid constraints ensure the consistency of cluster movement and avoid cluster dispersion leading to a decrease in binding efficiency. Therefore, in step S2, drug molecule connection logic is set. This connection logic is used to monitor the relative state between drug agent objects in real time and establish virtual rigid constraints when preset binding conditions are met, locking multiple drug agent objects into molecular cluster objects that move together. The specific implementation method is as follows: A real-time monitoring script was written in C++ and integrated into the Update function of Unity 3D (monitoring frequency is 100 times / second, i.e., acquiring the state of the drug agent every 10ms). The monitoring parameters include the Euclidean distance between the two drug agents and the charge distribution data of the contact area; the Euclidean distance is calculated using the formula for the distance between two points. x1, y1, z1 and x2, y2, z2 are the three-dimensional coordinates of the centroids of the two smart agents respectively. The charge distribution data is used to calculate the electrostatic potential (φ1, φ2) on the surface of each smart agent through the Amber 20 force field. The average electrostatic potential of the contact area (the area where the distance is less than the sum of the van der Waals radii) is taken as the basis for calculation. The preset combination conditions include the following: The sum of the van der Waals radii of the surface atoms of the two agents is determined by the van der Waals contact threshold value according to the atomic type of the drug agent, for example, the sum of the van der Waals radii of C and O atoms is 1.7+1.52=3.22Å, so the van der Waals threshold value of this type of contact is set to 3.2Å (considering the virtual space error, take the approximate value); The electrostatic potential product and energy threshold value are determined. Through the pre-experiment (using AutoDock Vina to calculate the binding energy of drug molecules under different electrostatic potential products), it is found that when the binding energy is lower than-5kcal / mol, the molecular binding stability meets the demand, and the electrostatic potential product (φ1xφ2) is-0.2eV 2 (Since φ1 and φ2 are of opposite signs, electrostatic attraction is generated, so the product is negative), so the "Euclidean distance < 3.2Å and electrostatic potential product ≤-0.2eV 2 " is set as the preset binding condition; The virtual rigid constraint is established, and the PhysX 5.1 physics engine is used to realize the constraint. When the binding condition is met, the "FixedJoint" function of the engine is called to define the two drug agents as a rigid combination: lock the relative position of the two agents relative to the common center of mass (the center of mass distance is fixed as the Euclidean distance when binding, such as 3.1Å), lock the relative angle (based on the optimal binding angle calculated by AutoDock Vina, such as the angle between C-O bond and the direction of hydrogen bond of the target site is 15°), and the force (such as electrostatic attraction) applied to the rigid combination in subsequent simulation will act on the two agents simultaneously, realizing the overall translation or rotation; Fracture judgment, in each simulation frame, calculate the external shear stress borne by the internal virtual rigid constraint of the molecular cluster object; according to the shear stress and contact area, calculate the work or energy accumulation of the constraint region, when the calculated energy value is greater than the binding energy threshold, the virtual rigid constraint is released, and the molecular cluster object is decomposed into independent drug agent objects, specifically: set every 10ps as a simulation frame, calculate the external shear stress borne by the internal virtual rigid constraint of the molecular cluster at the end of each frame: obtain the external force F (unit: kcal / (mol.Å)) at the constraint through the PhysX engine in real time, the contact area A (calculated according to the number of atoms in the binding region, the contact area of each atom is taken as 0.2Å 2 , such as A=0.6Å 2 when 3 atoms are in contact), calculate the current stress intensity according to the shear stress formula σ=F / A; the binding energy threshold is set to-5kcal / mol (consistent with the energy threshold of the binding condition), then the stress is converted into energy value for comparison: the calculation formula is (where, is the contact area, The characteristic length of the virtual constraint (representing the bond length or action distance), when calculated... When the energy is greater than 5 kcal / mol (i.e. the binding energy barrier has been overcome), the "BreakJoint" function is called to remove the virtual rigid constraint and decompose the molecular cluster into independent drug agents. Through the above technical means, the connection logic can accurately identify the binding conditions of drug intelligence, the rigid constraint ensures the stable movement of molecular clusters, and the fracture judgment simulates the situation of clusters breaking due to external forces in the real environment, avoiding invalid clusters (such as unstable clusters) from affecting subsequent optimization, and improving the effectiveness of clusters and the realism of simulation.

[0025] Because the movement of drug agents or clusters is controlled by a physical potential energy field, driving them to approach key sites on the target object and form an enclosure, the binding process of drug molecules and pathogens under real physiological conditions is simulated; at the same time, environmental interference is introduced to further enhance the realism of the simulation. Therefore, in S3, a dynamic simulation is run in the three-dimensional virtual space, and the position of the drug agent object or the molecular cluster object is updated based on the physical potential energy field, driving it to form a spatial enclosure around the surface of the target object; the specific implementation method is as follows: The construction of a physical potential field specifically includes: The physical potential field includes an electrostatic attraction component and a random perturbation component. For each drug agent or cluster, the net force vector F_net = F_static + F_perturbation is calculated. The electrostatic attraction component F_static is calculated based on Coulomb's law: F_static = k × q_1 × q_2 / r 2 Where k = 9e9N.m 2 / C 2 (Coulomb constant), q1 is the total charge of the drug agent / cluster (obtained by summing the AM1-BCC charges calculated by RDKit, e.g., the total charge of the molecule is +1e), q2 is the total charge of the key site of the target object (the total charge of the ACE2 binding site is -2e, calculated by the Amber force field), r is the distance between the centroid of the drug agent / cluster and the centroid of the key site; the direction of the force is from the drug agent / cluster to the key site (due to the opposite signs of q1 and q2, an attractive force is generated); the random perturbation component F_perturbation simulates the thermal motion in the physiological environment, generated based on a Gaussian distribution, with a mean μ=0 and a standard deviation σ=0.1kcal / (mol·Å) (calculated using the Boltzmann distribution at 37℃ to ensure that the perturbation intensity matches the thermal motion energy of small molecules), the direction is random, and it is updated every 1ps; the position update uses the Velocity Verlet numerical integration algorithm to update the spatial coordinates and velocity of the object, with the integration time step set to 1fs (determined according to the conventional time step of molecular dynamics simulation to balance calculation accuracy and efficiency), and the coordinate update formula is: The speed update formula is: where, (m is the mass of the drug agent / cluster, summed by atomic mass, e.g. the mass of this molecule is 250 u), is the 3D coordinates of the center of mass of the drug agent (or molecular cluster) at time t (unit: Å, Angstrom); is the solution of the 3D coordinates of the center of mass at time t+1 (i.e. the next time step); is the instantaneous velocity of the drug agent at time t (unit: Å / ps, Angstrom / picosecond), determined by the autonomous motion attribute of the drug agent and the resultant force; is the integral time step (in this invention, it is fixed as This value is chosen because the atomic vibration period of small molecule drugs is about 10 fs, less than 1 / 10 of the vibration period, in order to ensure the simulation accuracy and avoid distortion of the motion trajectory); The construction of the initialized simulation space also includes: generating a plurality of environmental disturbance particles in the three-dimensional virtual space, the environmental disturbance particles having a preset mass attribute and initial velocity vector; during the simulation process of step S3, when the coordinates of the environmental disturbance particles coincide with the coordinates of the drug agent or the molecular cluster object, the velocity vector after collision is calculated according to the law of conservation of momentum, thereby changing the motion trajectory of the drug agent or the molecular cluster object, specifically: generating H2O molecules as environmental disturbance particles in the three-dimensional virtual space, the number of which is 100 times the number of drug agents (to simulate the proportion of solvents in the physiological environment), the mass attribute of each disturbance particle is set to 18u (the molar mass of H2O molecule), and the initial velocity vector is based on the thermal motion velocity at 37℃ (about 500m / s, the direction is random, calculated by Maxwell velocity distribution); during the dynamics simulation, when the coordinate distance between the disturbance particle and the drug agent is <0.5Å (determined as collision), the velocity after collision is calculated according to the law of conservation of momentum: where , and , are the mass and pre-collision velocity of the drug agent, respectively, , are the mass and pre-collision velocity of the disturbance particle, respectively, The control of the drug agent object or the molecular cluster object for position updating further includes pose adjustment logic: when it is detected that the distance between the drug agent object or the molecular cluster object and the target object is less than a preset sensing threshold, an electrostatic torque acting on the drug agent object or the molecular cluster object is calculated based on the electric dipole moment of the drug agent object or the molecular cluster object and the potential gradient of the surface of the target object; a rotation transformation is performed on the drug agent object or the molecular cluster object by applying the torque, so that the binding site thereof faces the target object; specifically: The sensing threshold is set to 5 Å (determined based on the effective binding distance of the drug agent and the target site), and when the distance between the drug agent / cluster and the target object is < 5 Å, the torque relative to the potential gradient of the surface of the target object is calculated: first, the potential distribution of the surface of the target object is calculated by the Poisson-Boltzmann equation (x, y, z), the potential gradient is calculated, wherein the gradient direction points to the direction in which the potential increases, refers to the electrostatic potential (unit: eV, electron volt) of a point in a three-dimensional virtual space, is the potential with unchanged y and z coordinates, the rate of change of the potential along the x-axis direction (unit: eV / Å), if the value is positive, it means that the potential gradually increases along the positive direction of the x-axis; if it is negative, it means that the potential gradually decreases along the positive direction of the x-axis, respectively correspond to the potential change rates along the y-axis and z-axis directions, and the physical meaning is the same as ; then the gradient force (q is the charge of the drug agent) and the torque (r is the vector from the centroid of the drug agent to the centroid of the key site of the target object) are calculated; finally, a rotation transformation is performed on the drug agent / cluster by using a quaternion rotation algorithm, the rotation axis is the direction of the torque M, and the rotation angle is the angle corresponding to the size of the torque (until the binding site of the drug agent, such as the hydrogen bond donor site of the N-H bond, faces the hydrogen bond acceptor site of the key site of the target object, such as the O atom); Through the above technical means, the force control of the physical potential field enables the drug agent / cluster to continuously approach the target site, the random disturbance and the environmental interference simulate the uncertainty of the real physiological environment, and the pose adjustment improves the alignment accuracy of the binding site. The three work together to realize efficient spatial surrounding of the drug agent / cluster on the surface of the target object, and improve the initial surrounding rate of the key site.

[0026] In view of the quantitative evaluation of the binding effect of the molecular cluster on the key site of the target object, the stability of the cluster and the size of the molecule are also considered to avoid the optimization deviation caused by a single index, and a reliable evaluation basis is provided for subsequent evolutionary optimization. Therefore, in step S4, the surface coverage of the molecular cluster object on the key site of the target object is calculated, and a fitness score of the current drug molecule is generated based on the surface coverage; the specific implementation is as follows: Surface coverage calculation: Key site surface area determination: the surface area S0 of the key site (ACE2 binding site) of the target object is calculated using VMD 1.9.4 software, and the S0=200Å is obtained by calculating the "Accessible Surface Area" tool. 2 ; Cluster coverage area calculation: by using the ray detection function of Unity 3D, a ray is emitted from each atom on the surface of the key site of the target object to the three-dimensional space (the number of rays is 100 per atom, covering all directions), and when the ray is blocked by the atom of the molecular cluster (the length of the ray is <0.5Å, which is determined as coverage), the surface area element corresponding to the ray (0.01Å 2 corresponding to each ray is recorded, and the key area covered by the cluster S1 is obtained by summation; Surface coverage calculation: coverage C=(S1 / S0)×100%, for example, when S1=130Å 2 , C=65%; Fitness score generation: the fitness score F is calculated by using a weighted linear combination formula, and the formula is wherein a, b, and c are weight coefficients, which are determined by the analytic hierarchy process (AHP): five experts in the field of drug design are invited to score the importance of "binding effect (C)", "cluster stability (E)", and "molecular size (N)", a judgment matrix is constructed, and the weights are calculated, and finally a=0.5 (the binding effect is the most important), b=0.3 (the stability is the second), and c=0.2 (the molecular size needs to be controlled to avoid the decrease of the cell membrane penetration caused by the too large size); E is the total energy of all virtual rigid constraints in the molecular cluster: the energy of each virtual rigid constraint is set to -5kcal / mol (the energy when the binding is stable) based on the binding energy, and when the cluster contains constraints, E= , for example, =3, E=-15kcal / mol, for the convenience of calculation, E is normalized (E , wherein =-50kcal / mol (the maximum number of constraints is 10), =0 kcal / mol (unconstrained), E normalized to 0-1); N is the number of drug agents constituting the molecular cluster: set a reasonable range for N to be 1-10 (determined according to the size of the target site; too many will result in excessively large molecules), and normalize N ( , =10, =1, N normalization range is 0-1, the smaller the quantity, the higher the normalization value, and moderate molecule size is encouraged); for example, C=65% (i.e. 0.65), E normalization=0.7 (when E=-15kcal / mol, (-15-(-50)) / (0-(-50))=35 / 50=0.7), N normalization=0.8 (when N=3, (If we take 0.8), then .

[0027] Using the above-mentioned techniques, surface coverage can intuitively reflect the binding effect of clusters. The fitness score of weighted linear combination integrates binding ability, stability and molecular size, avoiding the problem of large but unstable clusters or small clusters with insufficient coverage caused by relying solely on coverage. The score results can guide the direction of subsequent evolution and optimization.

[0028] Since high-quality drug molecules are selected based on fitness scores, and better progeny molecules are generated through topological structure updates, the fitness of the molecules is gradually improved, achieving efficient search for the globally optimal drug molecule structure. Therefore, in step S5, based on the fitness scores, the topological structure of the drug molecules is updated using an evolutionary algorithm to generate a new generation of drug agent objects and return to step S2 until the preset convergence condition is met, outputting the optimal drug molecule structure; the specific implementation method is as follows: The evolutionary algorithm parameter settings include the following: A genetic algorithm was used as the evolutionary algorithm, with a population size of 50 (balancing computational efficiency and diversity), and 100 generations (preliminary experiments showed that the scores tended to stabilize after 100 generations). The selection operator was roulette wheel selection (molecules with higher fitness scores had a higher probability of being selected). , For the first (Fitness score of each molecule), the crossover operator probability is set to 0.7, and the mutation operator probability is set to 0.1; The evolutionary algorithm is used to update the topology of drug molecules, specifically including: converting the topology of high-fitness-score molecules into a molecular graph; identifying non-backbone nodes in the molecular graph, defined as side-chain atoms or terminal atoms that do not belong to the molecular pharmacophore backbone; randomly selecting one of the non-backbone nodes and replacing it with an atom or functional group fragment with a higher hydrophobicity parameter to enhance the binding ability between molecules, generating topology data of progeny molecules: specifically: The topology of the molecules with high fitness scores (select the top 20% of molecules in each generation, i.e. 10 molecules) is converted into a molecular graph, the molecular graph is constructed using the NetworkX 3.1 library, the nodes in the graph are the atoms of the molecules, and the edges are the chemical bonds between the atoms (single, double, and triple bonds are marked as 1, 2, and 3, respectively), thereby completing the molecular graph conversion; The "molecular pharmacophore skeleton" is defined as the core structure that maintains the pharmacophore in the molecule (such as the benzene ring and the amide bond "-CONH-" in this embodiment), and the non-skeleton nodes are side chain atoms or terminal atoms that do not belong to the skeleton, such as the terminal H atom of the methyl group (-CH3) and the H atom in the amino group (-NH2), thereby completing the identification of non-skeleton nodes; A non-skeleton node (such as the H atom in the methyl group) is randomly selected and replaced with an atom or functional group fragment with a higher hydrophobicity parameter: the hydrophobicity parameter uses the logP value (n-octanol-water partition coefficient, the higher the logP value, the stronger the hydrophobicity), the logP value of the original H atom is -0.5, and it is replaced with a Cl atom (logP value +0.3) or a trifluoromethyl group (-CF3, logP value +1.0); During the replacement process, the chemical bond must be reasonable (such as the H atom being a single bond, the replaced Cl atom still being a single bond, and avoiding the occurrence of unreasonable structures such as a pentavalent C atom), and the three-dimensional structure of the molecule is regenerated using RDKit after replacement, and the new charge, van der Waals radius, and other parameters are calculated as a new generation of drug agent object, thereby completing the replacement of non-skeleton nodes; The convergence condition judgment includes the following contents: Two convergence conditions are set, and any one condition is met to stop iteration: The change rate of the average fitness score of the population of the last 10 generations is less than 1% (i.e. , The average score of the nth generation), indicating that the score tends to be stable; There is a molecule in a certain generation with a fitness score ≥0.9 (preliminary experiments show that molecules with a score ≥0.9 have met the actual application requirements); When the convergence condition is met, the three-dimensional structure data (including atomic coordinates, bond types, and charge distribution) of the molecule with the highest fitness score in the population is output, which can be exported in PDB format for subsequent experimental synthesis and activity verification, thereby completing the output of the optimal molecule; Through the above technical means, the evolutionary algorithm can quickly improve the fitness score of the molecule through the cycle of selection, update, and screening, thereby improving the fitness score, surface coverage, and hydrophobicity parameter of the optimal molecule (enhancing cell membrane penetration), and achieving the optimization of the structure of the drug molecule.

[0029] Example 2 This example is a variant of Example 1, with the application scenario of penicillin binding protein of E. coli (PBP2a) optimized as the application scenario, the core steps are consistent with Example 1, and only the key parameters are adjusted according to the characteristics of the target pathogen and the drug molecule to verify the universality of the method: Step S1: the target pathogen is penicillin binding protein of E. coli (PBP2a, PDB ID 1VQQ), the key site is the region around the serine residue (Ser403) in its active center; the initial candidate drug molecule is penicillin G (SMILES string is “CC1(C)S[C@@H]2C@HC(=O)N2C1=O”), and the initial speed of the drug agent is set to 0.4 Å / ps (the molecular weight of penicillin G is larger than that of the small molecule in Example 1, and the speed is appropriately reduced); Step S2: The van der Waals contact threshold is adjusted according to the atom type of penicillin G, such as the threshold of S atom (van der Waals radius 1.85 Å) and O atom is set to 3.37 Å (1.85+1.52), and the electrostatic potential product threshold is set to -0.18 eV 2 (the binding energy threshold of penicillin G is -4.5 kcal / mol); Step S3: In the electrostatic attraction component, the total charge of the key site of PBP2a is -1e, and the standard deviation of the random disturbance component is set to 0.08 kcal / (mol.Å) (the molecular weight is large, and the thermal motion disturbance is reduced); in addition to H2O, a small amount of Mg 2 + ions (simulate the intracellular environment of bacteria) are added, the mass attribute is set to 24 u, and the initial speed is set to 400 m / s; Step S4: The surface area of the key site S0=150 Å 2 , the fitness score weight adjustment is a=0.55 (PBP2a key site is small, coverage is more important), b=0.25, c=0.2; Step S5: The non-backbone node is the terminal atom of the phenylacetyl side chain of penicillin G, which is replaced by naphthylacetyl (-COCH2C10H7, logP value is higher), and after 80 generations of evolution algorithm iteration, the convergence condition is met; This example successfully applies the method to antibacterial drug optimization through parameter adjustment, proving the universality of the method; the binding ability of the optimized drug molecule to PBP2a is significantly improved, providing a new technical path for drug research and development of drug-resistant bacteria.

[0030] The above shows and describes the basic principles, main features and advantages of the present application. Those skilled in the art should understand that the present application is not limited to the above-mentioned embodiments, and the above-mentioned embodiments and descriptions in the specification are only preferred examples of the present application and are not intended to limit the present application. Various changes and improvements can be made to the present application without departing from the spirit and scope of the present application, and these changes and improvements all fall within the scope of the claimed present application. The scope of protection of the present application is defined by the appended claims and their equivalents.

Claims

1. A method for three-dimensional modeling and optimization of drug molecule structures, characterized in that: Includes the following steps: S1. Construct a three-dimensional virtual space in the computing device, map the target pathogen as the target object, and map the initial candidate drug molecules as drug intelligent objects with autonomous movement attributes; S2. Set up drug molecule connection logic. The connection logic is used to monitor the relative state between drug intelligent objects in real time, and establish virtual rigid constraints when the preset binding conditions are met, locking multiple drug intelligent objects into molecular cluster objects that move together. S3. Run a dynamic simulation in the three-dimensional virtual space, control the position update of the drug intelligence object or the molecular cluster object based on the physical potential energy field, and drive it to form a spatial enclosure around the surface of the target object. S4. Calculate the surface coverage of the molecular cluster object to the key sites of the target object, and generate the fitness score of the current drug molecule based on the surface coverage. S5. Based on the fitness score, the topological structure of the drug molecule is updated using an evolutionary algorithm to generate a new generation of drug agent objects and return to step S2 until the preset convergence condition is met, and the optimal drug molecule structure is output.

2. The method for three-dimensional modeling and optimization of drug molecule structure according to claim 1, characterized in that: In step S2, the preset combination conditions specifically include: Extract the charge distribution data of the contact surface and calculate the electrostatic potential product of the contact area; When the Euclidean distance between two drug-initiated intelligent objects is less than the van der Waals contact threshold, and the electrostatic potential product is negative and meets the energy threshold requirement, the preset binding condition is determined to be satisfied.

3. The method for three-dimensional modeling and optimization of drug molecule structure according to claim 1, characterized in that, In step S2, establishing virtual rigid constraints specifically includes: In the physics engine, the first drug agent object and the second drug agent object are defined as a rigid body combination. The relative position and relative angle of the two objects with respect to their common center of mass are locked, so that the force applied to the rigid body combination in the subsequent simulation will produce an overall translational or rotational motion.

4. The method for three-dimensional modeling and optimization of drug molecule structure according to claim 1, characterized in that, In step S2, the connection logic also includes a breakage determination: In each simulation frame, the external shear stress borne by the virtual rigid constraint inside the molecular cluster object is calculated, and the corresponding equivalent deformation energy is calculated in combination with the geometric parameters of the constraint domain. When the equivalent deformation energy is greater than the binding energy threshold, the virtual rigid constraint is released, and the molecular cluster object is decomposed into independent drug intelligent agent objects.

5. The method for three-dimensional modeling and optimization of drug molecule structure according to claim 1, characterized in that, In step S3, the operational dynamics simulation specifically includes: For each simulation frame, the resultant force vector acting on each object is calculated. The resultant force vector includes at least an electrostatic attraction component pointing towards the target object and a random disturbance component simulating thermal motion. The spatial coordinates and orientation of the object are updated based on the resultant force vector using a numerical integration algorithm.

6. The method for three-dimensional modeling and optimization of drug molecule structure according to claim 1, characterized in that, In step S3, controlling the position update of the drug agent object or the molecular cluster object also includes attitude adjustment logic: When the distance between the drug agent object or the molecular cluster object and the target object is less than a preset sensing threshold, the electrostatic torque on the drug agent object or the molecular cluster object is calculated based on the electric dipole moment of the drug agent object or the potential gradient on the surface of the target object. The torque is applied to perform a rotational transformation on the drug agent object or the molecular cluster object, so that its binding sites face the target object.

7. The method for three-dimensional modeling and optimization of drug molecule structure according to claim 1, characterized in that, In step S1, constructing the initial simulation space further includes: Multiple sets of environmental interference particles are generated in the three-dimensional virtual space. The environmental interference particles have preset mass attributes and initial velocity vectors. During the simulation in step S3, when the coordinates of the environmental interference particles coincide with the coordinates of the drug agent or the molecular cluster object, the velocity vector after the collision is calculated according to the law of conservation of momentum, thereby changing the motion trajectory of the drug agent or the molecular cluster object.

8. The method for three-dimensional modeling and optimization of drug molecule structure according to claim 1, characterized in that, In step S4, the logic for generating the fitness score is as follows: The fitness score was calculated using a weighted linear combination method. The variables of the weighted linear function include at least the surface coverage, the sum of the bond energies of all virtual rigid constraints within the molecular cluster object, and the number of drug agents constituting the molecular cluster object.

9. The method for three-dimensional modeling and optimization of drug molecule structure according to claim 1, characterized in that, In step S5, the topological structure of the drug molecule is updated using an evolutionary algorithm, specifically including: The topological structure of high fitness score molecules is converted into a molecular graph; Identify non-skeleton nodes in the molecular diagram, whereby non-skeleton nodes are defined as side chain atoms or terminal atoms that do not belong to the molecular pharmacophore skeleton. A non-skeleton node is randomly selected and replaced with an atomic or functional group fragment with a higher hydrophobicity parameter to enhance the binding ability between molecules and generate the topological structure data of the progeny molecule.

Images