Method and system for predicting GPCR activation ability of agonist molecules

Abstract

A method and system for predicting the GPCR activation ability of agonist molecules, wherein the method comprises: using global molecular docking to blindly guess complex structures formed by a ligand structure binding to an active-state receptor structure and an inactive-state receptor structure, respectively; on the basis of an enhanced sampling algorithm, extracting an initial path capable of activating an inactive-state complex structure into an active-state complex structure; using an automatic path optimization algorithm to find a minimum free energy path closest to the initial path; and using umbrella sampling to calculate a free energy distribution curve along the minimum free energy path, and determining an energy barrier height and a free energy difference before and after activation, so as to determine the GPCR activation ability of the ligand structure.

Description

Method and system for predicting the GPCR activation ability of agonist molecules Technical Field

[0001] The present invention relates to the field of computational simulation research of biomolecular systems, and in particular to a method and system for predicting the ability of agonist molecules to activate GPCRs. Background Art

[0002] G protein-coupled receptors (GPCRs) are the largest membrane protein superfamily encoded in the human genome (comprising 800 members). As the most widely distributed intracellular receptors, GPCRs can be activated by a variety of ligands, including odorants, pheromones, hormones, neurotransmitters, and chemokines. Through complex signal transduction and regulatory mechanisms, they significantly influence nearly every aspect of human physiology and metabolism, playing a crucial role in numerous diseases, including cardiovascular disease, neurological disorders, inflammation, metabolic diseases, and cancer. Over 30% of approved drugs worldwide target proteins in GPCR-related signaling pathways.

[0003] GPCRs are composed of seven transmembrane α-helices (TM1-7), three extracellular loops (ECL), and three intracellular loops (ICL), as shown in Figure 1. When the extracellular end (typically composed of the extracellular loops and the extracellular end of the transmembrane region) binds to a ligand, the GPCR undergoes a conformational change, entering an activated state. Consequently, the binding strength of the intracellular end to the downstream G protein changes, inducing the G protein to exhibit properties of a guanine nucleotide exchange factor. This property facilitates the exchange of guanosine diphosphate (GDP) originally bound to the G protein for guanosine triphosphate (GTP), further promoting the separation of the G protein's α subunit from its β and γ subunits, each of which activates further downstream signaling pathways. These signaling pathways include adenylyl cyclase and the potassium ion channel GIRK, as shown in Figure 2.

[0004] Research on GPCRs remains a hot topic in academia and pharmaceutical companies. After nearly 20 years of dedicated effort, structural biologists have broadly delineated several key states during GPCR activation and successfully resolved high-resolution structures of multiple GPCRs in both activated and inactivated states. This plethora of known structures has significantly advanced the design and development of GPCR-targeted drugs. By comparing the resolved inactive and activated structures of class A GPCRs (such as rhodopsin and β2AR), key differences between GPCRs before and after activation have been identified. Despite the diverse nature of GPCR ligands, these changes are generally considered common hallmarks of GPCR activation.

[0005] There are a lot of detailed differences in the process of different agonists activating the corresponding GPCRs, and it is these detailed differences that determine the specific (selective) activation of GPCRs by GPCR agonists. For the design of drugs targeting GPCRs, ligand specificity is a core issue that cannot be avoided. Ligand specificity is directly related to whether it can activate a certain target GPCR while failing to activate other GPCR subtypes, and is closely related to issues such as its activity and toxicity. Ligand specificity is often rooted in the structural and kinetic details of its interaction with GPCRs. In fact, subtle differences between different ligands can lead to their activation of different subtypes in the same GPCR family.

[0006] To achieve rational design of GPCR agonists, researchers must be able to predict the specific activation ability of ligands for a given GPCR and possess the means to gain a detailed understanding of the mechanisms of ligand-GPCR binding and activation. Only with a mechanistic understanding of the activation process can we obtain the necessary information to optimize and improve the agonist's molecular structure.

[0007] However, current technologies struggle to meet this need. Functional experiments can only confirm activation. Structural biology, while providing atomic-level detail, reveals only a few states of the process and cannot pinpoint the ease or mechanism of activation. Identifying key protein residues that interact with ligands by mutating sites around the protein pocket one by one is effective but prohibitively expensive and cumbersome.

[0008] At the same time, the main technologies that can predict contemporary computer-aided drug design (CADD), such as molecular docking and free energy perturbation, all focus on the calculation of ligand binding strength and do not involve the activation process of the receptor after the ligand binds to the target receptor.

[0009] Therefore, developing algorithmic tools that can efficiently reveal the activation mechanisms and potency of candidate agonists on GPCRs is not only an intrinsic requirement for studying ligand activation specificity but also a necessary step for the in-depth development of CADD. Currently, there is no algorithmic framework that can quickly and accurately predict the activation potency of any ligand on any given GPCR, making the rational design of highly specific GPCR agonists difficult.

