CN113539381B - Molecular dynamics result analysis method based on residue interaction and PEN - Google Patents

Abstract

The invention relates to a molecular dynamics result analysis method based on residue interaction and PEN, which belongs to the technical field of pharmacokinetics, and comprises the following steps: 1) Constructing a residue interaction energy matrix; 2) Constructing and analyzing a protein energy network structure diagram; 3) And classifying and training the graph convolution neural network. According to the method, the medicine molecular dynamics analysis is more efficient and accurate through the interaction energy matrix, the protein energy network structure diagram, the molecular dynamics analysis and the graph convolution neural network classification training.

Description

Molecular dynamics result analysis method based on residue interaction and PEN

Technical Field

The invention belongs to the technical field of pharmacokinetics, and particularly relates to a method for analyzing the results of pharmacokinetics based on protein residue interaction and protein energy network analysis.

Background

The most important of pharmaceutical chemistry is to find the lead compound, one of two main supports of the source of the lead compound is a virtual screening technology, and as the final step of the serial technology of the virtual screening technology, the molecular dynamics plays an important role in selecting effective molecules for synthesis, so that the positive success rate of the target enzyme lead compound can be improved; in addition, the molecular dynamics can also be used for researching the known compounds, so as to obtain a space for further modification.

Molecular dynamics simulation has become a powerful tool and an indispensable research means for researching molecular conformational changes and functional analyses, and is widely used in the field of drug design. The method can predict, guide and explain experiments to a great extent, and the combination of computational simulation and experiments has become the main research means at present.

Molecular dynamics modeling yields a large amount of information related to protein dynamics. If properly analyzed, this information may bring new insights into the function of the protein. The main output of molecular dynamics simulation is the trajectory, which typically includes thousands of conformations of the biomolecular system.

In studies involving molecular dynamics simulation, the most complex part is the analysis of these data, not the actual simulation itself. Thus, obtaining an understandable conclusion from such multidimensional data can be a time-consuming task requiring the use of extensive and complex analytical methods and procedures.

The output of the molecular dynamics simulation analysis process is the use of individual amino acids as building blocks to determine the dynamics of the protein and to determine its relative (relevant) importance in terms of function and activity. However, it is far from straightforward to detect the effect of a single residue from a large amount of analog data. To obtain information about residue levels, the pair-wise non-bond interaction energy between amino acid residues can be calculated by using a single conformation or a set of conformations (e.g., conformations obtained from molecular dynamics simulation).

Protein Structure Network (PSN) refers to the use of a network in protein structure. PEN (protein energy network) can be considered a special type of "protein structure network" (PSN) in which interactions between individual residues can be used to construct IEMs (interaction energy matrices) that represent the "intensities" of edges between nodes (residues) in a network structure. These intensity values may be normalized to further represent "weights" or "information transfer costs" used from one node to another in the network analysis task, such as shortest path identification, etc.

The existing molecular dynamics result analysis method has a great improvement on speed, but has a plurality of defects.

1. In the analysis of molecular dynamics results, a method for calculating the binding free energy of a compound conjugate is generally adopted, but the method ignores the structural characteristics of the compound, can not comprehensively reflect the integral binding degree of a target point and the compound, and has certain unilaterality.

2. The main output of molecular dynamics simulation is the trajectory, typically comprising thousands of conformations of the biomolecular system, the most complex part of the study involving MD simulation being the analysis of the data rather than the simulation itself. Thus, obtaining an understandable conclusion from such multidimensional data requires the use of extensive and complex analytical methods and procedures, a time-consuming task.

Disclosure of Invention

Aiming at the defects existing in the prior art, the invention provides a molecular dynamics result analysis method based on residue interaction and PEN, which enables the analysis of the pharmacokinetics to be more efficient and accurate through Interaction Energy Matrix (IEM) and protein energy network structure diagram (PEN) and classification training of a graph convolution neural network.

In order to solve the technical problems, the invention adopts the following technical scheme:

a method of analysis of molecular dynamics results based on residue interactions and PEN, the method comprising the steps of:

1) Constructing a residue interaction energy matrix;

2) Constructing and analyzing a protein energy network structure diagram;

3) And classifying and training the graph convolution neural network.

The specific operation of the step 1) is as follows:

extracting residue-to-residue interaction energy from the simulated trajectory using a tool while calculating correlations between them, wherein interactions between individual residues can be used to construct an interaction energy matrix representing the "strength" of edges between residues in the network structure, the interaction energy matrix being abbreviated IEM;

