A protein ligand binding prediction method based on anisotropic coarse-grained force field

Abstract

The application discloses a protein ligand binding prediction method based on an anisotropic coarse-grained force field, and relates to the technical field of bioengineering; a candidate ligand molecule library is constructed; an amino acid sequence of a callosobruchus chinensis chemical sensing protein is acquired to construct a three-dimensional structure model; coarse-grained molecular dynamics simulation is performed on the callosobruchus chinensis chemical sensing protein and each candidate ligand in the candidate ligand molecule library by using the anisotropic coarse-grained force field, the binding affinity of each ligand to the callosobruchus chinensis chemical sensing protein is determined, and initial effective ligands are determined; the target free energy of each initial effective ligand is calculated; the initial effective ligand with the lowest target free energy is recorded as a target effective ligand; the target effective ligand and the callosobruchus chinensis chemical sensing protein constitute an optimal complex; the optimal complex is used to predict the binding mode, the binding stability and the key binding site of the protein ligand, and the protein ligand binding prediction result is determined. The efficiency and the accuracy of the protein ligand binding prediction are improved.

Description

Technical Field

[0001] This invention belongs to the field of biomedical engineering technology, specifically relating to a method for predicting protein ligand binding based on an anisotropic coarse-grained force field. Background Technology

[0002] From a bioengineering perspective, insects rely on their olfactory system to complete life activities such as mating, finding hosts, foraging, and avoiding predators. Understanding the olfactory recognition mechanism of insects is of great significance for the green control of pests. Chemosensory proteins (CSPs), as soluble olfactory proteins, participate in the recognition and transport of odor molecules and are important targets for developing insect behavior disruptors.

[0003] The bean weevil is a global pest of stored grains, primarily damaging legumes such as mung beans. Current control methods mainly rely on chemical insecticides, but long-term use can damage the ecological environment, harm natural enemies, and induce pesticide resistance in pests. Developing behaviorally disruptive green control agents by utilizing the interaction between the bean weevil's olfactory proteins and host volatiles is an important research direction. Predicting the binding of bean weevil chemoreceptor proteins to host volatiles and elucidating the molecular recognition mechanism can provide theoretical support for site-directed mutagenesis and the rational design of behaviorally disruptive agents.

[0004] Existing methods for predicting protein ligand binding mainly include molecular docking virtual screening, all-atom molecular dynamics simulation, traditional coarse-grained force field simulation, and Metadynamics enhanced sampling simulation. Molecular docking screening is fast, but it is difficult to fully consider the conformational flexibility of proteins, resulting in limited prediction accuracy. All-atom molecular dynamics simulation has high computational accuracy, but it is computationally expensive and has low conformational sampling efficiency, making it unsuitable for large-scale virtual screening. Traditional coarse-grained force field simulation has low computational overhead, but it generally uses isotropic Lennard-Jones potentials, which can easily lead to biases in the prediction of binding modes. Summary of the Invention

[0005] The purpose of this invention is to solve the problem that existing protein ligand binding prediction methods cannot balance computational efficiency and prediction accuracy, and to propose a protein ligand binding prediction method based on an anisotropic coarse-grained force field.

[0006] This invention proposes a protein ligand binding prediction method based on an anisotropic coarse-grained force field, the method comprising: The amino acid sequence of the chemoreceptor protein of the mung bean weevil was obtained, and a three-dimensional structural model was constructed based on the amino acid sequence. A library of candidate ligand molecules was constructed based on the volatile components of the host plants or food sources of the insect species to which the mung bean weevil chemosensory protein belongs. Based on the three-dimensional structural model, coarse-grained molecular dynamics simulations were performed on the chemosensory protein of the mung bean weevil and each candidate ligand in the candidate ligand library using an anisotropic coarse-grained force field to determine the binding affinity of each ligand to the chemosensory protein of the mung bean weevil. Each candidate ligand ranked by binding affinity as the top N ligands is used as the initial effective ligands; Calculate the target free energy for each initial effective ligand; The initial effective ligand with the lowest target free energy is denoted as the target effective ligand. Based on the formation of an optimal complex between the target effective ligand and the mung bean weevil chemoreceptor protein; The protein ligand binding prediction results are determined based on the optimal complex.