[0010] Summary of the Invention

[0011] In view of this, the present invention proposes a method and system for predicting the GPCR activation ability of agonist molecules. The specific scheme is as follows:

[0012] A method for predicting the GPCR activation ability of an agonist molecule, comprising:

[0013] Obtain the initial ligand structure from the agonist, obtain the inactive and activated receptor structures of a given GPCR, and determine whether the structure of the complex formed by the ligand and GPCR is known;

[0014] If the structure is unknown, the global molecular docking technology is used to blindly guess the activated state complex structure and the inactivated state complex structure formed by the ligand structure binding to the activated state receptor structure and the inactivated state receptor structure, and local structure optimization is performed;

[0015] Based on the preset enhanced sampling algorithm, the path that can activate the inactive state complex structure to the active state complex structure is extracted to obtain the initial path of the activation process;

[0016] Using a preset automatic path optimization algorithm to find the minimum free energy path closest to the initial path;

[0017] Umbrella sampling is used to calculate the free energy distribution curve along the minimum free energy path. The energy barrier height and the free energy difference before and after activation are determined from the free energy distribution curve to judge the activation ability of the ligand structure on the GPCR.

[0018] In one embodiment, partial structural information of the minimum free energy path is used to optimize the ligand structure to increase or decrease the ability of the agonist molecule to activate a given GPCR.

[0019] In a specific embodiment, the preset automatic path optimization algorithm includes: an automatic path search method based on the traveling salesman problem, a finite temperature chord method, a fast tomography method, and a transition path sampling method.

[0020] In a specific embodiment, the preset enhanced sampling algorithm includes a coarse-grained model, a high-temperature simulation algorithm, a stretching dynamics algorithm, a targeted dynamics algorithm, and a data-driven path algorithm.

[0021] In a specific embodiment, the automated path search method based on the traveling salesman problem specifically includes:

[0022] The RMSD between any two structures is defined by all heavy atoms of the complex structure and then PCV is calculated. PCV projects the high-dimensional space onto a temporary coordinate system constructed based on the current initial path.

[0023] Rapidly find the minimum free energy path fragment near the initial path by sampling in orthogonal directions;

[0024] Set up a virtual node with a distance of zero from the known path node, and perform a traveling salesman solution on the path nodes on all minimum free energy path segments including the virtual node. After removing the virtual node, complete the sorting and obtain a new initial path.

[0025] Correct the spatial resolution of the new initial path and complete a round of iteration;

[0026] When the change amplitude of the obtained initial path is lower than the preset amplitude, the iteration is stopped and the initial path is used as the minimum free energy path.

[0027] In a specific embodiment, for the molecular docking structures of the activated state complex and the inactivated state complex obtained by blind guessing, a short-term molecular dynamics simulation is used to perform local structure optimization to obtain an accurate two-state complex structure model.

[0028] In one embodiment, external forces are applied to all heavy atoms of the inactive state complex structure to pull the inactive state complex structure to the active state complex structure model, thereby generating an initial path for the activation process.

[0029] In a specific embodiment, the inactive state receptor structure and the activated state receptor structure of a given GPCR are obtained by known structure, homology modeling or AlphaFold.

[0030] In a specific embodiment, if the current ligand structure does not meet expectations in its ability to activate a GPCR, the ligand structure is optimized and modified based on the completed structure of the transition state and intermediate metastable state of the minimum free energy path to obtain a new ligand structure, and the activation ability of the new ligand structure on the GPCR is re-predicted.

[0031] A prediction system for the GPCR activation ability of agonist molecules, including the following modules:

[0032] A preparation module is used to obtain the initial ligand structure from the agonist, obtain the inactive receptor structure and the activated receptor structure of a given GPCR, and determine whether the structure of the complex between the ligand structure and the GPCR is known;

[0033] The blind guessing module is used to use global molecular docking technology to blindly guess the activated state complex structure and the inactivated state complex structure formed by the ligand structure binding to the activated state receptor structure and the inactivated state receptor structure respectively when the complex structure is unknown, and perform local structure optimization;

[0034] An initial path module is used to extract a path that can activate the inactive state complex structure to the active state complex structure based on a preset enhanced sampling algorithm, thereby obtaining an initial path of the activation process;

[0035] A free energy path module, configured to use a preset automatic path optimization algorithm to find a minimum free energy path closest to the initial path;

[0036] The activation ability prediction module is used to calculate the free energy distribution curve along the minimum free energy path using umbrella sampling, and to determine the energy barrier height and the free energy difference before and after activation through the free energy distribution curve, so as to judge the activation ability of the ligand structure on GPCR.