further, the interaction energy between individual residues is the sum of the non-bond interaction energies defined in the force field; non-bond interactions can be considered with van der Waals forces and electrostatic terms:

here, i and j each represent two residues, E _ij Is an array that calculates the interaction energy between residues i and j in all trajectory frames,for electrostatic interaction>Is van der waals interactions.

Further, interactions between individual residues can be used to construct IEMs (interaction energy matrices) that represent the "strength" of edges between nodes (residues) in a network structure.

S _ij ＝|E _ij ，N _ij |

S _ij ＝|E _ij ，N _ij I, here S _ij Representing an interaction energy matrix, E _ij Is the interaction energy between residues i and j in all trace frames, N _ij Is the number of residues in the trace frame.

The specific operation of the step 2) is as follows: the IEM of the "intensity" of edges between residues in the network structure is normalized to further represent the "weight" or "information transfer cost" used from one residue to another in the network analysis task, and the resulting Interaction Energy Matrix (IEM) is used to construct and analyze a protein energy network, abbreviated PEN, based on residue-based network metrics (e.g., degree, middle-center, and proximity-center) and the shortest path between selected residues in the structure.

Further, the "weight" attribute of the edge is determined using the following formula: wherein ω represents the edge weight between residues i and j, and if residues i and j are covalent bonds, the weight between them is 0.99, noThen the weight between them is x _ij ，x _ij Representing the average interaction energy between residues i and j.

The specific operation of the step 3) is as follows: analyzing a protein energy network graph, and using a residue interaction matrix as an input adjacency matrix of a graph convolution neural network, wherein the dimension of the matrix is N x N, and N is the number of protein residues; the other is a feature matrix H, the dimension of the matrix H is N x F, wherein N represents the number of residues in the graph, F is the feature dimension of each residue, and the dimension records the amino acid features in the residues.

Furthermore, in the topological structure of the protein energy network graph, the edges are weighted, the information of the edges is needed to be embedded into graph convolution, and the iteration and the circulation of the graph convolution are carried out twice, wherein the formula is as follows:

wherein X ε R ^N*C ，Θ∈R ^C*F’ ，Z∈R ^N*F’ N, C, F' represent the number of residues, the number of channels and the number of convolution kernels respectively; wherein the method comprises the steps ofAnd->Is further constrained by an adjacency matrix a and a degree matrix D;

and then using a sigmoid loss function to semi-supervise and classify the protein energy network diagram, wherein the formula is as follows:

wherein y is _lf Is the classified residueBase set, Y _lf Is a classification matrix, Z _lf The result is output by the graph convolution network, F represents the number of residues, and L is the classification result.

In the topological structure of the protein energy network graph, the edges are weighted, if the weight value is singly placed in the adjacent matrix (namely, the weight information is placed in the point), only the product of the adjacent matrix cannot highlight the information of the edges, and the information of the edges needs to be embedded into the graph convolution to perform iteration and circulation of the graph convolution twice. This is done to take the side information into account, support computation of sparse matrices, and semi-supervised classification of the graph from multiple angles (i.e. considering the sides and also considering the points) followed by using the sigmoid loss function.

Compared with the prior art, the invention has the beneficial effects that:

(1) The invention provides a group of molecular dynamics result analysis methods based on protein residue interaction and protein energy network analysis, which take structural characteristics of the compounds into consideration, reflect the integral combination degree of targets and the compounds, can analyze various characteristics of results and ensure the comprehensiveness of molecular dynamics result analysis.

(2) According to the method, the analysis and evaluation are carried out on the molecular dynamics result by fusing the method of the graph rolling neural network, so that the correlation can be greatly improved, the accuracy and the high efficiency of analysis of the molecular dynamics result are ensured, and the analysis of the molecular dynamics data by the graph rolling neural network is an attractive method with huge prediction capability.

Drawings

FIG. 1 is a flow chart of the method of the present invention;

FIG. 2 is a diagram of a graph roll-up neural network architecture;

FIG. 3 is a graph showing the results of the MM-PBSA method binding energy estimation;

fig. 4 is a graph showing the result of the method binding energy estimation of the GCN of the present invention.

Detailed Description

The present invention will be further described with reference to specific embodiments thereof, wherein it is apparent that the embodiments described are merely some, but not all, of the embodiments of the present invention. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.

Example 1

A method for analyzing molecular dynamics results based on residue interaction and PEN, as shown in fig. 1, wherein the method specifically comprises the following steps:

1) Construction of residue interaction energy matrix