[0007] Optionally, determining the binding affinity of each ligand to the bean weevil chemoreceptor protein includes: Each amino acid residue in the three-dimensional structural model is mapped to a coarse-grained particle, and a characteristic axis is predefined for the particles belonging to the side chain in the coarse-grained particles; the aromatic ring side chain site is characterized by the normal vector of the aromatic ring plane, the hydrophobic tail side chain site is characterized by the side chain extension direction, and the polar side chain site is characterized by the direction of the lone pair electrons of the heteroatom in the side chain. Each candidate ligand is subjected to coarse-grained mapping to obtain a coarse-grained representation of each candidate ligand. The coarse-grained representations of each candidate ligand were combined with the coarse-grained particle representation of the mung bean weevil chemosensory protein to construct a coarse-grained complex model. Based on the coarse-grained complex model, coarse-grained molecular dynamics simulations are performed to generate coarse-grained molecular dynamics simulation trajectories for each coarse-grained complex model. Based on the coarse-grained molecular dynamics simulation trajectory, the binding probability and centroid distance distribution of each candidate ligand with the mung bean weevil chemosensory protein were statistically analyzed. The binding affinity of each candidate ligand to the mung bean weevil chemosensory protein was calculated based on the binding probability and centroid distance distribution.

[0008] Optionally, coarse-grained molecular dynamics simulations include: Based on the relative distance between all coarse-grained particles, the isotropic Lennard-Jones-Gaussian basic potential is invoked, and based on the relative orientation between the side chain sites, the angle modulation factor is invoked to generate the coarse-grained molecular dynamics simulation trajectory of each coarse-grained complex model. The angle modulation factor includes the π-π stacking orientation modulation factor between aromatic rings, the hydrophobic tail orientation modulation factor, and the hydrogen bond direction modulation factor.

[0009] Optionally, the expression for the isotropic Lennard-Jones-Gaussian basis potential is: pass Calculate the isotropic Lennard-Jones-Gaussian basic potential; in, Let be the real-time centroid distance between coarse-grained particles i and j, n be the repulsive power exponent, and m be the attractive power exponent. Let be the potential well depth of the k-th particle pair (i, j). Let be the van der Waals diameter of the k-th particle pair (i, j). Let be the magnitude of the Gaussian term for the k-th particle pair (i, j), and be a dimensionless coefficient. The characteristic distance where the Gaussian term center is located is the Gaussian peak. This is the width parameter for the Gaussian term.

[0010] Optionally, calculating the target free energy for each initial effective ligand includes: A whole-atom complex model was constructed for each initial effective ligand and the mung bean weevil chemosensory protein, and the whole-atom molecular dynamics equilibrium simulation was performed to obtain the whole-atom molecular dynamics simulation trajectory. For each all-atom molecular dynamics simulation trajectory, multiple conformational snapshots are extracted at preset time intervals, and the free energy parameters of each conformational snapshot are calculated; the free energy parameters include gas-phase interaction energy, polar solvation energy, and nonpolar solvation energy; The target free energy is obtained by averaging and combining the free energy parameters of the all-atom molecular dynamics simulation trajectory.

[0011] By extracting conformational snapshots, three types of free energy parameters—gas-phase interaction energy, polar solvation energy, and nonpolar solvation energy—are calculated. The target free energy is then obtained by averaging and combining these parameters. Full-atom simulation restores the fine atomic interaction details of the protein ligand system, compensating for the accuracy loss caused by model simplification in coarse-grained simulations. Decomposing multiple free energy parameters and statistically averaging the trajectory snapshots reduces the interference of simulation random noise on the free energy calculation results, yielding more reliable binding free energy evaluation results. This enables high-precision rescreening of initial candidate ligands, improving the accuracy of subsequent optimal complex selection.

[0012] Optionally, the gas-phase interaction energy is the molecular mechanical gas-phase energy, including the van der Waals interaction energy and electrostatic interaction energy between protein ligands; the polar solvation energy is the polar solvation contribution calculated by the generalized Born model; and the nonpolar solvation energy is the nonpolar solvation contribution calculated based on the solvent accessible surface area.

[0013] Optionally, determining the protein ligand binding prediction results based on the optimal complex includes: Using the equilibrium structure obtained from the all-atom molecular dynamics equilibrium simulation of the optimal complex as the initial conformation, the simulated trajectory is obtained by performing all-atom molecular dynamics enhanced sampling simulation. Multiple conformational samples were collected based on the simulated trajectory. The dihedral angle of the backbone of each amino acid residue in each conformational sample was extracted as the input feature and then encoded to construct a training dataset. A hyperspherical variational autoencoder is trained based on the training dataset. The trained hyperspherical variational autoencoder is then embedded into a molecular dynamics enhanced sampling plugin. The three-dimensional latent space coordinates output by the hyperspherical variational autoencoder are used as collective variables to run a Metadynamics simulation and construct a landscape combining free energy. Based on the aforementioned combined free energy landscape, the conformations corresponding to the global energy minimum point and local minimum value are extracted; Calculate the root mean square deviation of the skeleton between each conformation and the relative free energy difference of each local minimum relative to the global minimum energy point; If the root mean square deviation of the skeleton is greater than the preset deviation threshold and there is a relative free energy difference that is less than the preset energy difference, then the target effective ligand is determined to have conformational polymorphism. The conformation of conformational polymorphism is used to determine the predicted protein ligand binding results.

