A method for three-dimensional modeling and optimization processing of a drug molecule structure

Abstract

The present application relates to the technical field of medicinal chemistry and bioinformatics, in particular to a kind of drug molecular structure three-dimensional modeling and optimization processing method, steps include: construct the three-dimensional virtual space including pathogen and drug intelligent agent;Set up the connection logic based on distance and energy, so that drug intelligent agent is self-assembled to form molecular cluster when meeting the condition;Run dynamics simulation, drive molecular cluster to form spatial encirclement on the surface of pathogen;Calculate the surface coverage of cluster to key site of pathogen and generate fitness score;Based on score, update the hydrophobicity modification of drug molecular topological structure using evolutionary algorithm, and iteratively generate optimal structure.The present application introduces multi-agent collaborative mechanism in the calculation model, and includes drug self-assembly and virus spatial encirclement into screening index, so that the method can find new drug structure by forming large volume cluster to physically block virus infection, effectively solve the problem of traditional method screening.

Description

Technical Field

[0001] This invention relates to the fields of medicinal chemistry and bioinformatics, and more specifically, to a method for three-dimensional modeling and optimization of drug molecular structures. Background Technology

[0002] In modern drug development, computer-aided drug design (CADD) has become a key technology for shortening development cycles and reducing costs. Current mainstream technologies mainly include structure-based drug design (SBDD), with molecular docking and molecular dynamics simulation (MD) being typical examples. The typical workflow of existing technologies involves obtaining the static crystal structure of the viral target protein, calculating the binding free energy of a single drug small molecule ligand to the target's active pocket using a scoring function, and then screening potential drugs. For example, software such as AutoDock or Schrödinger employs this type of "lock-and-key model" for screening.

[0003] However, existing technologies have the following obvious drawbacks:

[0004] First, existing methods largely ignore the swarm effect, predicting drug binding based on a one-to-one binding pattern. However, in real physiological environments with high-concentration local drug delivery, drug molecules often bind through weak intermolecular interactions (such as van der Waals forces, etc.). (Accumulation) self-assembles into dimers or clusters. These "molecular clusters" may wrap around the surface of the virus through a large steric hindrance volume, producing an inhibitory effect that monomers do not have. Existing technologies cannot simulate this dynamic assembly process, resulting in the failure to screen out drugs that have weak monomer binding force but strong group blocking ability.

[0005] Secondly, static optimization has a single path. Traditional structure-activity relationship (QSAR) analysis often relies on human experience to modify functional groups and lacks a closed-loop mechanism that can automatically reverse the evolution of molecular structures based on "spatial enclosement ability", resulting in low efficiency of structure optimization.

[0006] Therefore, a method for three-dimensional modeling and optimization of drug molecular structures is proposed. Summary of the Invention

[0007] The purpose of this invention is to provide a method for three-dimensional modeling and optimization of drug molecule structures, in order to solve the technical problems of existing drug design methods ignoring the synergistic encirclement effect of drug molecules forming clusters in real environments, and lacking an automated structural evolution mechanism based on steric hindrance.

[0008] To address the aforementioned technical problems, the present invention aims to provide a method for three-dimensional modeling and optimization of drug molecule structures, comprising the following steps:

[0009] 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;

[0010] 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.

[0011] 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.

[0012] 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.

[0013] 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.

[0014] As a further improvement to this technical solution, the preset combination conditions in step S2 specifically include:

[0015] Extract the charge distribution data of the contact surface and calculate the electrostatic potential product of the contact area;

[0016] 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.

[0017] As a further improvement to this technical solution, step S2, establishing virtual rigid constraints, specifically includes:

[0018] 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.

[0019] As a further improvement to this technical solution, in step S2, the connection logic further includes a breakage determination:

[0020] 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.

[0021] 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.

[0022] As a further improvement to this technical solution, in step S3, the operational dynamics simulation specifically includes:

[0023] 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.

[0024] The spatial coordinates and orientation of the object are updated based on the resultant force vector using a numerical integration algorithm.

[0025] As a further improvement to this technical solution, step S3, controlling the position update of the drug agent object or the molecular cluster object, further includes attitude adjustment logic:

[0026] 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.

[0027] 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.

[0028] As a further improvement to this technical solution, step S1, which involves constructing the initial simulation space, further includes:

[0029] 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.

[0030] 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.

[0031] As a further improvement to this technical solution, the logic for generating the fitness score in step S4 is as follows:

[0032] The fitness score was calculated using a weighted linear combination method.

[0033] 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.

[0034] As a further improvement to this technical solution, in step S5, the topological structure of the drug molecule is updated using an evolutionary algorithm, specifically including:

[0035] The topological structure of high fitness score molecules is converted into a molecular graph;

[0036] 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.

[0037] 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.

[0038] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0039] 1. In this method for three-dimensional modeling and optimization of drug molecular structures, a multi-agent collaborative mechanism is introduced into the computational model, and drug self-assembly and viral spatial encapsulation are incorporated into the screening indicators. This enables the method to discover novel drug structures that physically block viral infection by forming large-volume clusters, effectively solving the problem of missed screening by traditional methods.

[0040] 2. In this method for three-dimensional modeling and optimization of drug molecular structure, an automated closed loop from efficacy evaluation (coverage) to structural modification (topology update) is established. Based on the effect of synergistic encirclement, a molecular structure that is more conducive to grouping and binding is automatically evolved, thereby improving the efficiency of drug optimization.

[0041] 3. In the three-dimensional modeling and optimization method of drug molecular structure, by introducing environmental interference particles and momentum conservation calculations, the Brownian motion interference in the real body fluid environment is simulated, which makes the prediction results more robust and clinically valuable.

[0042] 4. In the three-dimensional modeling and optimization method of drug molecular structure, a weighted linear function based on coverage, bond energy and number is used as the fitness score, which not only ensures the pursuit of drug efficacy, but also takes into account the drug dosage cost and binding stability, and provides a multi-dimensional optimization direction. Attached Figure Description

[0043] Figure 1 This is a flowchart of the present invention. Detailed Implementation

[0044] 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.

[0045] Example 1

[0046] Currently, the typical process of existing technical solutions is as follows: obtain the static crystal structure of the viral target protein, calculate the binding free energy of a single drug small molecule ligand to the target active pocket using a scoring function, and screen potential drugs in this way. For example, software such as AutoDock or Schrödinger use this kind of "lock-and-key model" for screening. However, existing drug design methods have the problems of ignoring the synergistic enclosing effect after drug molecules form clusters in the real environment, and lacking an automated structural evolution mechanism based on the steric hindrance sealing effect.

[0047] In view of this, such as Figure 1 As shown, the purpose of this embodiment is to provide a method for three-dimensional modeling and optimization of drug molecule structures, using the optimization of small molecule inhibitors of the SARS-CoV-2 spike protein as an application scenario. This method for three-dimensional modeling and optimization of drug molecule structures includes:

[0048] Because constructing a virtual space that approximates the real physiological environment provides a foundation for subsequent simulations of the interaction between drug molecules and pathogens; mapping physical molecules to virtual objects facilitates dynamic control and monitoring through computer algorithms, therefore, step S1 involves constructing a three-dimensional virtual space in a computing device, mapping the target pathogen to the target object, and mapping the initial candidate drug molecules to drug intelligent agent objects with autonomous movement attributes; the specific implementation method is as follows:

[0049] First, the technical equipment selected was a workstation equipped with an Intel Xeon W-3275 processor (28 cores and 56 threads) and an NVIDIA A100 GPU (40GB of video memory). The operating system was Ubuntu 22.04 LTS, and the virtual space construction software used was Unity 3D 2022.3. Its physics engine module was used to build the spatial coordinate system (using the Cartesian coordinate system, with units in Å).

[0050] Secondly, the target pathogen was selected as the spike protein (S protein) of SARS-CoV-2. Its three-dimensional structure data was obtained from the PDB database (number 6VSB). The structure was preprocessed using PyMOL 2.5 software (removing water molecules and ligands, and optimizing the hydrogen bond network). The preprocessed structure was imported into Unity3D and mapped as a "target object". The size ratio of the target object in the virtual space was set to 1:1000 (i.e., 1 Å in the actual space corresponds to 1000 units in the virtual space). At the same time, its key sites were marked (ACE2 binding site, containing 20 amino acid residues such as Lys417, Tyr453, and His475, which was determined by the literature [S protein-ACE2 binding structure reported in Nature in 2020]).

[0051] 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°”.

[0052] 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.

[0053] 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:

[0054] 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.

[0055] The preset combination conditions include the following:

[0056] The van der Waals contact threshold is determined based on the atom type of the drug agent. The sum of the van der Waals radii of the two surface atoms of the agent is taken. For example, the sum of the van der Waals radii of C atoms and O atoms is 1.7 + 1.52 = 3.22 Å. Therefore, the van der Waals threshold for this type of contact is set to 3.2 Å (taking an approximate value considering virtual space error).

[0057] Determined by the electrostatic potential product and energy threshold, preliminary experiments (using AutoDock Vina to calculate the drug molecule binding energy under different electrostatic potential products) revealed that when the binding energy is below -5 kcal / mol, the molecular binding stability meets the requirements, corresponding to an electrostatic potential product (φ1×φ2) of -0.2 eV. 2 (Because electrostatic attraction occurs when φ1 and φ2 have opposite signs, their product is negative), therefore, it is set that "the Euclidean distance < 3.2 Å and the electrostatic potential product ≤ -0.2 eV". 2 "This represents the pre-defined combination condition;

[0058] The system is established using virtual rigid constraints and implemented using the PhysX 5.1 physics engine. When the binding conditions are met, the engine's "FixedJoint" function is called to define the two drug agents as a rigid body combination: locking the relative positions of the two agents with respect to their common centroid (the distance between the centroids is fixed to the Euclidean distance at the time of binding, such as 3.1 Å), and locking the relative angle (based on the optimal binding angle calculated by AutoDock Vina, such as the angle of 15° between the CO bond and the hydrogen bond direction at the target site). In subsequent simulations, the forces applied to this rigid body combination (such as electrostatic attraction) will act synchronously on the two agents to achieve overall translation or rotation.

