Disease marker binding pocket identification method

Abstract

The invention provides a disease marker binding pocket identification method, and belongs to the technical field of disease marker binding pockets, and the method comprises the following steps: firstly, constructing a simulation system containing a disease marker, and carrying out 100-500ns molecular dynamics simulation; a conformation change track is collected and subjected to normalization processing; carrying out decomposition calculation on the sampled conformation, and analyzing a conformation cluster with relatively high frequency; calculating structural parameters of amino acid residues, and screening a region with relatively high flexibility as a potential binding site; analyzing the solvent accessibility and electrostatic characteristics of the sites; identifying a temporary binding pocket through the water molecule distribution density; performing principal component analysis to extract a conformation change mode; constructing an amino acid residue interaction network, and screening a stable binding pocket based on indexes such as binding free energy, conformational stability and pocket comprehensive performance evaluation; the method solves the problem that most existing methods only consider static structural characteristics of molecules and neglect the influence of dynamic change of protein conformation on binding pocket formation.

Description

Technical Field

[0001] The present invention belongs to the technical field of disease biomarker binding pockets, and more particularly, relates to a method for identifying disease biomarker binding pockets. Background Art

[0002] Disease biomarkers are a class of biomolecules that are abnormally expressed during the occurrence and development of diseases and participate in the pathological processes of diseases through interactions with other molecules. Accurately identifying the binding pockets on disease biomarker molecules is of great significance for understanding disease mechanisms and developing targeted drugs. Currently, the methods for identifying binding pockets mainly include two categories: experimental methods and computational methods. Experimental methods such as X-ray crystallography, nuclear magnetic resonance spectroscopy, etc. can obtain high-resolution structural information, but they have a long experimental period, high cost, and it is difficult to capture the dynamic characteristics of molecules. Computational methods mainly identify based on geometric features and energy evaluation, such as algorithms like POCKET, LIGSITE, fpocket, etc. These methods have high computational efficiency, but limited accuracy.

[0003] Most of the existing computational methods only consider the static structural features of molecules and ignore the influence of the dynamic changes in protein conformation on the formation of binding pockets. In the physiological environment, disease biomarker molecules are in continuous conformational changes, and the binding pockets may appear, disappear, or change shape with conformational changes. This results in the prediction results may deviate from the actual situation under physiological conditions. Summary of the Invention

[0004] In view of this, the present invention provides a method for identifying disease biomarker binding pockets, which can solve the problem that most of the existing methods only consider the static structural features of molecules and ignore the influence of the dynamic changes in protein conformation on the formation of binding pockets.

[0005] The present invention is implemented as follows:

[0006] The present invention provides a method for identifying a disease biomarker binding pocket, comprising: constructing a simulation system containing a disease biomarker molecule, water molecules and ions, performing standardized preprocessing on the simulation system, and performing molecular dynamics simulation; collecting the conformational change trajectory of the disease biomarker during the molecular dynamics simulation, sampling the conformational change trajectory, and performing global normalization processing on the collected data; performing conformational decomposition calculation on the sampled disease biomarker conformations, separating stable conformations and variable conformations; calculating the root mean square fluctuation value and the secondary structure change rate of each amino acid residue in the disease biomarker conformation, establishing an amino acid residue structure parameter matrix, and determining the region where the flexibility of the amino acid residue is greater than 0.2 nanometers and the secondary structure change rate is greater than 15% as the potential binding site; performing solvent accessibility analysis on the potential binding site, identifying the high-density region of water molecules, and determining it as the temporary binding pocket; performing principal component analysis and clustering analysis on the temporary binding pocket, establishing a binding pocket feature vector; constructing binding site map data, calculating the binding free energy and conformational stability index of the temporary binding pocket, and selecting the stable binding pocket. Specifically, it includes the following operation steps:

[0007] S10. Construct a simulation system containing a disease biomarker molecule, water molecules and ions, perform standardized preprocessing on the simulation system, and perform molecular dynamics simulation for 100 nanoseconds to 500 nanoseconds through a molecular dynamics simulation software;

[0008] S20. Collect the conformational change trajectory of the disease biomarker during the molecular dynamics simulation, sample the conformational change trajectory at a time interval of 10 picoseconds to 50 picoseconds, and perform global normalization processing on the collected data;

[0009] S30. Perform conformational decomposition calculation on the sampled disease biomarker conformations, separate stable conformations and variable conformations, and select the conformation clusters where the occurrence frequency of the stable conformations is in the top 20% to 30% for structural feature analysis;

[0010] S40. Calculate the root mean square fluctuation value and the secondary structure change rate of each amino acid residue in the disease biomarker conformation, establish an amino acid residue structure parameter matrix, and determine the region where the flexibility of the amino acid residue is greater than 0.2 nanometers and the secondary structure change rate is greater than 15% as the potential binding site;

[0011] S50. Perform solvent accessibility analysis on the potential binding site, calculate the exposed surface area change value and the electrostatic potential distribution, and select the region where the fluctuation range of the exposed surface area is between 0.5 square nanometers and 2.0 square nanometers and the electrostatic potential fluctuation is greater than 3 kcal / mol;

[0012] S60. Analyze the distribution density of solvent molecules in the region, identify the high-density regions where the residence time of water molecules exceeds 500 picoseconds and the distribution density is more than 3 times that of the bulk phase, and determine them as temporary binding pockets;

[0013] S70. Conduct principal component analysis and clustering analysis on the temporary binding pockets, extract the main conformational change patterns with a contribution rate to the total variance of more than 90%, and establish the characteristic vectors of the binding pockets;

[0014] S80. Construct the binding site map data containing the interaction network of amino acid residues, with amino acid residues as nodes and the interactions between residues as edges, calculate the binding free energy and conformational stability index of the temporary binding pockets, and select the sites with a binding free energy lower than -20 kJ / mol and a conformational stability index greater than 0.75 as stable binding pockets; and conduct a comprehensive performance analysis of the candidate binding pockets based on the pocket comprehensive performance evaluation equation, and select the pockets with a comprehensive evaluation index greater than 0.6 as stable binding pockets.

[0015] Among them, the pocket comprehensive performance evaluation equation uses the comprehensive performance score P of the binding site as the evaluation index, and this score is obtained by normalizing and weighting four key parameters: binding free energy ΔG, root mean square fluctuation value RMSF of key residues, relative enrichment factor REF of water molecules, and sequence conservation score CN; among them, the binding free energy reflects the thermodynamic stability of the binding site, the root mean square fluctuation value characterizes the conformational flexibility of the pocket region, the relative enrichment factor of water molecules reflects the solvent accessibility characteristics of the pocket, and the sequence conservation score reflects the evolutionary importance of the pocket; each parameter is normalized using the corresponding standardized reference values ΔG0, RMSF0, REF0, CN0; the normalized parameters are weighted and summed through the weight coefficients α, β, γ, δ, and the sum of the weight coefficients is equal to 1; the range of the comprehensive performance score P is from 0 to 1, and the larger the value, the more stable and reliable the binding pocket is. When P is greater than 0.6, it can be determined as a stable binding pocket.