[0037] Beneficial effects: The present invention provides a method and system for predicting the ability of agonist molecules to activate GPCRs. Under the premise of assuming that the ligand structure is unknown, no information is input in the entire prediction process. The prediction is highly versatile and efficiently reveals the activation mechanism and activation ability of candidate agonists on GPCRs. The complex structure of the ligand bound to the activated and inactivated states of the GPCR is blindly guessed. All heavy atoms of the complex structure are used as input in the generation of the initial path, path optimization, and free energy curve calculation, avoiding interference of external selection on the prediction results. The resulting MFEP consists of a series of complete structures of the complex, which is sufficient to provide detailed mechanism information, provide guiding information for the optimization and improvement of agonist molecules, and complete its rational design iteration. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] FIG1 is a schematic diagram of the domain structure of GPCR;

[0039] Figure 2 is a schematic diagram of the GPCR signal transduction process;

[0040] FIG3 is a schematic flow chart of the prediction method of the present invention;

[0041] FIG4 is a schematic diagram of the free energy surface results along the minimum free energy path of the present invention;

[0042] FIG5 is a comparative schematic diagram of three key states along the path of the present invention;

[0043] FIG6 is a schematic diagram of the prediction system module of the present invention.

[0044] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments. It should be understood that the following drawings only illustrate certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other relevant drawings can be obtained based on these drawings without paying any creative work.

[0045] Figure numerals: 1-preparation module; 2-blind guessing module; 3-initial path module; 4-free energy path module; 5-activation ability prediction module; 6-ligand structure modification module. DETAILED DESCRIPTION

[0046] Hereinafter, various embodiments of the present disclosure will be described more fully. The present disclosure may have various embodiments, and modifications and variations may be made therein. However, it should be understood that there is no intention to limit the various embodiments of the present disclosure to the specific embodiments disclosed herein, but rather that the present disclosure should be construed to encompass all modifications, equivalents, and / or alternatives falling within the spirit and scope of the various embodiments of the present disclosure.

[0047] All-atom molecular dynamics (MD) simulations are a computational method based on Newtonian dynamics to study the microscopic properties of biomolecules. Their goal is to simulate the time-dependent structural fluctuations of biomolecules under thermal perturbations. After years of development, MD simulations have become a standard tool for studying complex biomolecules, and their application in studying the mechanisms of GPCR activation has also become widespread. However, predicting the GPCR activation ability of candidate agonists using conventional MD simulations remains challenging. MD simulations on standard high-performance computing devices (such as GPU clusters and CPU resources at supercomputing centers) typically only last nanoseconds to microseconds, while GPCR activation can require overcoming a high energy barrier, occurring on timescales of hundreds of microseconds or even longer. Furthermore, when predicting the activation ability of candidate agonists, the height of this barrier is unknown. Rational agonist design requires the assumption that both the agonist's molecular structure and its binding site on the GPCR are unknown. Furthermore, the resulting MFEP must consist of a series of complete complex structures. Otherwise, the MFEP is insufficient to provide detailed mechanistic information, guide agonist optimization, and enable rational design iterations. Because the ligand structure is unknown, the activated-state structure of the complex is also unknown. While the inactive structure of the complex can be inferred through molecular docking, the activated-state structure is typically generated by applying external forces to the inactive-state model using a known activated-state structure of a GPCR with high sequence similarity as a template. During this process, the specific atoms to which the forces are applied are unknown.

[0048] When running MD on ordinary high-performance computing equipment (such as GPU clusters and supercomputing center CPU resources), the total length of the simulation can usually only reach the nanosecond to microsecond level, while the process of GPCR activation may need to cross a higher energy barrier, which occurs on a time scale of hundreds of microseconds or even longer, and when predicting the activation ability of candidate agonist molecules, the height of the energy barrier is unknown in advance. The terms used in the various embodiments disclosed in the present invention are only used to describe the purpose of specific embodiments and are not intended to limit the various embodiments disclosed in the present invention. As used herein, the singular form is intended to also include the plural form, unless the context clearly indicates otherwise. Unless otherwise specified, all terms used herein (including technical terms and scientific terms) have the same meaning as that generally understood by ordinary technicians in the field of the various embodiments disclosed in the present invention. The terms (such as those defined in generally used dictionaries) will be interpreted as having the same meaning as the contextual meaning in the relevant technical field and will not be interpreted as having an idealized meaning or an overly formal meaning, unless clearly defined in the various embodiments disclosed in the present invention.

[0049] Example 1

[0050] Example 1 of the present invention discloses a method for predicting the GPCR activation ability of agonist molecules. No information is input throughout the prediction process, resulting in a highly versatile prediction method that efficiently reveals the activation mechanism and activation ability of candidate agonists on GPCRs. A flowchart of the prediction method is shown in Figure 3 of the specification. The specific scheme is as follows:

[0051] A method for predicting the GPCR activation ability of an agonist molecule comprises the following steps:

[0052] 101. Obtain the initial ligand structure from the agonist, obtain the inactive and activated receptor structures of a given GPCR, and determine whether the structure of the complex formed by the ligand and the GPCR is known;

[0053] If the structure is unknown, use global molecular docking technology to blindly guess the activated state complex structure and the inactivated state complex structure formed by the ligand structure binding to the activated state receptor structure and the inactivated state receptor structure, and perform local structure optimization;

[0054] 103. Extracting the path that can activate the inactive complex structure into the active complex structure based on the preset enhanced sampling algorithm, and obtaining the initial path of the activation process;

[0055] 104. Use the preset automatic path optimization algorithm to find the minimum free energy path closest to the initial path;

[0056] 105. Use umbrella sampling to calculate the free energy distribution curve along the minimum free energy path. The energy barrier height and the free energy difference before and after activation are determined through the free energy distribution curve to judge the activation ability of the ligand structure on GPCR.

[0057] In a preferred embodiment, after determining the activation ability of the ligand structure on a given GPCR, the following steps may be further included:

[0058] 106. Optimize the ligand structure using partial structural information of the minimum free energy path to increase or decrease the ability of the agonist molecule to activate a given GPCR.

[0059] In order to improve the versatility of the prediction method for ligand structures in this embodiment, no information is input in all steps, and global blind molecular docking is performed. In addition, all heavy atoms of the complex structure are used as input in the generation of the initial path, path optimization, and free energy curve calculation to minimize the influence of external factors.

[0060] This example predicts the ability of an agonist to activate a given GPCR, using the agonist as the ligand and the GPCR as the receptor. The receptor structures in both the activated and inactivated states can be obtained using known structures, homology modeling, or AlphaFold.

[0061] Among them, the structure of the agonist molecule can be known or unknown (such as an activator is given at random). When the ligand structure is unknown (the structure of the activator molecule is unknown), the activated state structure of the complex formed by the ligand and the receptor is also unknown. Although the inactivated state structure of the complex can be guessed through molecular docking, the activated state structure must usually be based on a known activated state structure of a GPCR with high sequence similarity as a template, starting from the inactivated state structure model and generating it by applying external force through tMD. In this process, it is unknown which atoms are applied by the external force, which makes it impossible to make subsequent predictions.

[0062] Therefore, in addition to obtaining the structures of the receptor in its activated and inactivated states in step 101, it is also necessary to determine whether the structure of the complex formed by the ligand and receptor is known. For example, structural biology research methods such as cryo-electron microscopy and crystal structure analysis can be used to obtain the complex formed by the ligand and receptor. If the complex structure is known, there is no need to use global molecular docking technology to make blind guesses, and the structures of the activated and inactivated complexes can be directly obtained.

[0063] Cryoelectron microscopy (cryo-EM) is a method that uses electron microscopy to observe frozen samples. The advantage of cryo-EM is that it can directly observe samples in liquid or semi-liquid states, such as biomacromolecules and cell membranes. Compared to other structural biology methods, cryo-EM can more easily obtain high-resolution structures.

[0064] X-ray crystallography is a method of determining the structure of molecules by irradiating a crystalline sample with X-rays. This method requires that the molecules be arranged in a crystal so that the X-rays can pass through the crystal and produce an analyzable diffraction pattern. The advantage of X-ray crystallography is that it can obtain high-resolution structures, which can reach Even higher.

[0065] In practical applications, the structure of the complex formed by the binding of a ligand to a receptor is more common. This embodiment addresses this situation by using global molecular docking technology to perform blind guesses without adding any input involving external influences. Molecular docking technology uses the characteristics of the receptor and the interaction between the receptor and the ligand to perform theoretical simulations, studying intermolecular interactions (such as ligand-receptor interactions) and predicting their binding mode and affinity.

[0066] The structures of the activated and inactivated complexes are formed by blindly guessing the ligand structures bound to the activated and inactivated receptor structures, respectively. The best of these blindly guessed structures are then selected for local structure optimization using short-time molecular dynamics simulations to improve simulation accuracy and reliability. Short-time molecular dynamics simulations are a computational biophysics method used to study the structural and dynamic properties of biomacromolecules. By simulating the motion of biomacromolecules in a computer, they provide information about their structure and dynamics over short timescales. In local structure optimization, short timesteps typically refer to using smaller timesteps to improve simulation accuracy and reliability. This helps avoid errors caused by vibration or jitter associated with longer timesteps, resulting in more accurate optimization. Therefore, it is recommended to use shorter timesteps when performing local structure optimization in molecular dynamics simulations to improve simulation accuracy and reliability.