[0059] Fracture determination involves calculating the external shear stress borne by the virtual rigid constraint within the molecular cluster object in each simulation frame. Based on the shear stress and contact area, the work done or energy accumulation in the constraint region is calculated. When the calculated energy value exceeds the binding energy threshold, the virtual rigid constraint is released, and the molecular cluster object is decomposed into independent drug agent objects. Specifically, each simulation frame is set to 10 ps. At the end of each frame, the external shear stress borne by the virtual rigid constraint within the molecular cluster is calculated. The external force F (unit: kcal / (mol·Å)) at the constraint is obtained in real-time using the PhysX engine, and the contact area A (calculated based on the number of atoms in the binding region, with the contact area of ​​each atom being 0.2 Å) is also calculated. 2 If three atoms are in contact, then A = 0.6 Å. 2 The current stress intensity is calculated using the shear stress formula σ=F / A; the binding energy threshold is set to -5kcal / mol (consistent with the energy threshold of the binding condition), and then the stress is converted into an energy value for comparison: the calculation formula is... (in, For 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.

[0060] 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.

[0061] 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:

[0062] The construction of a physical potential energy field specifically includes:

[0063] 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: ,in, (m is the mass of the drug agent / cluster, obtained by summing atomic masses, such as the molecular mass being 250u). for The three-dimensional coordinates of the centroid of the drug agent (or molecular cluster) at any given time (unit: Å, angstrom); To solve The three-dimensional coordinates of the centroid at time (i.e., the next time step); for The instantaneous velocity at any given moment (unit: Å / ps, angstrom / picosecond) is determined by the autonomous motion properties and the resultant force of the drug agent; The integration time step (fixed in this invention) is This value was chosen because the atomic vibration period of small molecule drugs is approximately 10 fs. The time must be less than 1 / 10 of the vibration period to ensure simulation accuracy and avoid distortion of the motion trajectory.

[0064] The initialization simulation space construction also includes: generating multiple sets of environmental interference particles in the three-dimensional virtual space, wherein the environmental interference particles have preset mass attributes and initial velocity vectors; during the simulation process 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. Specifically, H2O molecules are generated as environmental interference particles in the three-dimensional virtual space, with the number generated being 100 times the number of drug agents (simulating the solvent ratio in the physiological environment). The mass attribute of each interference particle is set to 18u (the molar mass of H2O molecules), and the initial velocity vector is based on the thermal motion velocity at 37℃ (approximately 500 m / s, random direction, calculated using Maxwell's velocity distribution); during the dynamic simulation process, when the coordinate distance between the interference particle and the drug agent is <0.5 Å (determined as a collision), the velocity after the collision is calculated according to the law of conservation of momentum. ,in , These represent the mass and pre-collision velocity of the drug-inspired agent, respectively. , These represent the mass and pre-collision velocity of the interfering particle, respectively. , The velocity after the collision is given, and energy loss is ignored (elastic collision assumption, which is consistent with the approximation of short-time collision).

[0065] 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 detected to be less than a preset sensing threshold, the electrostatic torque experienced by 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 rotation transformation on the drug agent object or the molecular cluster object so that its binding sites face the target object; specifically:

[0066] The sensing threshold is set to 5 Å (determined based on the effective binding distance between the drug agent and the target site). When the distance between the drug agent / cluster and the target object is < 5 Å, the torque relative to the potential gradient of the target object surface is calculated: first, the potential distribution on the target object surface is calculated using the Poisson-Boltzmann equation. Find the potential gradient (x, y, z). The gradient direction points in the direction of increasing electric potential. Electrostatic potential (unit: eV, electron volt) refers to the electrostatic potential at a point in a three-dimensional virtual space. With the y and z coordinates fixed, the electric potential The rate of change along the x-axis (unit: eV / Å) indicates that the potential gradually increases along the positive x-axis if the value is positive, and gradually decreases if the value is negative. The rates of change of electric potential corresponding to the y-axis and z-axis directions, respectively, have the following physical meanings and Consistent; then calculate the gradient force. (q is the charge of the drug agent), torque (r is the vector from the centroid of the drug agent to the centroid of the key site of the target object); finally, a quaternion rotation algorithm is used to perform a rotation transformation on the drug agent / cluster, with the rotation axis being the direction of the torque M and the rotation angle being the angle corresponding to the magnitude of the torque (until the binding site of the drug agent, such as the hydrogen bond donor site of the NH bond, and the hydrogen bond acceptor site facing the key site of the target object, such as the O atom);

[0067] Through the above-mentioned technical means, the resultant force control of the physical potential energy field enables the drug agent / cluster to continuously approach the target site. Random perturbation and environmental interference simulate the uncertainty of the real physiological environment, while attitude adjustment improves the alignment accuracy of the binding site. The synergistic effect of the three achieves efficient spatial encirclement of the target object surface by the drug agent / cluster, and improves the initial encirclement rate of key sites.

[0068] To quantitatively evaluate the binding effect of molecular clusters on key sites of the target, and to comprehensively consider the stability and size of the clusters to avoid optimization bias caused by a single indicator, thus providing a reliable evaluation basis for subsequent evolutionary optimization, step S4 involves calculating the surface coverage rate of the molecular cluster on the key sites of the target and generating a fitness score for the current drug molecule based on the surface coverage rate. The specific implementation method is as follows:

[0069] Surface coverage calculation:

[0070] Key site surface area determination: The surface area S0 of the key site (ACE2 binding site) of the target object was calculated using VMD 1.9.4 software. The calculation using the "Accessible Surface Area" tool yielded S0 = 200 Å. 2 ;

[0071] Cluster coverage area calculation: Using Unity 3D's ray detection function, rays are emitted into three-dimensional space from each atom on the surface of the target object's key sites (100 rays per atom, covering all directions). When a ray is blocked by atoms of a molecular cluster (ray length < 0.5 Å, considered coverage), the corresponding surface area element is recorded (0.01 Å per ray). 2 ), summing them up gives the area S1 of the critical region covered by the cluster;

[0072] Surface coverage calculation: Coverage C = (S1 / S0) × 100%, for example, S1 = 130 Å 2 At that time, C=65%;

[0073] Fitness score generation: The fitness score F is calculated using a weighted linear combination formula, the formula is as follows: Where a, b, and c are weighting coefficients, determined using the Analytic Hierarchy Process (AHP): Five experts in drug design were invited to score the importance of "binding effect (C)," "cluster stability (E)," and "molecular size (N)," constructing a judgment matrix and calculating the weights. Ultimately, a=0.5 (binding effect is the most important), b=0.3 (stability is secondary), and c=0.2 (molecular size needs to be controlled to avoid excessive size leading to decreased cell membrane permeability). E is the sum of the bond energies of all virtual rigid constraints within the molecular cluster: the bond energy of each virtual rigid constraint is set to -5 kcal / mol based on the binding energy (energy at which binding is stable). The cluster contains... Each constraint ,For example When the value is 3, E = -15 kcal / mol. To facilitate calculation, E is normalized. ,in =-50kcal / mol (maximum number of constraints: 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 .

[0074] 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.

[0075] 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:

[0076] The evolutionary algorithm parameter settings include the following:

[0077] 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;

[0078] 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:

[0079] The topological structure of molecules with high fitness scores (selecting 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 atoms (single bonds, double bonds, and triple bonds are labeled as 1, 2, and 3, respectively), thus completing the molecular graph conversion.

[0080] The “molecular pharmacophore skeleton” is defined as the core structure in the molecule that maintains the pharmacological effect (such as the benzene ring and amide bond “-CONH-” in this embodiment). 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) in the molecule, thereby completing the identification of non-skeleton nodes.

[0081] A non-skeleton node (such as the H atom in a methyl group) is randomly selected and replaced with an atom or functional group fragment with a higher hydrophobicity parameter. The hydrophobicity parameter is represented by the logP value (octanol-water partition coefficient; the higher the logP value, the stronger the hydrophobicity). The original H atom has a logP value of -0.5 and is replaced with a Cl atom (logP value +0.3) or trifluoromethyl (-CF3, logP value +1.0). During the replacement process, the rationality of chemical bonds must be ensured (e.g., if the H atom is a single bond, the replaced Cl atom should also be a single bond to avoid unreasonable structures such as pentavalent C atoms). After the replacement, the 3D structure of the molecule is regenerated using RDKit, and the new charge, van der Waals radius, and other parameters are calculated as a new generation of drug intelligence object, thus completing the replacement of the non-skeleton node.

[0082] Convergence criteria include the following:

[0083] Two convergence conditions are set, and iteration stops when either condition is met:

[0084] The average fitness score of the population changed by less than 1% over 10 consecutive generations (i.e.) , (The average score of the nth generation) indicates that the scores tend to stabilize;

[0085] In a certain generation, there are molecules with a fitness score ≥0.9 (preliminary experiments show that molecules with a score ≥0.9 meet the requirements for practical applications).

[0086] When the convergence condition is met, the three-dimensional structure data (including atomic coordinates, chemical bond type, and charge distribution) of the molecule with the highest fitness score in the population is output. This data can be exported in PDB format for subsequent experimental synthesis and activity verification, thus achieving the optimal molecule output.

[0087] Through the aforementioned techniques, the evolutionary algorithm, via a cycle of selection, updating, and filtering, can rapidly improve the fitness score of molecules, thereby enhancing the fitness score, surface coverage, and hydrophobicity parameter of the optimal molecule. (Enhancing cell membrane permeability) has led to the optimization of drug molecule structure.

[0088] Example 2

[0089] This embodiment is a variation of Example 1, using the optimization of β-lactam drugs targeting the penicillin-binding protein of Escherichia coli as the application scenario. The core steps are the same as in Example 1, only the key parameters are adjusted according to the characteristics of the target pathogen and the drug molecule to verify the universality of this method:

[0090] Step S1: The target pathogen is Escherichia coli penicillin-binding protein (PBP2a, PDB number 1VQQ), and the key site is the region surrounding 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 velocity 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, so the velocity is appropriately reduced);

[0091] Step S2: The van der Waals contact threshold is adjusted according to the atom type of penicillin G. For example, the threshold for S atoms (van der Waals radius 1.85 Å) and O atoms 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).