[0016] The construction of the simulation system containing disease biomarker molecules, water molecules and ions specifically places the disease biomarker molecules at the center of a cubic water box with a side length of 10 nanometers, adds water molecules using the TIP3P water model, neutralizes the charge using chloride ions or sodium ions according to the total charge of the system, adds 0.15 mol / L sodium chloride to simulate the physiological environment, and constructs the complete simulation system using periodic boundary conditions.

[0017] Perform standard preprocessing on the simulation system. Specifically, first perform energy minimization on the system, optimize it for 5000 steps using the steepest descent method, and set the convergence criterion to a maximum force less than 1000 kJ / mol / nm. Perform system equilibration for 100 picoseconds under constant volume and constant temperature conditions, control the temperature at 300 Kelvin, use the Berendsen temperature coupling algorithm, and set the time constant to 0.1 picosecond.

[0018] Perform the molecular dynamics simulation through the molecular dynamics simulation software. Specifically, use the GROMACS software package and perform the molecular dynamics simulation using the CHARMM force field. Set the time step to 2 femtoseconds, use the LINCS algorithm to constrain the bond lengths of all hydrogen-bonded bonds, use the PME method to handle long-range electrostatic interactions, set the cutoff radius for van der Waals interactions to 1.0 nm, and perform the simulation under the NPT ensemble. The pressure is maintained at 1 atmosphere through the Parrinello-Rahman method.

[0019] Collect the conformational change trajectories of the disease biomarker during the molecular dynamics simulation. Specifically, record the conformational information of the system every 20 picoseconds, including atomic coordinates, velocities, and energy data. Use periodic image correction to remove the influence of molecular cross-boundary motion, and perform structural alignment of the protein backbone through the least-squares superposition method to eliminate the influence of overall translation and rotation.

[0020] Perform conformational decomposition calculations on the conformations of the disease biomarker obtained by sampling. Specifically, use the principal component analysis method to perform dimensionality reduction on the conformational clusters, calculate the eigenvectors and eigenvalues of the covariance matrix, select the principal components with a cumulative contribution rate reaching 90% as the main dimensions of the conformational space, use the density peak clustering algorithm to classify the conformations, and determine the clustering centers through local density and distance criteria.

[0021] Calculate the root mean square fluctuation value and secondary structure change rate of each amino acid residue in the conformation of the disease biomarker. Specifically, calculate the standard deviation of the Cα atom position of each amino acid residue based on the aligned trajectory, use the DSSP algorithm to analyze the change in the secondary structure type of each amino acid residue during the simulation, and count the proportion of the time that each amino acid residue maintains its original secondary structure type under different conformations. Take the negative value of it as the secondary structure change rate.

[0022] Perform solvent accessibility analysis on the potential binding sites. Specifically, use the NACCESS program to calculate the solvent accessible surface area of each amino acid residue under different conformations, perform rolling ball calculations using a probe radius of 1.4 Å, use the APBS software to solve the Poisson-Boltzmann equation to calculate the molecular surface electrostatic potential distribution, and set the internal dielectric constant to 2 and the external dielectric constant to 80.

[0023] Perform solvent molecular distribution density analysis in the said region. Specifically, divide the simulation system into cubic grids with a side length of 0.5 Å, count the occurrence frequency of water molecule oxygen atoms in each grid, calculate the enrichment factor relative to the bulk water density, track the residence time of individual water molecules at specific positions, and evaluate the kinetic characteristics of water molecules through continuous residence time and exchange frequency.

[0024] Construct binding site map data containing the interaction network of amino acid residues. Specifically, establish network connections based on van der Waals contacts, hydrogen bonds, salt bridges, and hydrophobic interactions between the amino acid residues. The interaction strength determines the edge weights according to the contact frequency and energy contribution, and use the NetworkX library to calculate the degree distribution, clustering coefficient, and centrality topological parameters of the network.

[0025] Calculate the binding free energy and conformational stability index of the temporary binding pocket. Specifically, use the MM / PBSA method to calculate the polar and non-polar interaction energies in the binding pocket region, consider the contributions of solvation free energy and conformational entropy, calculate the stability index based on the topological parameters and kinetic characteristics of the amino acid residue network in the binding pocket, and comprehensively evaluate the thermodynamic and kinetic properties of the binding site.

[0026] Compared with the prior art, the beneficial effects of a method for identifying a binding pocket of a disease biomarker provided by the present invention are as follows: First, by integrating technologies such as molecular dynamics simulation, conformational analysis, and machine learning, the problems existing in the prior art are effectively solved. First, this method uses long-time-scale molecular dynamics simulation to comprehensively capture the conformational change characteristics of disease biomarker molecules, overcoming the limitations of static structure analysis. Through systematic sampling and analysis of the conformational change trajectories, the dynamic behavior of protein molecules is accurately reflected.

[0027] Secondly, this method innovatively introduces water molecule distribution density analysis. By identifying the residence time and distribution characteristics of water molecules, the solvation characteristics of potential binding sites are revealed. This method not only considers the interactions mediated by water molecules, but also can identify structurally significant water molecule sites with functional significance, providing important references for subsequent drug design.

[0028] The present invention establishes a complete binding pocket evaluation system, including multiple evaluation indicators such as binding free energy and conformational stability index. Through comprehensive analysis of multi-dimensional characteristics, the reliability of the identification results is ensured. In particular, by constructing an amino acid residue interaction network, the structural connectivity and stability of the binding site are deeply analyzed. In summary, the present invention solves the problem that most of the existing methods only consider the static structural characteristics of molecules and ignore the impact of the dynamic changes of protein conformations on the formation of binding pockets. Brief Description of the Drawings

[0029] Figure 1 It is a flowchart of the method provided by the present invention;

[0030] Figure 2 It is the cumulative contribution rate graph of principal component analysis in Example 2;

[0031] Figure 3 It is the RMSF distribution graph of amino acid residues in Example 2;

[0032] Figure 4 It is the relative enrichment factor distribution graph of water molecules in Example 2;

[0033] Figure 5 It is the binding site network analysis graph in Example 2. Specific embodiments

[0034] To make the objectives, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention.

[0035] As Figure 1 shown, it is a flowchart of a method for identifying a binding pocket of a disease biomarker provided by the present invention, and this method includes the following steps:

[0036] S10. Construct a simulation system containing disease biomarker molecules, water molecules and ions, perform standardized pretreatment on the simulation system, and perform molecular dynamics simulation for 100 nanoseconds to 500 nanoseconds through molecular dynamics simulation software;