[0067] The solution of this embodiment is to preset an enhanced sampling algorithm to generate an initial path. The initial path is the path from the inactive complex structure to the activated complex structure. Generating the initial path can better understand the structure and reaction mechanism of the molecule, and can also improve the efficiency and accuracy of our subsequent simulations. The enhanced sampling algorithm is an algorithm used to improve the efficiency of sampling in molecular dynamics simulations. The preset enhanced sampling algorithms include coarse-grained models, high-temperature simulation algorithms, Steered MD algorithms for stretching dynamics, targeted MD algorithms for targeted dynamics, and data-driven path algorithms.

[0068] Coarse-grained model: A method to simplify the structure of biological macromolecules by integrating atoms or small molecules in the original molecule into larger components, thereby reducing the amount of calculation and accelerating sampling efficiency.

[0069] High-temperature simulation: This algorithm accelerates molecular motion by increasing the temperature, thereby increasing sampling efficiency. Implementing high-temperature simulation requires setting parameters such as the simulation temperature, time step, and force constant.

[0070] Steered MD: This algorithm applies an external force to a specific atom in a molecule to control its motion, thereby accelerating sampling efficiency. Implementing Steered MD requires setting parameters such as the direction of the force, the force constant, and the time step.

[0071] Targeted MD: Using a specific structure as a target, starting from another structure, external forces are applied to certain atoms on the molecule to drive them toward the target structure, thereby accelerating sampling. Implementing Steered MD requires setting parameters such as the atoms in the molecule to apply the force, the force constant, and the time step.

[0072] Data-driven pathway algorithms: These algorithms use machine learning algorithms to predict possible reaction pathways based on known experimental data or computational results. Implementing data-driven pathway algorithms requires preparing training datasets, selecting appropriate machine learning models, and adjusting model hyperparameters.

[0073] After obtaining the initial path, it is necessary to find the minimum free energy path. The minimum free energy path of this embodiment is based on a high-dimensional space and consists of a series of complete complex structures, rather than the low-dimensional structure in existing path algorithms. Existing algorithms rely on prior guesses about the activation process mechanism, and the resulting path is defined in a pre-selected low-dimensional coordinate space. However, the complex of GPCR and ligand is precisely more complex, and the probability of selecting errors is high, resulting in its overall low efficiency. The MFEP defined in a low-dimensional space cannot provide inspiration for the next step of improving the ligand structure.

[0074] Preferably, the RMSD (Root Mean Square Distance) between any two structures is defined using all the heavy atoms in the complex. This is then used to calculate the Path-Collective-Variable (PCV). Starting from this initial path, an automated path optimization algorithm is then used to rapidly find the nearest minimum free energy path (MFEP). This is equivalent to defining the PCV using all the heavy atoms in the complex. Using this optimization algorithm, a high-dimensional minimum free energy path for the activation process can be obtained in one go, significantly reducing the computational effort and providing a foundation for subsequent multiple rounds of iterative calculations.

[0075] RMSD (Root Mean Square Deviation) is a common method to measure the deviation of atoms or atomic groups between two structures.

[0076] Path-Collective-Variable (PCV) is defined as follows:

[0077] where d x,i It is the RMSD distance between any conformation in the original high-dimensional space and the i-th node on the high-dimensional path (containing N high-dimensional node concepts), and λ is a constant determined by the average RMSD between adjacent node conformations on the path. Through PCV, the original high-dimensional space can be projected onto a temporary coordinate system based on the current path of this round.

[0078] Automated path optimization algorithms include automated path search methods based on the traveling salesman problem, string methods, fast projection methods, and transition state path sampling. These algorithms differ in their implementation and application scenarios, but all aim to find a minimum free energy path through optimization or sampling—a path connecting the starting state and the target state that minimizes energy or free energy.

[0079] The Traveling-salesman based Automated Path Searching method gradually approaches the minimum free energy path through multiple iterations. The Finite Temperature String Method represents the structural change path as a sequence of atomic coordinates in a low-dimensional space, and optimizes it through the minimum free energy principle to find the minimum free energy path in the low-dimensional space. The Fast Tomographic Method gradually approaches the true path by iteratively updating the initial path. It is suitable for studying the interaction between small molecules and large molecules. The Transition Path Sampling method generates a series of transition state trajectories through the Monte Carlo method, and then constructs the minimum free energy path. It is suitable for larger molecular systems.

[0080] Finally, the RMSD and PCV were defined again using all heavy atoms in the complex. Umbrella sampling was used to calculate the free energy profile along the MFEP (PCV-s component). The free energy difference before and after activation and the activation barrier height were used to determine the ligand's ability to activate the receptor. The free energy profile derived from umbrella sampling can determine the barrier height and the free energy difference before and after activation, which are key information for predicting the ability of a ligand to activate a GPCR. This information can be cross-correlated with the results of GPCR functional experiments (bioluminescence resonance energy transfer (BRET) assays, which measure the dissociation of downstream G proteins from the GPCR).