[0014] The hyperspherical variational autoencoder automatically learns collective variables from simulated trajectories, overcoming the subjectivity inherent in traditional Metadynamics simulations that rely on human experience to select collective variables. It employs a von Mises-Fisher distribution to adapt to the periodicity of dihedral angles, avoiding the boundary discontinuity problem of traditional Euclidean latent space. By constructing a binding free energy landscape, it can obtain the binding patterns of protein ligands, quantitatively analyze the dynamic conformational changes of binding pockets and the free energy differences between conformational states, and identify key binding sites to achieve a quantitative assessment of binding stability. This provides a comprehensive and quantitative basis for elucidating the ligand recognition mechanism, site-directed mutagenesis, and ligand molecular design of the bean weevil chemosensory protein.

[0015] Optionally, the dihedral angles of the backbone include the [missing information - likely referring to amino acid residues]. The angle and ψ angle; the encoding method adopts sin / cos encoding to convert each dihedral angle into a two-dimensional vector; the hyperspherical variational autoencoder adopts the von Mises-Fisher distribution as the prior distribution of the latent space, and the latent space is a three-dimensional unit sphere.

[0016] Optionally, the conformational effects of conformational polymorphism can be used to determine protein ligand binding predictions, including: For each conformation reduction to obtain the all-atom three-dimensional conformation, protein ligand interaction information is determined, and the amino acid residues involved in the interaction are recorded to obtain static key binding residues; the protein ligand interaction information includes hydrogen bonds, hydrophobic interactions, and π-π stacking interactions between protein ligands. The structural overlap comparison of the three-dimensional conformations of all atoms obtained by complete reduction was performed to determine the dynamic key residues. By identifying broken or newly formed hydrogen bonds and hydrophobic contact pairs between different conformations, a dynamic interaction network can be obtained. Key binding sites are obtained based on static and dynamic key binding residues; The binding mode of protein ligands is determined by combining the protein-ligand interaction information corresponding to each conformation and the dynamic interaction network.

[0017] By restoring the abstract free energy information of the potential space to the real three-dimensional atomic structure, we can quantitatively evaluate the thermodynamic stability of different conformational states and identify conformational polymorphism in the ligand binding process. At the same time, we can obtain the residues responsible for static binding as well as the dynamic key residues involved in pocket conformational transitions, distinguish between static interactions and dynamic interaction networks in the conformational transition process, and avoid mechanism misjudgments caused by relying solely on a single lowest energy conformation. This enables multi-dimensional quantitative analysis of the ligand binding mode and binding site of the chemoreceptor protein of the bean weevil, providing static and dynamic dual-layer residue target information for site-directed mutagenesis experiments and the rational design of behavioral disruptors.

[0018] The beneficial effects of this invention are: 1) This invention first obtains the amino acid sequence of the chemosensory protein of the mung bean weevil and constructs a three-dimensional structural model. Then, it constructs a candidate ligand molecular library of host volatiles and their analogues. Using an anisotropic coarse-grained force field, it conducts coarse-grained molecular dynamics simulation to complete the initial screening and obtain the initial effective ligands. The target free energy of each initial effective ligand is further calculated to screen the target effective ligands. After forming the optimal complex, the binding mode, binding stability and key binding sites are predicted. 2) This invention determines the connection relationship between proteins and ligands by constructing an anisotropic coarse-grained force field, enabling bioengineering schemes to quantitatively describe directional interactions such as π-π stacking and hydrogen bonding, balancing simulation efficiency and prediction accuracy, avoiding the subjective bias caused by the manual selection of collective variables in traditional enhanced sampling, effectively analyzing the dynamic conformational changes of protein binding pockets, and efficiently realizing the virtual screening and binding mechanism analysis of chemosensory proteins and ligands of the bean weevil, providing a reliable basis for the rational design of pest behavior interference agents and site-directed mutagenesis experiments. Attached Figure Description

[0019] The invention will now be further described with reference to the accompanying drawings.

[0020] Figure 1 A flowchart illustrating a protein ligand binding prediction method based on an anisotropic coarse-grained force field, provided for embodiments of the present invention; Figure 2 This is a flowchart of a protein ligand binding prediction method provided in an embodiment of the present invention. Detailed Implementation