[0037] S20. Collect the conformational change trajectory of the disease biomarker during the molecular dynamics simulation, sample the conformational change trajectory at time intervals of 10 picoseconds to 50 picoseconds, and perform global normalization processing on the collected data;

[0038] S30. Perform conformational decomposition calculation on the conformations of the disease biomarker obtained by sampling, separate the stable conformations and the variable conformations, and select the conformation clusters with the stable conformation occurrence frequency in the top 20% to 30% for structural feature analysis;

[0039] S40. Calculate the root mean square fluctuation value and the secondary structure change rate of each amino acid residue in the conformation of the disease biomarker, establish an amino acid residue structure parameter matrix, and determine the region where the flexibility of the amino acid residue is greater than 0.2 nanometers and the secondary structure change rate is greater than 15% as the potential binding site;

[0040] S50. Perform solvent accessibility analysis on the potential binding site, calculate the change value of the exposed surface area and the electrostatic potential distribution, and select the region where the fluctuation range of the exposed surface area is between 0.5 square nanometers and 2.0 square nanometers and the electrostatic potential fluctuation is greater than 3 kcal / mol;

[0041] S60. Analyze the distribution density of solvent molecules in the region, identify the high-density regions where the residence time of water molecules exceeds 500 picoseconds and the distribution density is more than 3 times that of the bulk phase, and determine them as temporary binding pockets;

[0042] S70. Conduct principal component analysis and clustering analysis on the temporary binding pockets, extract the main conformational change patterns with a total variance contribution rate of more than 90%, and establish the characteristic vectors of the binding pockets;

[0043] S80. Construct the binding site map data including the interaction network of amino acid residues. Using amino acid residues as nodes and the interactions between residues as edges, calculate the binding free energy and conformational stability index of the temporary binding pockets. Select the sites with a binding free energy lower than -20 kJ / mol and a conformational stability index greater than 0.75 as stable binding pockets; and conduct a comprehensive performance analysis of the candidate binding pockets based on the pocket comprehensive performance evaluation equation, and select the pockets with a comprehensive evaluation index greater than 0.6 as stable binding pockets.

[0044] The following is a detailed description of the specific implementation manners of the above steps:

[0045] The specific implementation manner of step S10 is to construct a suitable molecular dynamics simulation system and perform preprocessing and simulation. First, place the disease biomarker molecule to be studied at the center of a three-dimensional cubic water box with a side length of 10 nanometers. Select the TIP3P water model to add water molecules to simulate the real water environment. According to the total charge of the system, achieve charge neutralization by adding an appropriate amount of chloride ions or sodium ions. At the same time, add sodium chloride with a concentration of 0.15 mol / L to simulate the ionic strength of the physiological environment, and construct a complete simulation system using periodic boundary conditions. Subsequently, perform a standardized preprocessing process on the constructed simulation system, including two stages: energy minimization and system equilibration. In the energy minimization stage, use the steepest descent method to perform 5000 steps of optimization, and set the convergence criterion as the maximum force being less than 1000 kJ / mol / nm. In the system equilibration stage, perform a 100-picosecond equilibration process under constant volume and constant temperature conditions. Use the Berendsen temperature coupling algorithm to control the temperature at 300 Kelvin, and set the time constant to 0.1 picosecond. Finally, use the molecular dynamics simulation software GROMACS and the CHARMM force field to perform molecular dynamics simulation. Set the time step to 2 femtoseconds, use the linear constraint solver algorithm to limit the bond lengths of all hydrogen bonds, use the particle mesh Ewald method to handle long-range electrostatic interactions, set the van der Waals interaction cutoff radius to 1.0 nanometer, and perform the simulation under the constant temperature and constant pressure ensemble. Maintain the pressure at 1 atmosphere through the Parrinello-Rahman method. The main purpose of this step is to construct a simulation system that can accurately reflect the behavior of the disease biomarker in the physiological environment and provide reliable kinetic data for subsequent analysis.

[0046] The specific implementation of step S20 is to collect and preprocess molecular dynamics simulation trajectory data. During the molecular dynamics simulation, the conformational information of the system is recorded every 20 picoseconds, including the coordinates, velocities, and energy data of all atoms. To eliminate the influence of the molecular cross-boundary movement caused by the periodic boundary conditions, a periodic image correction method is needed to process the trajectory. The protein backbone is structurally aligned by the least-squares superposition method to eliminate the influence of overall translation and rotation, so that subsequent analysis can focus on the conformational changes inside the protein. The collected data is globally normalized to unify the data of different physical quantities into the same numerical range, facilitating subsequent analysis and comparison. The purpose of this step is to obtain high-quality conformational change trajectory data, laying a foundation for subsequent conformational analysis.

[0047] The specific implementation of step S30 is to deeply analyze the conformations of the disease markers obtained by sampling. First, the principal component analysis method is used to reduce the dimension of the conformational ensemble. By calculating the eigenvectors and eigenvalues of the covariance matrix, the principal components with a cumulative contribution rate reaching 90% are selected as the main dimensions of the conformational space. Subsequently, the density peak clustering algorithm is used to classify the conformations. The clustering centers are determined through local density and distance criteria, and the conformations are divided into two categories: stable conformations and variable conformations. The conformational clusters with occurrence frequencies in the top 20% to 30% are selected from the stable conformations for structural feature analysis, including features such as secondary structure composition, key residue positions, and interaction networks. The purpose of this step is to identify the main conformational states of the disease markers, providing a conformational basis for subsequent binding site analysis.

[0048] The specific implementation of step S40 is to analyze the kinetic characteristics of each amino acid residue of the disease marker. For each amino acid residue, first, the root-mean-square fluctuation value of the position of its Ca atom is calculated, which reflects the flexibility of the residue in space. The DSSP algorithm is used to analyze the change of the secondary structure type of each amino acid residue during the simulation, and the time ratio of each amino acid residue maintaining the original secondary structure type under different conformations is statistically calculated. The inverse value is taken as the secondary structure change rate. The region where the flexibility of the amino acid residue is greater than 0.2 nanometers and the secondary structure change rate is greater than 15% is determined as the potential binding site. The purpose of this step is to identify potential binding sites on the protein surface through kinetic characteristics.

[0049] The specific implementation of step S50 is to perform solvent accessibility analysis on potential binding sites. The NACCESS program is used to calculate the solvent-accessible surface area of each amino acid residue in different conformations, and a rolling sphere calculation is performed with a probe radius of 1.4 Å. The APBS software is used to solve the Poisson-Boltzmann equation to calculate the molecular surface electrostatic potential distribution, with the internal dielectric constant set to 2 and the external dielectric constant set to 80. Regions with an exposed surface area fluctuation range of 0.5 square nanometers to 2.0 square nanometers and an electrostatic potential fluctuation greater than 3 kcal / mol are selected as candidate binding sites. The purpose of this step is to further screen potential binding sites from the perspectives of structure and electrostatic properties.