[0081] Among them, MetaDynamics can be used to calculate the free energy curve. MetaDynamics is a method for calculating free energy curves. It uses a strategy called "evolutionary dynamics" to calculate free energy. It is often used to calculate the free energy of proteins or other biological molecules to understand their behavior in biological environments. In MetaDynamics, the system is subjected to a small force, and then the free energy is calculated by calculating the evolutionary behavior of the system. This method usually involves simulations at multiple time scales, from microscopic to macroscopic simulations, and systems with multiple degrees of freedom. MetaDynamics can be combined with other methods, such as molecular docking, molecular dynamics simulation, Monte Carlo simulation, etc., to obtain more accurate free energy curves. This method can provide information on fields such as drug design, protein engineering, biophysics, etc., and help understand the structure and function of biomolecules.

[0082] The high-dimensional minimum free energy pathway consists of a series of complete complex structures, where the completed structures of the transition state and intermediate metastable states can provide insights into modifying the ligand structure to increase or decrease its ability to activate a given GPCR. If the current ligand structure does not meet expectations for GPCR activation, the ligand structure is optimized and modified based on the completed structures of the transition state and intermediate metastable states along the minimum free energy pathway to obtain a new ligand structure, and its GPCR activation ability is re-predicted. Due to the low computational complexity required in this embodiment, the integration of consumer-grade GPU servers makes it possible to perform multiple rounds of iterative calculations and rational design of ligand molecules, which is of groundbreaking significance.

[0083] This example provides a method for predicting the GPCR activation ability of agonist molecules. Under the assumption that the ligand structure is unknown, no information is input throughout the prediction process. The prediction is highly versatile and can efficiently reveal the activation mechanism and activation ability of candidate agonists on GPCRs. The structures of the complexes of the ligand bound to the activated and inactivated GPCRs are blindly guessed. All heavy atoms of the complex structure are used as input in initial path generation, path optimization, and free energy curve calculation, avoiding interference from external selection on the prediction results. The resulting MFEP, consisting of a series of complete complex structures, is sufficient to provide detailed mechanistic information, provide guidance for the optimization and improvement of agonist molecules, and complete their rational design iterations.

[0084] Example 2

[0085] Example 2 of the present invention discloses a preferred solution of Example 1. The solution is as follows:

[0086] A method for predicting the GPCR activation ability of an agonist molecule comprises the following steps:

[0087] 201. Use AlphaFold or homology modeling to obtain the receptor structures of the inactive and activated states of a given GPCR. Then use global molecular docking to blindly guess the structural models of the complexes of the ligand bound to the receptor in the two states, and select the best one from them to perform local structure optimization using short-term molecular dynamics simulations.

[0088] 202. Apply external forces to all heavy atoms of the complex to pull the inactive complex structure to the active complex structure model, thus generating the initial path of the activation process;

[0089] 203. The RMSD (Root Mean Square Distance) between any two structures is defined using all the heavy atoms in the complex, and the path-collective-variable (PCV) is calculated. TAPS, an automated path optimization algorithm based on the traveling salesman problem, is used to quickly find the nearest minimum free energy path (MFEP) starting from the initial path.

[0090] 204. Again, all heavy atoms of the complex are used to define RMSD and PCV. Umbrella sampling is used to calculate the free energy distribution curve along the MFEP (PCV-s component). The activation ability of the ligand on the receptor is determined by the free energy difference before and after activation and the activation energy barrier height.

[0091] 205. Use the structural information of transition states and intermediate states on MFEP to optimize the chemical structure of the ligand. Return to 201 to predict the activation ability of the new ligand structure.

[0092] The Traveling-salesman-based Automated Path Searching (TAPS) algorithm, based on the traveling salesman problem, rapidly finds minimum free energy path segments near the initial path through orthogonal sampling. It then uses a traveling salesman solver to reorder the resulting nodes and adjust the spatial resolution of the new path, repeating this process multiple times until the path no longer changes significantly. In short peptide systems, TAPS's search efficiency is at least eight times that of other methods. The traveling salesman rearrangement process involves creating a virtual node with zero distance from a known node and performing a traveling salesman solver on all nodes, including the virtual node. The sorting is then completed after the virtual node is removed.

[0093] The solution of this embodiment adopts non-local sampling to quickly search for MFEP fragments in the vicinity of the current path in all directions orthogonal to the current path, thereby speeding up the optimization process. Under the premise that the String Method guesses the static coordinate space in advance, the computational efficiency of TAPS in the short peptide system is 8 times that of the String Method. At the same time, since the new path nodes obtained by non-local sampling in the orthogonal direction may not follow the arrangement order in the previous round, this embodiment introduces the traveling salesman algorithm to reorder the new path nodes. Since a static coordinate system is not introduced in TAPS, in each round of iteration, the calculation method of the distance (Root-Mean-Square-Distance, RMSD) between the two conformations is defined (which atoms are selected to align the structure, and which atoms are used to calculate the RMSD), and then the path variable (Path Collective Variable, PCV) is calculated. It is defined as follows:

[0094] where d x,i It is the RMSD distance between any conformation in the original high-dimensional space and the i-th node on the high-dimensional path (containing N high-dimensional node concepts), and λ is a constant determined by the average RMSD between adjacent node conformations on the path.

[0095] Through PCV, the original high-dimensional space can be projected onto a temporary coordinate system based on the current path, thus bypassing the difficulty of selecting a static coordinate system. This significantly improves the automation and overall efficiency of path optimization, successfully analyzing the functional conformational change mechanisms of multiple complex systems at a low cost. Unlike other methods that use prior assumptions about the mechanism of the process under investigation as input, this embodiment directly obtains high-dimensional minimum free energy paths based on path variables without pre-defining a physically meaningful coordinate space. This minimizes the separation of the sampling process and mechanistic interpretation, significantly improving the algorithm's automation. Using the solution of this embodiment, a minimum free energy path for sphingosine-1-phosphate (S1P) activation of the receptor S1PR1 was obtained after 12 days of calculations on a server with 8 consumer-grade GPUs (NVIDIA RTX2080Ti). The minimum free energy path for S1PR activation by the ligand S1P is shown in Figures 4 and 5. Figure 4 shows the free energy surface (umbrella sampling) along the obtained minimum free energy path. Figure 5 shows a comparison of the structures of the three key states along the pathway (inactive state, transition state, and pre-activation intermediate state) in the binding pocket, the middle of the receptor, and the intracellular G protein binding region. The key residues that interact with S1P and reflect the main changes in the receptor activation process are marked in the inactive state.

[0096] Example 3

[0097] This example provides a system for predicting the ability of agonist molecules to activate GPCRs. Based on Example 1, the method of Example 1 is systematized. The specific structure is shown in Figure 6 of the specification. The specific scheme is as follows:

[0098] Preparation module 1 is used to obtain the initial ligand structure from the agonist, obtain the inactive state receptor structure and the activated state receptor structure of a given GPCR, and determine whether the structure of the complex between the ligand structure and the GPCR is known;

[0099] Blind guessing module 2 is used to blindly guess the activated state complex structure and the inactivated state complex structure formed by the ligand structure binding to the activated state receptor structure and the inactivated state receptor structure respectively, when the structure is unknown, and perform local structure optimization;

[0100] An initial path module 3 is used to extract a path that can activate the inactive state complex structure to the active state complex structure based on a preset enhanced sampling algorithm, thereby obtaining an initial path of the activation process;

[0101] Free energy path module 4, used to use a preset automatic path optimization algorithm to find the minimum free energy path closest to the initial path;

[0102] The activation ability prediction module 5 is used to calculate the free energy distribution curve along the minimum free energy path using umbrella sampling, and determine the energy barrier height and the free energy difference before and after activation through the free energy distribution curve, so as to judge the activation ability of the ligand structure on GPCR.

[0103] The ligand structure modification module 6 is used to optimize the ligand structure using partial structural information of the minimum free energy path to increase or decrease the ability of the agonist molecule to activate a given GPCR.

[0104] This example provides a prediction system for the GPCR activation ability of agonist molecules, which systematizes the method of Example 1 to make it more practical.

[0105] The present invention provides a method and system for predicting the ability of agonist molecules to activate GPCRs. Under the assumption that the ligand structure is unknown, no information is input during the entire prediction process. The prediction is highly versatile and effectively reveals the activation mechanism and activation ability of candidate agonists on GPCRs. The structure of the complexes of the ligand bound to the activated and inactivated states of the GPCR is blindly guessed. All heavy atoms of the complex structure are used as input during initial path generation, path optimization, and free energy curve calculation, avoiding interference with the prediction results due to external selection. The resulting MFEP, consisting of a series of complete complex structures, is sufficient to provide detailed mechanistic information, provide guidance for the optimization and improvement of agonist molecules, and complete their rational design iterations.

[0106] Those skilled in the art will understand that the accompanying drawings are only schematic diagrams of a preferred implementation scenario, and the modules or processes in the accompanying drawings are not necessarily necessary for the implementation of the present invention. Those skilled in the art will understand that the modules in the devices in the implementation scenario can be distributed in the devices of the implementation scenario according to the description of the implementation scenario, or can be changed accordingly and located in one or more devices different from the present implementation scenario. The modules of the above-mentioned implementation scenarios can be combined into one module, or can be further split into multiple sub-modules. The above-mentioned serial numbers of the present invention are only for description and do not represent the advantages and disadvantages of the implementation scenarios. Disclosed above are only a few specific implementation scenarios of the present invention, but the present invention is not limited thereto, and any changes that can be thought of by those skilled in the art should fall within the scope of protection of the present invention.