[0021] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0022] It should be noted that all formula calculations in the scheme are purely numerical calculations.

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

[0024] This invention provides a method for predicting protein ligand binding based on an anisotropic coarse-grained force field. See also... Figure 1 The method includes the following steps: S101, Obtain the amino acid sequence of the chemoreceptor protein of the mung bean weevil, and construct a three-dimensional structural model based on the amino acid sequence; S102, a candidate ligand molecule library was constructed based on the volatile components of the host plants or food sources of the insect species to which the mung bean weevil chemosensory protein belongs; S103, based on a three-dimensional structural model, uses an anisotropic coarse-grained force field to perform coarse-grained molecular dynamics simulations on the chemosensory protein of the bean weevil and each candidate ligand in the candidate ligand molecular library to determine the binding affinity of each ligand to the chemosensory protein of the bean weevil. S104, obtain each candidate ligand ranked in the top N by binding affinity as the initial effective ligand; S105, Calculate the target free energy of each initial effective ligand; S106, the initial effective ligand with the lowest target free energy is denoted as the target effective ligand; S107 is based on the formation of an optimal complex between the target effective ligand and the green bean weevil chemoreceptor protein. S108, based on the optimal complex, determines the predicted protein ligand binding results.

[0025] In one implementation, AlphaFold2 is used to predict the three-dimensional structural model of the chemoreceptor protein of the bean weevil, and key aromatic residues and hydrophobic residues in the inner wall of the binding pocket of the three-dimensional structural model are identified; the volatile components from the host plant or food source include at least one of hexanal, methyl heptenone, n-nonanal, decanal, and trans-3-hexen-1-ol.

[0026] In one implementation, the top N positions are determined by technical personnel.

[0027] In one implementation, determining the isotropic Lennard-Jones-Gaussian basic potential includes: pass Calculate the isotropic Lennard-Jones-Gaussian basic potential; in, Let be the real-time centroid distance between coarse-grained particles i and j, n be the repulsive power exponent, and m be the attractive power exponent. Let be the potential well depth of the k-th particle pair (i, j). Let be the van der Waals diameter of the k-th particle pair (i, j). Let be the magnitude of the Gaussian term for the k-th particle pair (i, j), and be a dimensionless coefficient. The characteristic distance where the Gaussian term center is located is the Gaussian peak. The width parameter of the Gaussian term. It is a Gaussian function.

[0028] In one implementation, n is greater than m, and the specific value is determined by the technicians; the total potential energy of the anisotropic coarsened force field is determined by multiplying the isotropic Lennard Jones Gaussian basic potential by the angle modulation factor.

[0029] In one implementation, the complete form of the anisotropic generalized Lennard-Jones-Gaussian potential function includes: (via the formula...) Indicates; among which, Let be the distance between particle i and particle j. Let be the set of relative orientations of particle i and particle j, representing the orientation relationship between particles. The weighting coefficients and , Let be the angle modulation factor for the k-th type of interaction. It is an isotropic Lennard-Jones-Gaussian basic potential.

[0030] In one implementation, an anisotropic coarse-grained force field represents the atomic set in the molecular system as coarse-grained particles through a coarse-grained mapping rule. The bonding interaction energy, non-bonding interaction energy (isotropic Lennard-Jones-Gaussian basis potential), and electrostatic interaction energy between these coarse-grained particles are calculated using a bonding potential function, a non-bonding potential function, and an electrostatic potential function, respectively. An anisotropic coarse-grained force field is constructed based on the coarse-grained particles, their bonding interaction energy, non-bonding interaction energy, and electrostatic interaction energy. Then, coarse-grained molecular dynamics simulation determines the coarse-grained molecular dynamics simulation trajectory using the anisotropic coarse-grained force field. The non-bonding potential function is an anisotropic generalized Lennard-Jones-Gaussian potential function. The bonding potential function is used to describe the interaction between two or more coarse-grained particles connected by chemical bonds within the same molecule. The electrostatic potential function is used to describe the long-range electrostatic interaction between partially charged coarse-grained particles, which is an important long-range force in protein-ligand binding.