[0050] The specific implementation of step S60 is to analyze the distribution characteristics of solvent molecules at candidate binding sites. The simulation system is divided into cubic grids with a side length of 0.5 Å, the occurrence frequency of water molecule oxygen atoms in each grid is counted, and the enrichment factor relative to the bulk water density is calculated. The residence time of individual water molecules at specific positions is traced, and the kinetic characteristics of water molecules are evaluated through the continuous residence time and exchange frequency. High-density regions where the residence time of water molecules exceeds 500 picoseconds and the distribution density is more than 3 times that of the bulk phase are identified and determined as temporary binding pockets. The purpose of this step is to identify binding sites on the protein surface through the distribution characteristics of water molecules.

[0051] The specific implementation of step S70 is to conduct in-depth kinetic analysis of the temporary binding pocket. The principal component analysis method is used to extract the main conformational change patterns with a total variance contribution rate of more than 90%, and the eigenmotion characteristics of the binding pocket region are obtained. The conformational clustering analysis method is used to classify the conformations of the binding pocket at different times, and its conformational diversity and conversion characteristics are studied. A feature vector of the binding pocket containing geometric features, electrostatic properties, and kinetic parameters is established for subsequent pocket evaluation and screening. The purpose of this step is to obtain the kinetic characteristics of the binding pocket and provide a basis for evaluating its stability.

[0052] The specific implementation of step S80 is to construct and analyze the interaction network of binding sites. Network connections are established based on van der Waals contacts, hydrogen bonds, salt bridges, and hydrophobic interactions between amino acid residues, and the interaction strength determines the edge weight according to the contact frequency and energy contribution. Graph theory methods are used to calculate topological parameters such as the degree distribution, clustering coefficient, and centrality of the network. The binding free energy of the binding pocket is calculated by the molecular mechanics Poisson-Boltzmann surface area method, considering the contributions of solvation free energy and conformational entropy. Sites with a binding free energy lower than -20 kJ / mol and a conformational stability index greater than 0.75 are selected as stable binding pockets. Based on the comprehensive performance evaluation equation of the pocket, the comprehensive evaluation index is calculated, and pockets with a comprehensive evaluation index greater than 0.6 are selected as the final stable binding pockets. The purpose of this step is to evaluate the stability of the binding pocket from the perspectives of thermodynamics and structural network.

[0053] The following is a detailed description of the calculation process involved in the present invention.

[0054] 1. Principal Component Analysis (PCA) equation for conformational decomposition calculation:

[0055] The calculation of the principal component analysis covariance matrix is expressed as follows:

[0056]

[0057] Where C ij is the covariance matrix element; N is the number of conformations; x ki is the coordinate of the i-th atom in the k-th conformation; μ i is the average value of the coordinates of the i-th atom.

[0058] The eigenvalue decomposition is expressed as follows:

[0059] Cv i = λ i v i ;

[0060] Where λ i is the eigenvalue; v i is the corresponding eigenvector.

[0061] 2. The calculation of the Root Mean Square Fluctuation (RMSF) is expressed as follows:

[0062]

[0063] Where RMSF i is the root mean square fluctuation value of the i-th residue; T is the number of sampling frames; r i (t) is the position of the i-th residue in the t-th frame; is the position of the i-th residue in the reference conformation.

[0064] 3. The calculation of the secondary structure change rate is expressed as follows:

[0065]

[0066] Where SSC i is the secondary structure change rate of the i-th residue; is the number of frames maintaining the original secondary structure type; T total is the total number of frames.

[0067] 4. Solvent accessibility analysis equation:

[0068] The calculation of the solvent accessible surface area is expressed as follows:

[0069]

[0070] Where, SASA i is the solvent accessible surface area of the i-th residue; r p is the probe radius; M is the number of surface sampling points; Δθ j is the solid angle increment of the j-th sampling point.

[0071] 5. Electrostatic potential distribution calculation (Poisson - Boltzmann equation):

[0072]

[0073] Where, ∈(r) is the dielectric constant; φ(r) is the electrostatic potential; ρ f (r) is the fixed charge density; z i is the ionic valence; e is the elementary charge; is the ionic concentration; k B is the Boltzmann constant; T is the temperature.

[0074] 6. Water molecule distribution density analysis:

[0075] The relative enrichment factor is calculated as follows:

[0076]

[0077] Where, REF ijk is the relative enrichment factor at the grid point (i, j, k); ρ ijk is the water molecule density at this point; ρ bulk is the bulk water density.

[0078] 7. Residue interaction network construction:

[0079] The interaction strength is calculated as follows:

[0080] W ij = αf ij + βE ij ;

[0081] Where, W ij is the interaction weight between residues i and j; f ij is the contact frequency; E ij is the energy contribution; α, β are weight coefficients.

[0082] 8. Binding free energy calculation (MM / PBSA method):

[0083] The total binding free energy is expressed as follows:

[0084] ΔG bind = ΔE MM + Δ Gsol - TΔS;

[0085] where ΔG bind is the binding free energy; ΔE MM is the change in molecular mechanics energy; ΔG sol is the change in solvation free energy; T is the temperature; ΔS is the entropy change.

[0086] The calculation of the change in molecular mechanics energy is expressed as follows:

[0087] ΔE MM = ΔE ele + ΔE vdw + ΔE int ;

[0088] where ΔE ele is the change in electrostatic interaction energy; ΔE vdw is the change in van der Waals interaction energy; ΔE int is the change in internal energy.

[0089] The calculation of the change in solvation free energy is expressed as follows:

[0090] ΔG sol = ΔG polar + ΔG nonpolar ;

[0091] where ΔG polar is the change in polar solvation free energy; ΔG nonpolar is the change in non-polar solvation free energy.

[0092] 9. The equation for evaluating the comprehensive performance of the pocket is expressed as follows:

[0093]

[0094] Parameter description: In the equation, P is the comprehensive performance score of the pocket, with a value range of 0 - 1. When P is greater than 0.6, it can be determined as a stable binding pocket. The higher the score, the more stable and reliable the binding pocket. ΔG is the binding free energy, with the unit of kJ / mol, which is calculated by the MM / PBSA method after molecular dynamics simulation, considering the comprehensive contributions of molecular mechanics energy, solvation free energy, and conformational entropy. ΔG0 is its standard reference value. RMSF is the root mean square fluctuation value of key residues, with the unit of nm, which is calculated based on the molecular dynamics simulation trajectory to obtain the standard deviation of the Cα atom positions of residues, used to reflect the conformational flexibility of the binding pocket region. RMSF0 is its standard reference value. REF is the relative enrichment factor of water molecules, dimensionless, obtained by statistically analyzing the distribution density of water molecules in the pocket region, and using the ratio with the bulk water density to characterize solvent accessibility. REF0 is its standard reference value. CN is the sequence conservation score, dimensionless, calculated based on sequence alignment and evolutionary analysis, reflecting the degree of conservation of the pocket region during evolution. CN0 is its standard reference value. α, β, γ, and δ are the weight coefficients of each parameter, and their sum is equal to 1. The specific weight values are determined according to the contribution degree of each parameter to the stability of the binding pocket. This equation comprehensively considers thermodynamic stability, conformational dynamic characteristics, solvent effects, and evolutionary conservation, realizing a comprehensive evaluation of the binding pocket of disease markers. Each parameter can be directly obtained through molecular dynamics simulation and structural analysis, ensuring the reliability and operability of the evaluation results.