Tankyrase polymerase, which is useful for the treatment of colorectal cancer and has a known value of maximum half inhibition concentration (IC 50), is interfaced with 57 small molecules, and after molecular docking, the highest scored complex conformation is taken as input, molecular dynamics simulation is performed, and the whole simulation time is set to a simulation trajectory of 10ns to ensure better energy convergence and state space sampling.

In the simulated file, the solvent molecules were removed. And extracting residue-residue interaction energy from the simulated trajectory obtained after simulation by using a tool, and simultaneously calculating correlations between the residue-residue interaction energy and the residue interaction energy to obtain paired residue interaction energy and correlations between all possible residue pairs in the structure.

The interaction energy between residues i and j is the sum of the non-bond interaction energies defined in the force field. Non-bond interaction energy is generally considered by van der Waals forces and electrostatic terms:

here, i and j each represent two residues, E _ij Is an array that calculates the interaction energy between residues i and j in all trajectory frames,for electrostatic interaction>Is Van der Waals force interactionActing as a medicine. Wherein interactions between individual residues can be used to construct IEMs (interaction energy matrices) that represent the "intensities" of edges between nodes (residues) in a network structure.

S _ij ＝|E _ij ，N _ij |

Here, S _ij Representing an interaction energy matrix, E _ij Is the interaction energy between residues i and j in all trace frames, N _ij Is the number of residues in the trace frame.

2) Construction of protein energy network structure diagram

The IEM of the "strength" of edges between nodes (residues) in a network structure may be normalized to further represent the "weight" or "cost of information transfer" used from one node to another in a network analysis task, and the resulting Interaction Energy Matrix (IEM) is used to construct and analyze a Protein Energy Network (PEN) based on residue-based network metrics (e.g., degree, mid-center, and proximity-center) and the shortest path between selected residues in the structure.

The protein structure network of complex modeling data is constructed using the energy of interactions of pairs of residues from the trace, with individual residues as nodes, and the average energy of interactions between each pair of residues as "weights" added to the edges between the nodes of the residues.

The information contained in IEM (interaction energy matrix) is used to construct a protein energy network of complex modeling data, with each residue in the structure being considered a node. If there is a non-zero average interaction energy between two corresponding residues in the IEM, an edge is added between the two nodes (residues). Using the value of the average IE, the "weight" attribute of the edge is determined using the following formula: wherein ω represents the edge weight between residues i and j, and if residues i and j are covalent bonds, the weight between them is 0.99, otherwise the weight between them is x _ij ，x _ij Representing the average interaction energy between residues i and j;

3) Analysis using graph convolution neural network classification training

As shown in fig. 2, a protein energy network graph is analyzed, using a residue interaction matrix as an input adjacency matrix for graph convolution neural network, the dimension of the matrix is n×n, where N is the number of protein residues. The other is an input feature matrix H, the dimension of the matrix H is N x F, wherein N represents the number of nodes in the graph, F is the feature dimension of each node, and the dimension records the amino acid features in the residues.

In the topological structure of the protein energy network graph, the edges are weighted, if the weight value is singly placed in the adjacent matrix (namely, the weight information is placed in the point), only the product of the adjacent matrix cannot highlight the information of the edges, the information of the edges needs to be embedded into the graph convolution, and the iteration and the circulation of the graph convolution are carried out twice, wherein the formula is as follows:

wherein X ε R ^N*C ，Θ∈R ^C*F’ ，Z∈R ^N*F’ N, C, F' represent the number of nodes, the number of channels and the number of convolution kernels respectively; wherein the method comprises the steps ofAnd->Is further constrained by the adjacency matrix a and the degree matrix D.

This is done to take the information of the edges into account, support the computation of sparse matrices, and take into account from multiple angles (i.e. consider edges and also consider points). And then semi-supervised classification is carried out on the graph by using a sigmoid loss function, wherein the formula is as follows:

wherein y is _lf Is a categorized residue set, Y _lf Is a classification matrix, Z _lf The result is output by the graph convolution network, F represents the number of residues, and L is the classification result.

In order to compare the molecular dynamics result analysis method based on protein residue interactions and protein energy network analysis with the widely used alternative method in virtual screening, two types of comparison methods were selected in this example. Another approach is the MM-PBSA method, which represents a more stringent but computationally expensive binding energy estimation method.

The comparison index is selected from the Spearman coefficients and Pearson correlation coefficients.

I.e. the ratio of the covariance of the two variables X, Y to the product of the standard deviations of the two variables, the standard deviation of the two variables X, Y cannot be zero. Spearman (sample correlation coefficient):