[0031] In one embodiment, based on a three-dimensional structural model, coarse-grained molecular dynamics simulations are performed on the chemoreceptor protein of the bean weevil and each candidate ligand in the candidate ligand library using an anisotropic coarse-grained force field to determine the binding affinity of each ligand to the bean weevil chemoreceptor protein, including: Each amino acid residue in the three-dimensional structural model is mapped to a coarse-grained particle, and characteristic axes are predefined for particles belonging to the side chain in the coarse-grained particles; the aromatic ring side chain site is characterized by the normal vector of the aromatic ring plane, the hydrophobic tail side chain site is characterized by the direction of side chain extension, and the polar side chain site is characterized by the direction of the lone pair electrons of the heteroatom in the side chain. Each candidate ligand is subjected to coarse-grained mapping to obtain a coarse-grained representation of each candidate ligand. The coarse-grained representations of each candidate ligand were combined with the coarse-grained particle representation of the mung bean weevil chemosensory protein to construct a coarse-grained complex model. Based on the coarse-grained complex model, coarse-grained molecular dynamics simulations are performed. The isotropic Lennard-Jones-Gaussian basis potential is invoked according to the relative distance between all coarse-grained particles, and the angle modulation factor is invoked according to the relative orientation between side chain sites to generate the coarse-grained molecular dynamics simulation trajectory of each coarse-grained complex model. The angle modulation factor includes the π-π stacking orientation modulation factor between aromatic rings, the hydrophobic tail orientation modulation factor, and the hydrogen bond directionality modulation factor. Based on coarse-grained molecular dynamics simulation trajectories, the binding probability and centroid distance distribution of each candidate ligand with the chemoreceptor protein of the bean weevil were statistically analyzed. The binding affinity of each candidate ligand to the chemoreceptor protein of the bean weevil was calculated based on the binding probability and centroid distance distribution.

[0032] In one implementation, each amino acid residue in the three-dimensional structural model is mapped to a coarse-grained particle. Each coarse-grained particle includes a main chain site and a side chain site, with glycine only mapped to the main chain site. The coarse-grained mapping is performed based on the functional group type of each candidate ligand: aromatic ring structures in the candidate ligand molecule are mapped to aromatic ring particles, alkyl chains or hydrophobic groups are mapped to hydrophobic tail particles, and polar groups such as hydroxyl, carbonyl, or amino groups are mapped to polar particles. After this mapping, each candidate ligand is represented by several coarse-grained particles, each particle possessing a type attribute corresponding to the protein's side chain site, thereby ensuring that matching angle modulation factors can be invoked in subsequent simulations.

[0033] In one implementation, each candidate ligand is placed at the entrance of the binding pocket of the bean weevil chemosensory protein or inside the pocket as the initial conformation. Then, each complex model is placed in a simulation box, coarse-grained solvent particles are added and ions are added to neutralize the system charge. Energy minimization, NVT ensemble equilibrium, and NPT ensemble equilibrium are performed sequentially to eliminate unreasonable contact and bring the system to the target temperature and pressure. Based on each coarse-grained complex model after equilibrium, coarse-grained molecular dynamics production simulation is performed to generate the coarse-grained molecular dynamics simulation trajectory of each coarse-grained complex model.

[0034] In one embodiment, calculating the target free energy for each initial effective ligand includes: For each initial effective ligand, a whole-atom complex model was constructed with the mung bean weevil chemosensory protein, and the whole-atom molecular dynamics equilibrium simulation was performed to obtain the whole-atom molecular dynamics simulation trajectory. For each all-atom molecular dynamics simulation trajectory, multiple conformational snapshots are extracted at preset time intervals, and the free energy parameters of each conformational snapshot are calculated; the free energy parameters include gas-phase interaction energy, polar solvation energy, and nonpolar solvation energy; The target free energy is obtained by averaging and combining the free energy parameters of the all-atom molecular dynamics simulation trajectory.

[0035] In one implementation, the preset time interval is determined by technical personnel.

[0036] In one implementation, the gas-phase interaction energy is the molecular mechanical gas-phase energy, which includes the van der Waals interaction energy and electrostatic interaction energy between protein ligands; the polar solvation energy is the polar solvation contribution calculated by the generalized Born model; and the nonpolar solvation energy is the nonpolar solvation contribution calculated based on the solvent accessible surface area.

[0037] In one embodiment, see Figure 2 A flowchart of a protein ligand binding prediction method is provided, specifically including: S1081, using the equilibrium structure obtained from the all-atom molecular dynamics equilibrium simulation of the optimal complex as the initial conformation, the simulated trajectory is obtained by all-atom molecular dynamics enhanced sampling simulation; S1082: Multiple conformational samples were collected based on simulated trajectories. The dihedral angles of the backbone of each amino acid residue in each conformational sample were extracted as input features and then encoded to construct a training dataset. S1083: Train a hyperspherical variational autoencoder based on the training dataset, embed the trained hyperspherical variational autoencoder into the molecular dynamics enhanced sampling plugin, use the three-dimensional latent space coordinates output by the hyperspherical variational autoencoder as the collective variable, run Metadynamics simulation, and construct a landscape combining free energy. S1084, based on the combination of free energy landscape, extracts the conformations corresponding to the global energy minimum point and local minimum value; S1085, calculate the root mean square deviation of the skeleton between each conformation and the relative free energy difference of each local minimum relative to the global minimum energy point; S1086, If the root mean square deviation of the skeleton is greater than the preset deviation threshold and there is a relative free energy difference that is less than the preset energy difference, then the target effective ligand is determined to have conformational polymorphism. S1087, the conformation of conformational polymorphism is used to determine the predicted results of protein ligand binding.