[0095] Explanation of the construction principle of each equation:

[0096] 1) The principal component analysis equation is calculated based on the covariance matrix, and extracts the main conformational change patterns through eigenvalue decomposition, which can effectively reduce the dimension and retain key information;

[0097] 2) The RMSF calculation considers the fluctuations of residue positions over time, reflecting local structural flexibility;

[0098] 3) The secondary structure change rate is in the form of a ratio, intuitively reflecting structural stability;

[0099] 4) The solvent accessibility analysis uses a geometric method, simulating water molecule contact through a probe sphere;

[0100] 5) The Poisson - Boltzmann equation considers the physical nature of the electrostatic potential distribution, including the dielectric constant and ionic effects;

[0101] 6) The water molecule distribution density analysis uses the relative enrichment factor, which is convenient for comparison with the bulk phase;

[0102] 7) The residue interaction network comprehensively considers the contact frequency and energy contribution through a weighted method;

[0103] 8) The binding free energy calculation is based on a thermodynamic cycle and takes into account the molecular mechanics energy and solvation effects.

[0104] Specifically, the principle of the present invention is based on the basic principles of protein molecular dynamics and structural biology. First, by constructing a complete simulation system containing disease marker molecules, water molecules, and ions, the periodic boundary conditions are used to ensure the continuity of the system. The TIP3P water model and appropriate ion concentration are used to simulate the physiological environment, providing a reliable physicochemical basis for subsequent analysis.

[0105] During the molecular dynamics simulation, the CHARMM force field is used to describe the interactions between various atoms in the system, including bonded and non-bonded interactions. By setting a reasonable time step and using the LINCS algorithm to restrict the vibrations containing hydrogen bonds, the stability of the simulation is ensured. The PME method is used to handle long-range electrostatic interactions, ensuring the accuracy of the calculation. The simulation is carried out in the NPT ensemble, and the stability of the system under physiological conditions is maintained through the coupling algorithms of temperature and pressure.

[0106] Conformational analysis is performed using principal component analysis for dimensionality reduction. This is based on the principle of eigenvalue decomposition in mathematics and can extract the main motion modes of the system. The conformational classification is carried out using the density peak clustering algorithm, which is based on the principle of local density and distance criteria and can effectively identify the main conformational clusters in the conformational space.

[0107] During the binding site identification process, the root mean square fluctuation value and the secondary structure change rate of amino acid residues are calculated to reflect the flexibility of local regions of the protein. The solvent accessibility analysis is based on the rolling ball algorithm, and the accessibility of potential binding sites is evaluated by calculating the exposed surface area. The calculation of the electrostatic potential distribution uses the Poisson-Boltzmann equation, taking into account the differences in the dielectric environment inside and outside the molecule.

[0108] The water molecule distribution density analysis is based on the principles of statistical mechanics. By calculating the enrichment degree of local water molecule density relative to the bulk phase and combining the residence time of water molecules, the hydration sites with functional significance are identified. These sites are often potential ligand binding sites because they have specific physicochemical properties that are conducive to forming stable intermolecular interactions.

[0109] The construction of the amino acid residue interaction network is based on graph theory principles. By analyzing various interactions between residues, a topological structure describing the local environment of the binding site is established. The binding free energy is calculated using the MM / PBSA method, comprehensively considering the contributions of molecular mechanics energy, electrostatic interaction, non-polar solvation free energy, and conformational entropy.

[0110] The application of the cross - validation method is based on statistical principles. By randomly splitting the data set and performing multiple rounds of validation, the statistical significance and reliability of the model prediction results are ensured.

[0111] A specific embodiment 1 of the present invention is provided below. The specific implementation of each step in this embodiment 1 is described in detail as follows: First, for the specific implementation of step S10, a simulation system containing disease biomarker molecules, water molecules, and ions needs to be constructed. This process can be described using the following formula:

[0112] Construct the simulation system:

[0113] r protein =(x protein ,y protein ,z protein );

[0114] r water =(x water ,y water ,z water );

[0115] r ion =(x ion ,y ion ,z ion );

[0116] r total ={r protein ,r water ,r ion};

[0117] V box =10nm×10nm×10nm;

[0118] C NaCl =0.15mol / L;

[0119] Among them, r protein , r water , r ion respectively represent the coordinate vectors of disease biomarker molecules, water molecules, and ions, r total is the total coordinate set of the simulation system, V box is the volume of the simulation box, and c NaCl is the added sodium chloride concentration.

[0120] For the standardized pre - processing, the energy minimization and system equilibration processes can be described using the following formula:

[0121] Energy minimization:

[0122]

[0123] ||F i ||<1000 kJ / (mol·nm);

[0124] System equilibrium:

[0125] T = 300 K;

[0126] p = 1 atm;

[0127]

[0128] where F i is the force on the i-th atom, U is the total potential energy function, v i is the velocity of the i-th atom, Δt is the time step, T0 and p0 are the target temperature and pressure respectively, τ T and τ p are the temperature and pressure coupling time constants.

[0129] For molecular dynamics simulations, the following equations can be used to describe the relevant parameter settings:

[0130] Molecular dynamics simulation:

[0131] Δt = 2 fs;

[0132] r cutoff = 1.0 nm;

[0133]

[0134] where Δt is the time step, r cutoff is the cutoff radius for van der Waals interactions, m i is the mass of the i-th atom.

[0135] Step S20 involves sampling the conformational change trajectory and global normalization, which can be described by the following equations:

[0136] Sampling of conformational change trajectory:

[0137] r protein (t) = {r i (t)|i = 1, 2, …, N atom};

[0138] t = n·Δt s , n = 0, 1, 2, …, N sample ;

[0139] Δt s = 20 ps;

[0140] Global normalization:

[0141]

[0142] where, r protein (t) is the coordinate of the disease marker at time t, N atom is the number of atoms, Δt s is the sampling time interval, N sample is the number of sampling frames, <·> represents taking the average value.

[0143] The conformational decomposition calculation involved in step S30 includes principal component analysis and clustering analysis, and can be described by the following formula:

[0144] Principal component analysis:

[0145]

[0146] Cv i = λ i v i ;

[0147]