the data of the two variables (X, Y) are ordered (unified ascending order or descending order), the position of each variable after the ordering is the rank order (X ', Y'), the difference value of the rank orders X ', Y' of the X, Y with the same original position is di, and n is the number of the variables.

Both correlation coefficients (Pearson, spearman) reflect the direction and extent of the trend of the change between the two variables, with values ranging from-1 to +1,0 indicating that the two variables are uncorrelated, positive values indicating positive correlation, negative values indicating negative correlation, and larger absolute values indicating stronger correlation.

The graphs of the experimental results are shown in fig. 3 and 4.

The Spearman and Pearson correlation coefficients between the inhibitor pIC50 values and the affinity estimates are summarized in table 1.

TABLE 1 correlation coefficient results comparison Table

Method	MM-PBSA	GCN
			Spearman	-0.49	0.67
Pearson	-0.34	0.71

As can be seen from the table, the correlation based on the MM-PBSA method can obtain a good correlation coefficient, compared with the correlation based on the GCN method, the correlation can be greatly improved, and compared with the traditional energy analysis method, the combination energy prediction accuracy of the molecular dynamics result analysis method based on protein residue interaction and protein energy network analysis is higher.

Claims (5)

1. A method for analyzing molecular dynamics results based on residue interactions and PEN, characterized in that the method comprises the steps of:

1) Constructing a residue interaction energy matrix;

2) Constructing and analyzing a protein energy network structure diagram;

3) Classifying and training a graph convolution neural network;

the specific operation of the step 1) is as follows:

extracting residue-to-residue interaction energy from the simulated trajectory using a tool while calculating correlations between them, wherein interactions between individual residues can be used to construct an interaction energy matrix representing the "strength" of edges between residues in the network structure, the interaction energy matrix being abbreviated IEM;

the specific operation of the step 2) is as follows: the IEMs of the "intensities" of the edges between residues in the network structure are normalized to further represent the "weights" or "information transfer costs" used from one residue to another in the network analysis task, and the resulting interaction energy matrix is used to construct and analyze a protein energy network, abbreviated PEN, based on the residue-based network metrics and the shortest paths between selected residues in the structure;

the specific operation of the step 3) is as follows: analyzing a protein energy network graph, and using a residue interaction matrix as an input adjacency matrix of the graph-convolution neural network, wherein the adjacency matrix has a dimension of N x N, and N is the number of protein residues; the other is a feature matrix H, the dimension of the feature matrix H is N x F, wherein N represents the number of residues in the protein energy network diagram, F is the feature dimension of each residue, and the feature dimension records the amino acid features in the residues.

2. A method of analysis of molecular dynamics based on residue interactions and PEN according to claim 1, characterized in that the interaction energy between individual residues is the sum of the non-bond interaction energies defined in the force field; non-bond interactions can be considered with van der Waals forces and electrostatic terms:

here, i and j each represent two residues, E _ij Is an array that calculates the interaction energy between residues i and j in all trajectory frames,is in electrostatic phaseInteraction, I/O (administration of the drug)>Is van der waals interactions.

3. A method of analysis of molecular dynamics based on residue interactions and PEN according to claim 1, characterized in that the interactions between individual residues can be used to construct IEMs representing the "strength" of edges between residues in the network structure;

S _ij ＝|E _ij ，N _ij i, here S _ij Representing an interaction energy matrix, E _ij Is the interaction energy between residues i and j in all trace frames, N _ij Is the number of residues in the trace frame.

4. The method of claim 1, wherein step 2) uses the following formula to determine the "weight" attribute of the edge:

wherein ω represents the edge weight between residues i and j, and if residues i and j are covalent bonds, the weight between them is 0.99, otherwise the weight between them is x _ij ，x _ij Representing the average interaction energy between residues i and j.

5. The method for analyzing molecular dynamics results based on residue interactions and PEN according to claim 1, wherein step 3) in the topological structure of the protein energy network graph, edges are weighted, and information of the edges is needed to be embedded into graph convolution to perform iteration and loop of two graph convolution, and the formula is as follows:

wherein X ε R ^N*C ，Θ∈R ^C*F’ ，Z∈R ^N*F’ N, C, F' represent the number of residues, the number of channels and the number of convolution kernels respectively; wherein the method comprises the steps ofAnd->Is further constrained by an adjacency matrix a and a degree matrix D;

and then using a sigmoid loss function to semi-supervise and classify the protein energy network diagram, wherein the formula is as follows:

wherein y is _lf Is a categorized residue set, Y _lf Is a classification matrix, Z _lf The result is output by the graph convolution network, F represents the number of residues, and L is the classification result.