Claims

1. A method for predicting the ability of an agonist molecule to activate GPCR, characterized in that: include: Obtain the initial ligand structure from the agonist, obtain the inactive receptor structure and activated receptor structure of a given GPCR, and determine whether the structure of the complex formed by the ligand structure and the GPCR is known; If it is unknown, the global molecular docking technology is used to blindly guess the activated complex structure and the inactivated complex structure formed by the ligand structure combined with the activated receptor structure and the inactivated receptor structure, and local structure optimization is performed; Based on the preset enhanced sampling algorithm, a path that can activate the inactive state complex structure into an activated state complex structure is extracted to obtain the initial path of the activation process; Using a preset automatic path optimization algorithm to find a minimum free energy path that is closest to the initial path; The free energy distribution curve along the minimum free energy path is calculated using umbrella sampling. The energy barrier height and the free energy difference before and after activation are determined through the free energy distribution curve, so as to judge the activation ability of the ligand structure on GPCR.

2. The prediction method according to claim 1, characterized in that: The partial structural information of the minimum free energy path is used to optimize the ligand structure in order to increase or decrease the ability of the agonist molecule to activate a given GPCR.

3. The prediction method according to claim 1, characterized in that: The preset automatic path optimization algorithms include: an automatic path search method based on the traveling salesman problem, a finite temperature chord method, a fast tomography method, and a transition path sampling method.

4. The prediction method according to claim 1, characterized in that: The preset enhanced sampling algorithms include coarse-grained model, high-temperature simulation algorithm, stretching dynamics algorithm, targeted dynamics algorithm, and data-driven path algorithm.

5. The prediction method according to claim 3, characterized in that: The automated path search method based on the traveling salesman problem specifically includes: The RMSD between any two structures is defined by all the heavy atoms of the complex structure, and then the PCV is calculated. The high-dimensional space is projected onto a temporary coordinate system based on the current initial path through PCV. The minimum free energy path fragment near the initial path is quickly found through orthogonal direction sampling; a virtual node with a distance of zero from the known path node is set up, and the path nodes on all the minimum free energy path fragments including the virtual node are solved by the traveling salesman method, and the new initial path is obtained after the virtual node is removed and sorted; Correct the spatial resolution of the new initial path and complete a round of iteration; When the change amplitude of the obtained initial path is lower than the preset amplitude, the iteration is stopped and the initial path is taken as the minimum free energy path.

6. The prediction method according to claim 1, characterized in that: For the molecular docking structures of the activated state complex and the inactivated state complex obtained by blind guessing, short-term molecular dynamics simulation was used to perform local structural optimization to obtain an accurate two-state complex structure model.

7. The prediction method according to claim 1, characterized in that: External forces are applied to all heavy atoms of the inactive state complex structure to pull the inactive state complex structure to the active state complex structure model, thereby generating the initial path of the activation process.

8. The prediction method according to claim 1, characterized in that: The inactive and activated receptor structures of a given GPCR are obtained using known structures, homology modeling or AlphaFold.

9. The prediction method according to claim 2, characterized in that: If the current ligand structure does not meet expectations in its ability to activate GPCR, the ligand structure is optimized and modified based on the transition state of the minimum free energy path and the completed structure of the intermediate metastable state to obtain a new ligand structure, and the activation ability of the new ligand structure on GPCR is re-predicted.

10. A prediction system for the ability of agonist molecules to activate GPCR, characterized in that: Includes the following modules: A preparation module is used to obtain the initial ligand structure from the agonist, obtain the inactive receptor structure and the activated receptor structure of a given GPCR, and determine whether the structure of the complex formed by the ligand structure and the GPCR is known; The blind guessing module is used to blindly guess the activated state complex structure and the inactivated state complex structure formed by the ligand structure combined with the activated state receptor structure and the inactivated state receptor structure respectively by using the global molecular docking technology when the complex structure is unknown, and perform local structure optimization; An initial path module, used for extracting a path that can activate the inactive state complex structure into an active state complex structure based on a preset enhanced sampling algorithm, and obtaining an initial path of the activation process; A free energy path module, using a preset automatic path optimization algorithm to find a minimum free energy path that is closest to the initial path; The activation capacity prediction module is used to calculate the free energy distribution curve along the minimum free energy path using umbrella sampling, and determine the energy barrier height and activation front through the free energy distribution curve. The free energy difference after binding is used to determine the activation ability of the ligand structure on GPCR.

Images