[0148] where, C is the covariance matrix, x k is the coordinate vector of the k-th conformation, μ is the average value of the coordinate vectors, v i and λ i are the i-th eigenvector and eigenvalue respectively.

[0149] Clustering analysis:

[0150]

[0151] ρ i > ρ c , d i > d c ;

[0152] where, ρ i is the local density of the i-th conformation, d ij is the distance between the i-th and j-th conformations, δ is the kernel function bandwidth parameter, ρ c and d c are the thresholds of density and distance.

[0153] The calculation of the amino acid residue structure parameters involved in step S40 can be described by the following formula:

[0154] RMSF calculation:

[0155]

[0156] Calculation of the secondary structure change rate:

[0157]

[0158] Among them, RMSF i is the root mean square fluctuation value of the i-th residue, and r i (t) is the coordinate of the i-th residue in the t-th frame, is the coordinate of the i-th residue in the reference conformation, T is the total number of sampling frames, and SSC i is the secondary structure change rate of the i-th residue, is the number of frames maintaining the original secondary structure type, and T total is the total number of frames.

[0159] The solvent accessibility analysis involved in step S50 can be described by the following formula:

[0160] Solvent accessible surface area calculation:

[0161]

[0162] Electrostatic potential distribution calculation:

[0163]

[0164] Among them, SASA i is the solvent accessible surface area of the i-th residue, and r p is the probe radius, M is the number of surface sampling points, and Δθ j is the solid angle increment of the j-th sampling point, ∈(r) is the dielectric constant, φ(r) is the electrostatic potential, and ρ f (r) is the fixed charge density, and z i is the ionic valence state, e is the elementary charge, is the ionic concentration, and k B is the Boltzmann constant, and T is the temperature.

[0165] The water molecule distribution density analysis involved in step S60 can be described by the following formula:

[0166] Relative enrichment factor calculation:

[0167]

[0168] Among them, REF ijk is the relative enrichment factor at the grid point (i, j, k), and ρ ijk is the water molecule density at this point, and ρ bulk is the bulk water density.

[0169] The principal component analysis and clustering analysis involved in step S70 are similar to step S30 and will not be elaborated here.

[0170] The construction of the amino acid residue interaction network and the calculation of the binding free energy involved in step S80 can be described by the following formula:

[0171] Construction of residue interaction network:

[0172] W ij = αf ij + βE ij ;

[0173] Calculation of binding free energy:

[0174] ΔG bind = ΔE MM + ΔG sol - TΔS;

[0175] ΔE MM = ΔE ele + ΔE vdw + ΔE int ;

[0176] ΔG sol = ΔG polar + ΔG nonpolar ;

[0177] Where, W ij is the interaction weight between residues i and j, f ij is the contact frequency, E ij is the energy contribution, α and β are weight coefficients, ΔG bind is the binding free energy, ΔE MM is the change in molecular mechanics energy, ΔG sol is the change in solvation free energy, T is the temperature, ΔS is the entropy change, ΔE ele is the change in electrostatic interaction energy, ΔE vdw is the change in van der Waals interaction energy, ΔE int is the change in internal energy, ΔG polar is the change in polar solvation free energy, ΔG nonpolar is the change in non-polar solvation free energy.

[0178] Next, based on the pocket comprehensive performance evaluation equation, the comprehensive performance analysis of the candidate binding pockets is carried out:

[0179] First, perform MM / PBSA calculation on the temporary binding pocket to obtain the binding free energy ΔG. When calculating, consider the contributions of molecular mechanics energy, polar and non-polar solvation free energies, and conformational entropy, and select the conformation after system equilibrium for calculation;

[0180] Second, calculate the root mean square fluctuation value RMSF of the key residues of the pocket based on the molecular dynamics simulation trajectory. Use the least squares superposition method to perform structural alignment on the protein backbone to eliminate the influence of overall motion, and statistically analyze the standard deviation of the Cα atom position;

[0181] Next, analyze the water molecule distribution characteristics in the pocket region to obtain the relative enrichment factor REF. Divide the simulation system into cubic grids with a side length of 0.5 Å, count the occurrence frequency of water molecule oxygen atoms, and calculate the enrichment degree relative to the bulk water density;

[0182] Then, calculate the conservation score CN of the pocket region through sequence alignment and evolutionary analysis, considering the similarity of the physicochemical properties of amino acid residues, the conservation of secondary structure, and the conservation of surface exposure;

[0183] Finally, substitute the above four parameters into the pocket comprehensive performance evaluation equation:

[0184] Among them, ΔG0, RMSF0, REF0, and CN0 are the standard reference values of each parameter, and α, β, γ, and δ are weight coefficients and their sum is 1; the calculated comprehensive performance score P ranges from 0 to 1. Select the sites with P greater than 0.6 as stable binding pockets; for the selected stable binding pockets, they are used as the final recognition results.

[0185] To better understand and implement the present invention, the following provides Example 2 of a specific application scenario of the present invention: A research team used the disease biomarker binding site recognition method proposed by the present invention to analyze an important lung cancer biomarker, EGFR (epidermal growth factor receptor). EGFR is a transmembrane receptor protein, and its extracellular domain can bind specific growth factors, thereby activating downstream signal pathways and promoting the proliferation and differentiation of tumor cells. Accurately identifying the key binding sites in the extracellular domain of the EGFR protein is of great significance for deeply understanding its mechanism of action in the occurrence and development of lung cancer and developing targeted drugs.

[0186] First, the research team constructed a simulation system containing the extracellular domain of EGFR, water molecules, and ions. Specifically, place the atomic coordinates of the extracellular domain of EGFR (PDB ID: 1M17) at the center of a cubic water box with a side length of 10 nm, add water molecules using the TIP3P water model, neutralize the charge with chloride ions according to the total charge of the EGFR domain, and add 0.15 mol / L sodium chloride to simulate the physiological environment. The final constructed simulation system contains 1 extracellular domain of EGFR, 22493 water molecules, and 164 chloride ions, with a total number of atoms of approximately 72000.

[0187] The simulation system was subjected to standardized preprocessing. First, the steepest descent method was used for 5000 steps of energy minimization optimization, and the convergence criterion was set as the maximum force being less than 1000 kJ / mol / nm. Then, under the constant volume and constant temperature conditions, that is, in the NPT ensemble, the system was equilibrated for 100 picoseconds, the temperature was controlled at 300 Kelvin, the Berendsen temperature coupling algorithm was used, and the time constant was set to 0.1 picosecond. The pressure was maintained at 1 atmosphere through the Parrinello-Rahman method.

[0188] Next, the research team used the GROMACS software package and the CHARMM36 force field to perform a 300-nanosecond molecular dynamics simulation on this simulation system. The time step was set to 2 femtoseconds, the LINCS algorithm was used to constrain the bond lengths of all hydrogen bonds, the PME method was used to handle long-range electrostatic interactions, and the cutoff radius for van der Waals interactions was set to 1.0 nanometer. Throughout the simulation process, the temperature and pressure were maintained at 300 Kelvin and 1 atmosphere.