[0038] In one implementation, the Metadynamics simulation employs a well-tempered scheme with an initial Gaussian height of 1.0-2.0 kJ / mol, a Gaussian width of 0.06, a bias factor of 5-15, and a Gaussian potential added every 400 steps. The simulation time is 50-100 ns. The convergence criteria include: the mean Gaussian height drops to less than 10% of the initial height within the last 10,000 steps; the standard deviation of the diffusion coverage of the three collective variables changes by less than 5% within the last 20 ns; and the root mean square deviation of the free energy between time-block landscapes is less than 2 kJ / mol.

[0039] In one implementation, local minima are determined by using a valley-finding algorithm, such as gradient descent or density-based spatial clustering, to scan the extreme points of the combined free energy landscape and extract the CV coordinates corresponding to all minimum points with positive curvature. Based on the CV coordinates, the corresponding three-dimensional conformations are obtained by back-applying them from the simulated trajectory. Subsequently, RMSD-based clustering analysis is performed on each conformation, and the minimum values ​​with RMSD between conformations less than a preset threshold of 1 Å are merged into the same conformation cluster. The preset deviation threshold is determined by technical personnel.

[0040] In one implementation, the dihedral angles of the backbone include the amino acid residues. Angles and ψ; the encoding method uses sin / cos encoding, converting each dihedral angle into a two-dimensional vector (sinθ, cosθ).

[0041] In one implementation, the hyperspherical variational autoencoder uses the von Mises-Fisher distribution as the prior distribution of the latent space, which is a three-dimensional unit sphere.

[0042] In one implementation, the root mean square deviation of the skeleton between each conformation and the relative free energy difference between each local minimum and the global minimum energy point are calculated. If the root mean square deviation of the skeleton between at least one conformation and the global minimum conformation is greater than a preset structural deviation threshold, and the relative free energy difference is less than or equal to a preset energy threshold, then it is determined that the binding of the target effective ligand to the chemoreceptor protein of the bean weevil exhibits conformational polymorphism. Otherwise, it is determined that there is no conformational polymorphism, and only a single all-atom three-dimensional conformation corresponding to the global minimum energy point exists. Hydrogen bond, hydrophobic interaction, and π-π stacking interaction are directly identified for this conformation to obtain static key binding residues, and then the subsequent key binding sites and binding modes are analyzed.

[0043] In one embodiment, determining the protein ligand binding prediction result by acting on the conformation of conformational polymorphism includes: For each conformation reduction, the all-atom three-dimensional conformation is obtained, the protein-ligand interaction information is determined, and the amino acid residues involved in the interaction are recorded to obtain the static key binding residues; the protein-ligand interaction information includes hydrogen bonds, hydrophobic interactions, and π-π stacking interactions between protein ligands. The structural overlap comparison of the three-dimensional conformations of all atoms obtained by complete reduction was performed to determine the dynamic key residues. By identifying broken or newly formed hydrogen bonds and hydrophobic contact pairs between different conformations, a dynamic interaction network can be obtained. Key binding sites are obtained based on static and dynamic key binding residues; The binding mode of protein ligands is determined by combining the protein-ligand interaction information corresponding to each conformation and the dynamic interaction network.

[0044] In one implementation, the preset energy is 4 kJ / mol; the structural overlap comparison of multiple three-dimensional conformations of all atoms obtained by reduction is performed, and residues with main chain or side chain displacement greater than 2.0 Å are identified as dynamic key residues; Å is a unit of length; 1 Å = 0.1 nm.

[0045] In one implementation, after reconstructing the binding free energy landscape from the enhanced sampling technique used for molecular simulations (Metadynamics), the global minimum energy point is located on the three-dimensional latent space grid of this binding free energy landscape. At the same time, all local minima with an energy difference from the global minimum point not greater than a preset energy are selected, and the three-dimensional latent space coordinates corresponding to each minima are extracted. The three-dimensional latent space coordinates of each minima are input into the decoder of the trained hyperspherical variational autoencoder. Through the mapping operation of the decoder, the abstract coordinates of the three-dimensional latent space are restored and output as the three-dimensional conformation of the protein ligands corresponding to each minima.