[0092] Step S3: In the electrostatic attraction component, the total charge at the key site of PBP2a is -1e, and the standard deviation of the random perturbation component is set to 0.08 kcal / (mol·Å) (larger molecular mass reduces thermal perturbation); in addition to H2O, a small amount of Mg is added to the environmental interference particles. 2 + Ions (simulating the intracellular environment of bacteria), with a mass attribute set to 24u and an initial velocity set to 400m / s;

[0093] Step S4: Surface area of ​​key sites S0 = 150 Å 2 The fitness score weights were adjusted to a=0.55 (the PBP2a key site is small, so coverage is more important), b=0.25, and c=0.2;

[0094] Step S5: The non-skeleton node is the terminal atom of the phenylacetyl side chain of penicillin G, which is replaced with a naphthylacetyl group (-COCH2C10H7, with a higher logP value). After 80 iterations, the evolutionary algorithm meets the convergence condition.

[0095] This embodiment successfully applied the method to the optimization of antibacterial drugs by adjusting parameters, demonstrating the method's versatility; the optimized drug molecules showed a significantly improved binding ability to PBP2a, providing a new technical path for the development of drugs against drug-resistant bacteria.

[0096] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely preferred examples and are not intended to limit the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention 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. 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. 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.

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 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.

6. 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.

7. 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