[0189] During the molecular dynamics simulation, the research team recorded the conformational information of the EGFR extracellular domain every 20 picoseconds, including atomic coordinates, velocities, and energy data. To eliminate the influence of overall translation and rotation, the least squares superposition method was used to align the EGFR backbone structure. Finally, a total of 15000 frames of conformational trajectory data were obtained. Global normalization was performed on these data to eliminate the systematic errors introduced by different scales.

[0190] Next, the research team performed principal component analysis and clustering analysis on the sampled EGFR extracellular domain conformations. First, the covariance matrix was calculated and eigenvalue decomposition was performed on it to extract the first 20 eigenvectors, with the cumulative contribution rate reaching 98.4%. The density peak-based clustering method was used to classify the conformations, and after multiple attempts, 11 representative conformational clusters were finally determined. Among them, 5 conformational clusters with the top 30% occurrence frequencies were considered relatively stable conformations. As Figure 2 shown, it presents the change in the cumulative contribution rate of the principal component analysis of the EGFR extracellular domain conformation. The horizontal axis represents the number of principal components, and the vertical axis represents the cumulative contribution rate (%). The blue dotted line in the figure shows the change trend of the cumulative contribution rate of the first 20 principal components, and the final cumulative contribution rate value (98.4%) is marked at the 20th principal component. This figure intuitively shows that most of the conformational changes can be explained by using a small number of principal components.

[0191] Then, the research team calculated the RMSF and SSC values of each amino acid residue in the extracellular domain of EGFR. As shown in Table 1, the amino acid residues with RMSF values greater than 0.2 nm and SSC values greater than 15% are mainly distributed in the connecting region between domains I and II, and some surface regions of domain III. These regions may be involved in the binding process of EGFR and are potential binding sites.

[0192] Table 1 RMSF and SSC values of key amino acid residues in the extracellular domain of EGFR

[0193] Residue number Amino acid type RMSF (nm) SSC (%) 222 Asp 0.23 16.7 223 Ser 0.25 17.4 224 Tyr 0.27 18.2 311 Glu 0.21 15.8 312 Arg 0.24 16.5 313 Leu 0.22 15.4 445 Asp 0.26 17.9 446 Gly 0.28 18.7 447 Phe 0.25 17.1

[0194] As Figure 3 shown, the root mean square fluctuation (RMSF) distribution of each amino acid residue in the extracellular domain of EGFR is presented. The horizontal axis represents the residue number (200 - 500), and the vertical axis represents the RMSF value (nm). The red solid line in the figure indicates the change of RMSF with the residue number, and the positions of three important regions (222 - 224, 311 - 313, 445 - 447) are specifically marked. These regions show relatively high conformational flexibility.

[0195] To further evaluate the properties of these potential binding sites, the research team conducted solvent accessibility analysis and electrostatic potential distribution calculation on them. Using the NACCESS program, the solvent accessible surface area (SASA) of each residue in the extracellular domain of EGFR in different conformations was calculated. It was found that the SASA fluctuation range of residues 222 - 224 and 445 - 447 is 0.6 - 1.8 square nanometers, while that of 311 - 313 is 0.5 - 1.5 square nanometers. At the same time, using the APBS software, the electrostatic potential distribution on the surface of the extracellular domain of EGFR was calculated, and it was found that the electrostatic potential fluctuation in these regions is greater than 3 kcal / mol. The above results indicate that these regions have strong affinity and solubility and are very likely to be involved in the binding process of EGFR.

[0196] On this basis, the research team further analyzed the distribution characteristics of water molecules near these potential binding sites. The simulation system was divided into cubic grids with a size of 0.5 Å, and the occurrence frequency of water molecule oxygen atoms in each grid was statistically analyzed to calculate the enrichment factor (REF) relative to the bulk water density. As shown in Table 2, in the regions where residues 222 - 224 and 445 - 447 are located, there are high-density water molecule regions with REF greater than 3, and the average residence time of water molecules exceeds 500 picoseconds. In the region where 311 - 313 is located, although there are also high-density water molecule regions, the average residence time is relatively short, only about 300 picoseconds. This indicates that the first two regions may form relatively stable temporary binding pockets, while the latter has relatively weak affinity.

[0197] Table 2 Distribution characteristics of water molecules near potential binding sites in the extracellular domain of EGFR

[0198] Site Highest REF Average residence time (ps) 222-224 3.4 540 311-313 2.8 290 445-447 3.7 610

[0199] As Figure 4 shown, the relative enrichment factor (REF) distribution of water molecules on the surface of the EGFR extracellular domain was presented in the form of a heat map. The darker the color, the higher the REF value, indicating a higher degree of water molecule enrichment. Three main water molecule enrichment regions (222 - 224, 311 - 313, 445 - 447) were marked on the figure, and the highest REF value of each region was noted.

[0200] For the above three transient binding pockets, the research team further used principal component analysis and clustering analysis to extract their key conformational change characteristics. For the 222 - 224 region, the first 3 principal components explained 91.3% of the total variance; for the 311 - 313 region, the first 3 principal components explained 88.4% of the total variance; for the 445 - 447 region, the first 3 principal components explained 93.2% of the total variance. This indicates that the conformational changes in these three regions can be described by relatively few principal components, suggesting that they have relatively stable structural characteristics.

[0201] Subsequently, the research team constructed a binding site map containing the amino acid residue interaction network based on these transient binding pockets. Specifically, first, weighted network connections were established according to van der Waals contacts, hydrogen bonds, salt bridges, and hydrophobic interactions between residues. As shown in Table 3, the network in the 222 - 224 region has a relatively high degree distribution and clustering coefficient, indicating relatively strong interactions between these residues; the network connectivity in the 311 - 313 region is relatively weak; the network in the 445 - 447 region is between the two. In addition, the binding free energy and conformational stability index of these three regions were calculated according to the MM / PBSA method. The results showed that the binding free energy of the 222 - 224 region was -23.7 kJ / mol, and the conformational stability index was 0.82, far superior to the other two regions.

[0202] Table 3 Network topological parameters and thermodynamic indices of transient binding pockets in the EGFR extracellular domain

[0203] Site Degree distribution Clustering coefficient Binding free energy (kJ / mol) Conformational stability index 222-224 5.8 0.71 -23.7 0.82 311-313 4.2 0.48 -18.5 0.64 445-447 5.1 0.59 -21.3 0.73

[0204] Taking into account the analysis results of the above aspects, the research team finally determined that the main binding site of the EGFR extracellular domain is located in the 222-224 region. This region not only has relatively high local structural flexibility and stability, but also its surface properties are highly relevant to the binding process. It also forms a relatively stable temporary binding pocket, and there is a tight interaction network among the residues, and the binding free energy is also relatively superior. Although the 311-313 region and the 445-447 region also show a certain binding tendency, their overall performance is significantly weaker than that of the 222-224 region.