[0046] In one implementation, the hyperspherical variational autoencoder comprises an encoder and a decoder. The encoder is responsible for compressing and mapping the high-dimensional protein backbone dihedral conformation information to the latent space of a three-dimensional unit sphere to obtain the three-dimensional latent space coordinates. The decoder is a symmetric inverse network structure of the encoder, which receives the three-dimensional latent space coordinates as input, re-decodes the low-dimensional latent space representation in reverse to restore the high-dimensional conformation features of the protein backbone dihedral, and then recovers the full-atom three-dimensional conformation of the protein ligand corresponding to each minimum point.

[0047] In one implementation, the specific method for determining the binding stability of protein ligands based on the relative free energy difference is as follows: the relative free energy difference represents the thermodynamic energy cost of the local conformation compared to the globally most stable conformation; the smaller the relative free energy difference, the higher the proportion of this conformation in the system, and the stronger the conformation stability; the larger the relative free energy difference, the lower the probability of this conformation appearing, and the weaker the stability. If the relative free energy difference of the local minimum conformation is greater than or equal to a preset energy, then this local conformation is thermodynamically unstable, and the system mainly uses the conformation corresponding to the lowest global energy point as the dominant binding conformation.

[0048] In one implementation, the specific operation of determining dynamic key residues involves structural overlap comparison of all reduced full-atom three-dimensional conformations: All reduced protein ligand full-atom three-dimensional conformations are spatially aligned using the protein backbone atoms as a reference. After alignment, the atomic spatial coordinates of each amino acid residue in different conformations are compared, and the spatial displacements of the main chain and side chain atoms are calculated. Amino acid residues with a main chain or side chain atomic spatial displacement greater than 2.0 Å are identified as dynamic key residues. These residues undergo significant spatial position changes during conformational transitions, participating in the opening and closing conformational transitions of the protein binding pocket. By comparing the interactions between different conformations, newly generated hydrogen bonds and hydrophobic contact pairs, as well as broken hydrogen bonds and hydrophobic contact pairs, are identified during conformational switching, resulting in a dynamic interaction network. This dynamic interaction network reflects the reconstruction and disappearance patterns of intermolecular interactions during conformational transitions.

[0049] In one implementation, the key binding site is obtained based on static and dynamic key binding residues in the following way: all identified static and dynamic key binding residues are merged, and duplicate amino acid residues are removed; the resulting set of residues is the key binding site; the static key binding residues embody the static basis of ligand binding, the dynamic key residues embody the structural basis of pocket conformational changes, and the key binding site contains residue information related to both ligand binding and conformational changes; by comprehensively considering the protein-ligand interaction information corresponding to each conformation, the dynamic interaction network corresponding to conformational changes, and the relative free energy difference between different conformational states, the binding mode of the protein ligand is determined, and the mechanism of the protein-ligand binding system is elucidated.

[0050] The foregoing has described one embodiment of the present invention in detail, but this content is merely a preferred embodiment and should not be considered as limiting the scope of the present invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the scope of the claims of this invention.

Claims

1. A method for predicting protein ligand binding based on anisotropic coarse-grained force fields, characterized in that, The method includes: The amino acid sequence of the chemoreceptor protein of the mung bean weevil was obtained, and a three-dimensional structural model was constructed based on the amino acid sequence. A library of candidate ligand molecules was constructed based on the volatile components of the host plants or food sources of the insect species to which the mung bean weevil chemosensory protein belongs. Based on the three-dimensional structural model, coarse-grained molecular dynamics simulations were performed on the chemosensory protein of the mung bean weevil and each candidate ligand in the candidate ligand library using an anisotropic coarse-grained force field to determine the binding affinity of each ligand to the chemosensory protein of the mung bean weevil. Each candidate ligand ranked by binding affinity as the top N ligands is used as the initial effective ligands; Calculate the target free energy for each initial effective ligand; The initial effective ligand with the lowest target free energy is denoted as the target effective ligand. Based on the formation of an optimal complex between the target effective ligand and the mung bean weevil chemoreceptor protein; The protein ligand binding prediction results are determined based on the optimal complex.