[0205] As Figure 5 shown, the amino acid residue interaction networks of three potential binding site regions are presented. Nodes of different colors represent amino acid residues in different regions, and the connecting lines indicate the interactions between residues. The network topology parameters (degree distribution and clustering coefficient) of each region are also marked in the figure, visually showing that the 222-224 region has the tightest interaction network.

[0206] To further verify the reliability of the recognition results, the research team evaluated the obtained results by means of cross-validation. The conformational dataset of the EGFR extracellular domain was randomly divided into a training set (75%) and a test set (25%). A binding site feature model was constructed using the training set and evaluated on the test set. After 5 rounds of cross-validation, the average accuracy of the model on the test set was 91.3%, the average precision was 87.5%, and the average recall was 89.2%. The high values of these indicators indicate that the main binding site of the EGFR extracellular domain identified by the method of the present invention has strong reliability and biological significance.

[0207] It should be noted that the detailed explanations of the variables involved in the present invention are shown in Table 4 below.

[0208] Table 4 Variable Explanation Table

[0209]

[0210]

[0211] The above is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention can easily think of changes or substitutions, which should all be covered by the protection scope of the present invention.

Claims

1. A method for identifying a disease marker binding pocket, characterized in that, Including: Construct a simulation system containing disease biomarker molecules, water molecules and ions, perform standardized preprocessing on the simulation system, and conduct molecular dynamics simulations; collect the conformational change trajectories of disease biomarkers during the molecular dynamics simulations, sample the conformational change trajectories, and perform global normalization processing on the collected data; Perform conformational decomposition calculations on the sampled disease biomarker conformations to separate stable conformations and fluctuating conformations; calculate the root mean square fluctuations and secondary structure change rates of each amino acid residue in the disease biomarker conformations, establish an amino acid residue structure parameter matrix, and determine the region where the flexibility of amino acid residues is greater than 0.2 nanometers and the secondary structure change rate is greater than 15% as potential binding sites; Conduct solvent accessibility analysis on the potential binding sites to identify high-density water molecule regions and determine them as temporary binding pockets; perform principal component analysis and clustering analysis on the temporary binding pockets to establish binding pocket eigenvectors; construct binding site map data, calculate the binding free energy and conformational stability index of the temporary binding pockets, and select stable binding pockets.

2. The method for identifying a disease marker binding pocket according to claim 1, wherein To construct a simulation system containing disease biomarker molecules, water molecules and ions, place the disease biomarker molecules at the center of a cubic water box with a side length of 10 nanometers, add water molecules using the TIP3P water model, neutralize the charge using chloride ions or sodium ions according to the total charge of the system, add 0.15 moles per liter of sodium chloride to simulate the physiological environment, and use periodic boundary conditions to construct a complete simulation system.

3. The method for identifying a disease marker binding pocket according to claim 2, wherein To perform standardized preprocessing on the simulation system, first perform energy minimization on the system, optimize it for 5000 steps using the steepest descent method, set the convergence criterion to a maximum force less than 1000 kilojoules per mole per nanometer, perform a 100 picosecond system equilibration under constant volume and constant temperature conditions, control the temperature at 300 Kelvin, and use the Berendsen temperature coupling algorithm with a time constant set to 0.1 picosecond.

4. The method for identifying a disease biomarker binding pocket according to claim 3, wherein, To conduct molecular dynamics simulations, use the GROMACS software package, perform molecular dynamics simulations using the CHARMM force field, set the time step to 2 femtoseconds, use the LINCS algorithm to constrain the bond lengths of all hydrogen bonds, use the PME method to handle long-range electrostatic interactions, set the van der Waals interaction cutoff radius to 1.0 nanometer, and conduct simulations in the NPT ensemble, maintaining the pressure at 1 atmosphere using the Parrinello-Rahman method.

5. The method for identifying a disease marker binding pocket according to claim 4, wherein, To collect the conformational change trajectories of disease biomarkers during the molecular dynamics simulations, record the system conformational information, including atomic coordinates, velocities, and energy data, every 20 picoseconds, use periodic image correction to remove the influence of molecular cross-boundary motion, and perform structural alignment of the protein backbone using the least squares superposition method to eliminate the influence of overall translation and rotation.

6. The method for identifying a disease biomarker binding pocket according to claim 5, wherein To perform conformational decomposition calculations on the sampled disease biomarker conformations, use the principal component analysis method to reduce the dimensionality of the conformational clusters, calculate the eigenvectors and eigenvalues of the covariance matrix, select the principal components with a cumulative contribution rate reaching 90% as the main dimensions of the conformational space, use the density peak clustering algorithm to classify the conformations, and determine the clustering centers through local density and distance criteria.

7. The method for identifying a disease marker binding pocket according to claim 6, wherein Calculate the root-mean-square fluctuation value and the secondary structure change rate of each amino acid residue in the conformation of the disease biomarker. The standard deviation of the Cα atom position of each amino acid residue is calculated for the aligned trajectory. The DSSP algorithm is used to analyze the change of the secondary structure type of each amino acid residue during the simulation, and the proportion of the time that each amino acid residue maintains its original secondary structure type in different conformations is statistically analyzed. The negative value of it is taken as the secondary structure change rate.

8. A method for identifying a disease biomarker binding pocket according to claim 7, characterized in that, Perform solvent accessibility analysis on potential binding sites. Calculate the solvent accessible surface area of each amino acid residue in different conformations using the NACCESS program, perform rolling sphere calculations with a probe radius of 1.4 Å, and use the APBS software to solve the Poisson-Boltzmann equation to calculate the molecular surface electrostatic potential distribution, with the internal dielectric constant set to 2 and the external dielectric constant set to 80.

9. A method for identifying a disease marker binding pocket according to claim 8, characterized in that, Before identifying the high-density region of water molecules, it is necessary to perform solvent molecule distribution density analysis in the region. The simulation system is divided into cubic grids with a side length of 0.5 Å, the occurrence frequency of water molecule oxygen atoms in each grid is statistically analyzed, the enrichment factor relative to the bulk water density is calculated, the residence time of a single water molecule at a specific position is traced, and the dynamic characteristics of water molecules are evaluated through the continuous residence time and exchange frequency.

10. A method for identifying a disease marker binding pocket according to claim 9, characterized in that, Construct the binding site map data, that is, construct the binding site map data containing the amino acid residue interaction network. Specifically, network connections are established based on the van der Waals contacts, hydrogen bonds, salt bridges, and hydrophobic interactions between amino acid residues. The interaction strength determines the edge weight according to the contact frequency and energy contribution. The NetworkX library is used to calculate the degree distribution, clustering coefficient, and centrality topological parameters of the network.

Images