2. The protein ligand binding prediction method based on anisotropic coarse-grained force field according to claim 1, characterized in that, The binding affinity of each ligand to the chemoreceptor protein of the bean weevil was determined as follows: Each amino acid residue in the three-dimensional structural model is mapped to a coarse-grained particle, and a characteristic axis is predefined for the particles belonging to the side chain in the coarse-grained particles; the aromatic ring side chain site is characterized by the normal vector of the aromatic ring plane, the hydrophobic tail side chain site is characterized by the side chain extension direction, and the polar side chain site is characterized by the direction of the lone pair electrons of the heteroatom in the side chain. Each candidate ligand is subjected to coarse-grained mapping to obtain a coarse-grained representation of each candidate ligand. The coarse-grained representations of each candidate ligand were combined with the coarse-grained particle representation of the mung bean weevil chemosensory protein to construct a coarse-grained complex model. Based on the coarse-grained complex model, coarse-grained molecular dynamics simulations are performed to generate coarse-grained molecular dynamics simulation trajectories for each coarse-grained complex model. Based on the coarse-grained molecular dynamics simulation trajectory, the binding probability and centroid distance distribution of each candidate ligand with the mung bean weevil chemosensory protein were statistically analyzed. The binding affinity of each candidate ligand to the mung bean weevil chemosensory protein was calculated based on the binding probability and centroid distance distribution.

3. The protein ligand binding prediction method based on anisotropic coarse-grained force field according to claim 2, characterized in that, Performing coarse-grained molecular dynamics simulations includes: Based on the relative distance between all coarse-grained particles, the isotropic Lennard-Jones-Gaussian basic potential is invoked, and based on the relative orientation between the side chain sites, the angle modulation factor is invoked to generate the coarse-grained molecular dynamics simulation trajectory of each coarse-grained complex model. The angle modulation factor includes the π-π stacking orientation modulation factor between aromatic rings, the hydrophobic tail orientation modulation factor, and the hydrogen bond direction modulation factor.

4. The protein ligand binding prediction method based on anisotropic coarse-grained force field according to claim 3, characterized in that, The expression for the isotropic Lennard-Jones-Gaussian basic potential is: pass Calculate the isotropic Lennard-Jones-Gaussian basic potential; in, Let be the real-time centroid distance between coarse-grained particles i and j, n be the repulsive power exponent, and m be the attractive power exponent. Let be the potential well depth of the k-th particle pair (i, j). Let be the van der Waals diameter of the k-th particle pair (i, j). Let be the magnitude of the Gaussian term for the k-th particle pair (i, j), and be a dimensionless coefficient. The characteristic distance where the Gaussian term center is located is the Gaussian peak. This is the width parameter for the Gaussian term.

5. The protein ligand binding prediction method based on anisotropic coarse-grained force field according to claim 1, characterized in that, Based on the optimal complex, the protein ligand binding prediction results include: Using the equilibrium structure obtained from the all-atom molecular dynamics equilibrium simulation of the optimal complex as the initial conformation, the simulated trajectory is obtained by performing all-atom molecular dynamics enhanced sampling simulation. Multiple conformational samples were collected based on the simulated trajectory. The dihedral angle of the backbone of each amino acid residue in each conformational sample was extracted as the input feature and then encoded to construct a training dataset. A hyperspherical variational autoencoder is trained based on the training dataset. The trained hyperspherical variational autoencoder is then embedded into a molecular dynamics enhanced sampling plugin. The three-dimensional latent space coordinates output by the hyperspherical variational autoencoder are used as collective variables to run a Metadynamics simulation and construct a landscape combining free energy. Based on the aforementioned combined free energy landscape, the conformations corresponding to the global energy minimum point and local minimum value are extracted; Calculate the root mean square deviation of the skeleton between each conformation and the relative free energy difference of each local minimum relative to the global minimum energy point; If the root mean square deviation of the skeleton is greater than the preset deviation threshold and there is a relative free energy difference that is less than the preset energy difference, then the target effective ligand is determined to have conformational polymorphism. The conformation of conformational polymorphism is used to determine the predicted results of protein ligand binding.

6. The protein ligand binding prediction method based on anisotropic coarse-grained force field according to claim 5, characterized in that, The results of conformational polymorphism analysis to determine protein ligand binding predictions include: For each conformation reduction to obtain the all-atom three-dimensional conformation, protein ligand interaction information is determined, and the amino acid residues involved in the interaction are recorded to obtain static key binding residues; the protein ligand interaction information includes hydrogen bonds, hydrophobic interactions, and π-π stacking interactions between protein ligands. The structural overlap comparison of the three-dimensional conformations of all atoms obtained by complete reduction was performed to determine the dynamic key residues. By identifying broken or newly formed hydrogen bonds and hydrophobic contact pairs between different conformations, a dynamic interaction network can be obtained. Key binding sites are obtained based on static and dynamic key binding residues; The binding mode of protein ligands is determined by combining the protein-ligand interaction information corresponding to each conformation and the dynamic interaction network.

Images