Method for dynamically distributing RNA LJ parameters through GBeta-PM applying model

Abstract

The invention discloses a method for dynamically distributing RNA (Ribonucleic Acid) LJ (Language Junction) parameters by a GBeta-PM applying model, which comprises the following steps of: acquiring all RNA three-dimensional structures from a PDB (Proportional Data Block) database, and processing to form two kinds of RNA nucleotide fragments, namely paired RNA nucleotide fragments and stacked RNA nucleotide fragments; making a high-precision RNA quantum chemistry level data set, and carrying out quantitative calculation on the two RNA nucleotide fragments to obtain intermolecular interaction energy; preparing a data set for training a GBeta model, and constructing a deep learning model GBeta based on a graph neural network GraphSAGE model; based on a particle swarm optimization algorithm, an LJ mapping function theory and molecular dynamics simulation, obtaining an optimal mapping function model PM (mapping); and obtaining an RNA structure of a to-be-allocated LJ parameter, converting the RNA structure of the to-be-allocated parameter into a to-be-predicted matrix, predicting the to-be-predicted matrix through the trained deep learning model to obtain an electron cloud density attenuation coefficient (beta) of each atom in an RNA molecule, and transmitting the beta into a PMapping model to calculate the LJ parameter of each atom.

Description

Technical Field

[0001] The invention belongs to the technical field of biochemistry, and in particular relates to a method for dynamically allocating LJ parameters of RNA in a molecular force field. Background Art

[0002] Ribonucleic acid (RNA) plays a crucial role in living organisms, including catalysis, gene expression regulation, and disease treatment. RNA research can also serve as a breakthrough for understanding protein structure and function, as well as the genetic information contained in DNA sequences. While experimental methods such as nuclear magnetic resonance (NMR), X-ray diffraction, and cryo-electron microscopy have been used to uncover the secrets of RNA, more tools are needed to further elucidate the mechanisms by which RNA performs its various functions in living organisms. Computational chemistry methods are an important tool for studying the molecular properties of RNA.

[0003] Common computational chemistry methods include quantum mechanics (QM) and molecular mechanics (MM). Quantum mechanics methods are accurate. Among them, density functional theory (DFT) is the mainstream theoretical method for quantum mechanics calculations. However, first-principles calculations are complicated and difficult to balance computational accuracy and cost, making them only applicable to small-scale molecular systems. Molecular dynamics methods greatly save computational time and extend the computational scale by introducing molecular force field potential functions. They are widely used in studies of RNA molecule folding and metal ion-RNA interaction mechanisms.

[0004] The use of molecular dynamics methods to study the relationship between RNA structure and function still faces many limitations. Among them, traditional molecular force fields, including Amber, GAFF, and CHARMM, use fixed parameters to simulate different molecular systems, and their parameters cannot be dynamically adjusted according to changes in molecular structure. RNA molecular structures are diverse and complex, and the electron cloud density distribution of each atom shows specific differences as the environment in which it is located changes. Fixed molecular force field parameters cannot be applied to all RNA molecules, making it difficult to obtain accurate results during the simulation process, and even large deviations may occur. Therefore, it is necessary to assign different force field parameters to different RNA molecular structures to obtain more accurate simulation results.

[0005] Currently, conventional molecular force fields classify atoms according to their chemical environment and assign specific Lennard-Jones (LJ) parameters to each type. This approach effectively reduces computational effort and can accurately reproduce some experimental data. However, it also has the following shortcomings: (1) It cannot accurately describe the electron cloud density distribution of each atom. (2) Existing molecular force field LJ parameters are usually obtained by empirically fitting limited experimental data, which limits their application to a wider range of experimental data. Summary of the Invention

[0006] To solve the above technical problems, the present invention proposes a method for dynamically allocating RNA LJ parameters using a GBeta-PMapping model to solve the above problems in the prior art.

[0007] To achieve the above object, the present invention provides a method for dynamically allocating RNA LJ parameters using a GBeta-PMapping model, comprising:

[0008] The three-dimensional structure of all RNAs obtained from the PDB database was processed to form paired and stacked RNA nucleotide fragments;

[0009] Produce a high-precision RNA quantum chemistry level data set, perform quantitative calculations on the two RNA nucleotide fragments, and obtain intermolecular interaction energy;

[0010] Prepare a dataset for training the GBeta model and build a deep learning model GBeta based on the graph neural network GraphSAGE model;

[0011] Based on particle swarm optimization algorithm, LJ mapping function theory and molecular dynamics simulation, the optimal mapping function model PMapping is obtained;

[0012] The RNA structure to which the LJ parameters are to be assigned is obtained, and the RNA structure to which the parameters are to be assigned is converted into a matrix to be predicted. The matrix to be predicted is predicted by a trained deep learning model to obtain the electron cloud density attenuation coefficient (beta) of each atom in the RNA molecule. Beta is passed into the PMapping model to calculate the LJ parameters of each atom.

[0013] Optionally, the processing of the RNA three-dimensional structure includes two processes: cutting and saturation. The cutting process uses the Python script program written in this experiment to cut, and the cutting is done from top to bottom in sequence without skipping some nucleotides in the middle. The paired nucleotide fragment refers to the two base-paired dimers after the double-stranded RNA is cut, and the stacked nucleotide fragment refers to the two base stacked dimers after the single-stranded RNA is cut (or the double-stranded RNA is first split into two single strands); the saturation process uses Discovery Studio software for automatic saturation. For the stacked nucleotide fragments, because the distance is too close, the phosphate group of the nucleotide below is removed and then hydrogen atoms are added for saturation. For the paired nucleotide fragments, hydrogen atoms are added directly at the bond break point where the program cuts for saturation.

[0014] Optionally, the high-precision RNA quantum chemical level data set is prepared by performing quantum chemical calculations using ORCA software, and the stacked nucleotide fragments are directly calculated using the calculation program provided by the software to obtain the intermolecular interaction energy at the QM level of the stacked fragments; in order to ensure the high accuracy of the paired fragments, BSSE correction is added, and ORCA is used to first calculate the single point energy of three systems, namely the single point energy of the entire paired nucleotide fragment, the single point energy of the left nucleotide shielded, and the single point energy of the right nucleotide shielded. The intermolecular interaction energy of the paired fragments is calculated using the BSSE formula:

[0015] E QM =E AB -E A(AB) -E B(AB)

[0016] Where E_QM is the intermolecular interaction energy of the paired fragments, E_AB is the single-point energy of the paired nucleotide fragments, E_A(AB) is the single-point energy of shielding the A nucleotide, and E_B(AB) is the single-point energy of shielding the B nucleotide;

[0017] As the target data for PMapping model optimization, the QM-level intermolecular interaction energy of stacked and paired nucleotide fragments is compared with the molecular mechanics (MM)-level intermolecular interaction energy calculated by molecular dynamics simulation in each iteration of the particle swarm algorithm as the true value.

[0018] Optionally, the two RNA nucleotide fragment structure vectorizations include atomic feature vectorization, bond feature vectorization between atoms, and construction of an edge adjacency matrix;

[0019] Among them, atomic feature vector quantization is performed by using the RDKit toolkit to convert atomic feature vectors for atoms in the nucleotide fragment. The atomic feature vectors include features such as atomic coordinates, charge, mass, element category, hybridization type, and whether it is an aromatic atom;

[0020] The bond feature vectors between atoms are quantified by using the RDKit toolkit to convert the structural connection relationship between atoms in the nucleotide fragment into a bond feature vector, which includes three features: bond type, whether conjugated, and whether in a ring;

[0021] The edge adjacency matrix construction process includes:

[0022] Construct the adjacency matrix based on the structural connection relationship between atoms:

[0023]

[0024] Among them, the structural connection relationship between atoms, if there is a connection relationship between the i-th atom and the j-th atom in the same nucleotide fragment, the matrix element A(i,j) in the corresponding correlation matrix in the adjacency matrix is set to 1, otherwise it is 0, where i and j represent the labels of different atoms;

[0025] The preparation of the dataset for training the GBeta model involves two steps: vectorization of the two RNA nucleotide fragment structures and calculation of the electron cloud density attenuation coefficient beta, which is the label of each nucleotide fragment;

[0026] The electron cloud density attenuation coefficient beta is calculated based on the Minimal Basis MBIS method of the atomic theory of molecules to obtain the beta value of each atom in the nucleotide fragment.

[0027] Optionally, the process of constructing the GBeta model includes:

[0028] The GBeta model is constructed based on the graph neural network GraphSAGE model, which includes two GCNConv layers and an output layer. The number of nodes in the GCNConv layer and the activation function used in each layer are modified.

[0029] The process of training the GBeta model includes:

[0030] The dataset prepared for training the GBeta model was shuffled using PyTorch's built-in shuffle() function, and the training set and test set were divided into two sets in a ratio of 8:2. Then, the GBeta model constructed above was trained.

[0031] Optionally, the LJ mapping function theory is:

[0032]

[0033]

[0034] Among them, σ i and ε i represents the LJ parameter of the i-th atom, β i represents the electron cloud density attenuation coefficient beta of the i-th atom, and The mapping function parameters of the mapping function theory are related to the type of elements, e i Represents the element type of the i-th atom (such as oxygen).

[0035] The optimal mapping function model PMapping refers to a mapping function theory that determines mapping function parameters for all elements of RNA, wherein the mapping function parameters are determined by a particle swarm optimization algorithm calling a molecular dynamics simulation.

[0036] Optionally, the matrix to be predicted includes the sum of all atomic eigenvectors, all bond eigenvectors, and edge adjacency matrices described in point 4 of the RNA.

[0037] Compared with the prior art, the present invention has the following advantages and technical effects:

[0038] This paper proposes a Gbeta-PMapping model for dynamically allocating RNA LJ parameters. This method integrates a high-precision RNA quantum chemistry dataset and quantifies the interaction energy between paired and stacked RNA nucleotide fragments. Combining graph neural networks (GraphSAGE), particle swarm optimization (PSO) algorithms, LJ mapping function theory, and molecular dynamics simulations, it achieves an optimal mapping function model, PMapping. Using the electron cloud density attenuation coefficient (beta) of each atom in the RNA molecule, the PMapping model is used to calculate the LJ parameters of each atom. Specifically:

[0039] (1) Based on a graph neural network, the present invention vectorizes the structures of paired and stacked RNA nucleotide fragments, obtains the characteristics of RNA nucleotide fragments, prepares a data set, and constructs a deep learning model GBeta. Through this model, the electron cloud density attenuation coefficient (beta) of each atom in the RNA molecule is obtained.

[0040] (2) The present invention builds an optimal mapping function model PMapping based on particle swarm optimization algorithm, LJ mapping function theory and molecular dynamics simulation, takes the beta generated by the deep learning model as input, and obtains the LJ parameters of each atom. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] The accompanying drawings, which constitute part of this application, are intended to provide a further understanding of this application. The exemplary embodiments and descriptions of this application are intended to explain this application and do not constitute an improper limitation on this application. In the accompanying drawings:

[0042] Figure 1 This is an overall flow chart of a method for dynamically allocating RNA LJ parameters using the GBeta-PMapping model based on AIM theory and graph deep learning simulation according to an embodiment of the present invention;

[0043] Figure 2 A linear regression diagram of the GBeta model prediction value and the beta value calculated by the MBIS method according to an embodiment of the present invention;

[0044] Figure 3 This is a linear regression plot of the intermolecular interaction energy calculated using the LJ parameters of the PMapping model and the LJ parameters of the GAFF force field and the intermolecular interaction energy at the QM level according to an embodiment of the present invention. DETAILED DESCRIPTION

[0045] It should be noted that, in the absence of conflict, the embodiments and features of the embodiments in this application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0046] It should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and that, although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0047] In response to the current state of the art, existing RNA force fields mostly assign LJ parameters based on atom type. This method can only assign different LJ parameters to different atom types, thus failing to accurately describe the chemical environments of different atoms, which affects the quality of RNA molecular dynamics simulations. The present invention provides the following technical solution. The goal of this invention is to assign unique LJ parameters to each atom in an RNA structure, breaking through the limitations of atom type and fully considering the distribution of electron cloud density around atoms. Based on AIM theory, mapping functions, and deep learning methods, this method assigns LJ parameters, significantly improving the quality of RNA molecular dynamics simulations.

[0048] The present invention proposes a method for dynamically allocating RNA LJ parameters using the GBeta-PMapping model, and the technical solution adopted is as follows:

[0049] Acquire RNA structural data; perform segmentation and saturation processing on the RNA structural data to form two types of stacked and paired RNA nucleotide fragments; perform feature vectorization on the processed RNA structural data to construct a feature matrix, and create a graph dataset, i.e., the dataset for training the GBeta model; the processed RNA structural data is also used to obtain a high-precision quantum chemical level dataset. Quantum chemical calculations are performed using ORCA quantization software to obtain the QM-level intermolecular interaction energy, which serves as the target data for optimizing the PMapping model parameters;

[0050] Based on the GraphSAGE deep learning model and regression task learning strategy, a GBeta model for predicting the atomic electron cloud density attenuation coefficient beta was constructed;

[0051] Based on particle swarm optimization algorithm, mapping function theory, and molecular dynamics simulation, a model for optimizing mapping function parameters is constructed;

[0052] Input the feature matrix into the GBeta model, train the GBeta model, predict the electron cloud density attenuation coefficient of each atom in the RNA structure, and evaluate the model performance;

[0053] Using high-precision quantum chemical level data sets, we optimize the mapping function parameters and obtain the best mapping function parameters, that is, the best mapping function model PMapping.

[0054] The RNA structure to which the LJ parameters are to be assigned is obtained, and the RNA structure to which the parameters are to be assigned is converted into a matrix to be predicted. The matrix to be predicted is predicted by a trained deep learning model to obtain the electron cloud density attenuation coefficient (beta) of each atom in the RNA molecule. Beta is passed into the PMapping model to calculate the LJ parameters of each atom.

[0055] In some embodiments, obtaining RNA structural data and processing the structural data to form stacked paired two nucleotide fragments include the following steps:

[0056] All RNA structure data (.mol2 files) were obtained from the PDB database;

[0057] The Python script program written in this experiment was used for cutting, cutting from top to bottom in sequence without skipping some nucleotides in the middle. Paired nucleotide fragments refer to two base-paired dimers after double-stranded RNA cutting, and stacked nucleotide fragments refer to two base stacked dimers after single-stranded RNA cutting (or double-stranded RNA can be split into two single strands first);

[0058] Automatic saturation was performed using Discovery Studio software. For stacked nucleotide fragments, because the distance was too close, the phosphate group of the lower nucleotide was removed and hydrogen atoms were added for saturation. For paired nucleotide fragments, hydrogen atoms were added directly at the bond-breaking point of the program for saturation.

[0059] In some embodiments, the processed RNA structure data is subjected to feature vector quantization to produce a graph dataset and quantum chemical calculations are performed to obtain a high-precision quantum chemical level dataset, including the following steps:

[0060] Use the RDKit toolkit to convert atomic feature vectors for the atoms in the nucleotide fragment, where the atomic feature vectors include the atomic coordinates, charge, mass, degree, number of valence electrons, element category, hybridization type, and whether it is an aromatic atom, totaling 17 dimensions;

[0061] The RDKit toolkit is used to convert the structural connection relationship between atoms in the nucleotide fragment into a bond feature vector. The bond feature vector includes three features, namely, bond type, whether it is conjugated, and whether it is in a ring, with a total of 6 dimensions. The feature of a nucleotide is the sum of all its atomic features and bond features;

[0062] The adjacency matrix of each nucleotide fragment is constructed according to the connection relationship of the atoms of the nucleotide fragment:

[0063]

[0064] Among them, the structural connection relationship between atoms, if there is a connection relationship between the i-th atom and the j-th atom in the same nucleotide fragment, the matrix element A(i,j) in the corresponding correlation matrix in the adjacency matrix is set to 1, otherwise it is 0, where i and j represent the labels of different atoms;

[0065] The ORCA quantum chemical calculation software was used to perform quantum chemical calculations for the two types of nucleotide fragments. The structural data files (.mol2 files) of the two nucleotide fragments were input into the ORCA software. The single point energy and wave function files (.gbw files) of the system were obtained at the BLYP D3 / def2-TZVPP basis set theoretical level. The intermolecular interaction energy at the QM level was then obtained by single point energy calculation.

[0066] In order to ensure high accuracy of the data, BSSE correction was added to the paired fragments, and the intermolecular interaction of the paired fragments was calculated using the BSSE formula:

[0067] E QM =E AB -E A(AB) -E B(AB)

[0068] E_QM is the intermolecular interaction energy of the paired fragments, E_AB is the single-point energy of the paired nucleotide fragments, E_A(AB) is the single-point energy of shielding nucleotide A, and E_B(AB) is the single-point energy of shielding nucleotide B. The structural data of the nucleotide fragments and their corresponding QM-level intermolecular interaction energies constitute a high-precision quantum chemical level data set;

[0069] Use the script written in this experiment based on the Multiwfn software to convert the ORCA wave function file (.gbw file) into a Gaussian wave function file (.Fchk file);

[0070] The .Fchk file of each nucleotide fragment is used as input, and the MBIS method based on AIM theory is used to calculate the beta value of each atom in the nucleotide fragment, which is the label in the graph dataset. The label of a nucleotide is the sum of the beta values of all its atoms;

[0071] The sum of the atomic features, labels, and adjacency matrices of each nucleotide constitutes a nucleotide fragment graph, and the graphs of all nucleotide fragments constitute the graph dataset for training the GBeta model.

[0072] As some embodiments, based on the GraphSAGE deep learning model and the regression task learning strategy, constructing a GBeta model for predicting the atomic electron cloud density attenuation coefficient beta includes the following steps:

[0073] Based on GraphSAGE, it includes two GCNConv layers and an output layer. The number of nodes in the GCNConv layer and the activation function used in each layer are modified. GraphSAGE can sample the nodes in each nucleotide graph and then use the aggregation function to aggregate the information of the neighboring nodes connected to the node to obtain the embedding vector of each node.

[0074] Based on the regression task learning strategy, the mean squared error function is selected as the loss function, the loss of labels and predicted values is calculated as the optimization target, and the Adam optimizer is selected to obtain the target model with the minimum loss.

[0075] As some embodiments, based on a particle swarm optimization algorithm, mapping function theory, and molecular dynamics simulation, constructing a model for optimizing mapping function parameters includes the following steps:

[0076] Based on the particle swarm optimization algorithm, a certain number of particle positions and the particle's speed are initialized. Each particle position represents a set of mapping function parameters, and the particle's speed represents the particle's own exploration inertia. Each particle must calculate the fitness of the current particle position in each iteration. The algorithm records the historical optimal solution of each particle (the position with the minimum fitness value) and the current optimal solution of all particles. After calculating the fitness in each iteration, the particle updates its position based on its own historical optimal solution, the current global optimal solution, and the particle's speed.

[0077] Based on mapping function theory and molecular dynamics simulation, the LJ parameters corresponding to the particle position are calculated and written back to the molecular file corresponding to the nucleotide fragment so that GROMACS can perform MD simulation. The fitness calculation method is modified, and the sum of the intermolecular interaction energy at the MM level calculated by molecular dynamics simulation of all nucleotide fragments and the error between the intermolecular interaction energy at the QM level is used as the fitness; for each particle position, each round of iteration requires calling GROMACS to perform MD simulation to calculate the intermolecular interaction energy at the MM level.

[0078] In some embodiments, the feature matrix, adjacency matrix, and label matrix of the nucleotide fragments are input into the deep learning model, beta values are predicted, and the model performance is evaluated, including the following process:

[0079] Input the above matrix into the GBeta deep learning model to obtain the feature vector embedded in each node;

[0080] The feature vector of each node is passed to the output layer. The output layer is a linear layer that finally outputs a value, which is the predicted beta value, using the relu activation function.

[0081] The model was trained and tested by dividing the data set in an 8:2 ratio, where 80% of the data was randomly selected as the training set and the remaining 20% as the test set. The predicted values and labels corresponding to the test set were compared to verify the performance of the model. The predicted values and labels were compared to calculate the coefficient of determination, root mean square error (RMSE), and mean absolute error (MAE) indicators, and the prediction performance of the model was evaluated based on the indicators.

[0082] As some embodiments, using a high-precision quantum chemical level data set, optimizing mapping function parameters to obtain the best mapping function parameters, that is, obtaining the best mapping function model PMapping, includes the following process:

[0083] Initialize a 20*12 matrix to store the positions of 20 particles. The 12 dimensions refer to 12 mapping function parameters (two for each element, and hydrogen is divided into polar hydrogen and non-polar hydrogen).

[0084] Initialize 20*12 and 1*12 matrices to store the historical optimal and global optimal values of each particle, initialize the particle speed and the learning coefficient of the particle swarm optimization algorithm;

[0085] The iteration termination condition is set to automatically stop the algorithm iteration when the fitness change corresponding to the last five optimal solutions is less than a certain threshold. At this time, the global optimum recorded is the optimal mapping function parameter;

[0086] The above technical solution of the present invention is described in detail with reference to the relevant drawings:

[0087] like Figure 1 As shown, the overall process of the GBeta-PMapping model dynamic allocation RNA LJ parameter method provided by the present invention includes three parts:

[0088] a) Dataset preparation and model construction: obtain RNA structural data, process RNA structural data to obtain nucleotide fragments, perform quantum chemical calculations on the nucleotide fragments to produce high-precision QM level data, perform feature vector quantization on the nucleotide fragments to produce a graph dataset for training deep learning models, construct a GBeta model for predicting beta values based on the GraphSAGE model, and construct a model for optimizing the mapping function parameters of the PMapping model based on the particle swarm optimization algorithm; b) Model training: i) construct the data of the graph dataset into a feature matrix and input it into the graph neural network to obtain the embedded feature vector of the node, and then calculate the loss and backpropagation to obtain a trained model; ii) input the high-precision QM level data as the target data into the model for optimizing the mapping function parameters, calculate the algorithm fitness, continuously update the particle position, and obtain the optimal mapping function parameters; c) Model prediction: i) predict the electron cloud density attenuation coefficient beta of the RNA structure to be assigned LJ parameters; ii) input beta into the optimal mapping function model PMapping to calculate the LJ parameters of each atom.

[0089] In the first part, regarding the assignment of LJ parameters, unlike previous studies that were limited by atom type, this paper constructs a LJ parameter assignment method based on AIM theory that fully considers the electron cloud density distribution around atoms and breaks through the atom type restriction. The construction of the method and the acquisition of the data required for the method include the following steps:

[0090] All RNA structural data were obtained from the PDB protein database. The RNA structural data were cut and saturated to form two types of stacked and paired RNA nucleotide fragments. The processed RNA structural data were feature vectorized and constructed into a feature matrix to create a graph dataset, i.e., the dataset for training the GBeta model. The processed RNA structural data were also used to obtain a high-precision quantum chemical dataset. Quantum chemical calculations were performed using the ORCA quantization software to obtain the molecular interaction energy at the QM level, which served as the target data for optimizing the PMapping model parameters.

[0091] Based on the GraphSAGE deep learning model and regression task learning strategy, a GBeta model for predicting the atomic electron cloud density attenuation coefficient beta was constructed;

[0092] Based on particle swarm optimization algorithm, mapping function theory and molecular dynamics simulation, a model for optimizing mapping function parameters is constructed.

[0093] In the second part, for model training, the graph dataset and the QM level dataset are used to train the models to obtain the GBeta model that accurately predicts the beta value and the PMapping model that obtains the best mapping function. The following steps are included:

[0094] The feature matrix, adjacency matrix, and label matrix of the nucleotide fragments in the graph dataset are input into the GraphSAGE model. After two layers of GCNConv convolutional layers, the embedding vector of each atom in the nucleotide fragment is passed to the output layer to obtain the predicted atomic beta value. The predicted value is compared with the label, and then the loss and backpropagation are calculated. The loss function is calculated as follows:

[0095]

[0096] in, The beta value of the model predicted for the i-th atom, y i is the label value of the i-th atom, n is the total number of atoms in all nucleotide fragments, and the model performance is evaluated by loss. The model with the smallest loss is selected as the GBeta model;

[0097] The beta value and QM-level dataset calculated by the MBIS method are input into the model of the optimized mapping function parameters constructed based on the particle swarm algorithm, and 20 particles and their related properties and hyperparameters are initialized; each round of iteration of each particle calls the MD simulation to obtain the MM-level intermolecular interaction energy of the LJ parameter represented by the particle (calculated by the particle mapping function parameters and beta value), and compared with the QM-level dataset, the QM-level data is used as the target data; the iteration is automatically terminated until the difference between the MM-level intermolecular interaction energy and the QM-level intermolecular interaction energy changes less than the threshold, and the currently recorded global optimum is output.

[0098] The third part, for model prediction, includes the following steps:

[0099] The RNA structure data to be predicted is obtained from the test set, and it is vectorized and converted into the corresponding feature matrix and input into the target model (GBeta model) to obtain the beta value of each atom;

[0100] To demonstrate the prediction performance of the GBeta model of the present invention, the present invention randomly selected 5 structures from the test set according to the type of paired nucleotide fragments (AG, AU, CG, GU, UU), one for each type, and drew a linear regression graph of the beta value calculated by the MBIS method for each structure and the predicted value of the GBeta model, as shown in the figure. Figure 2 As shown:

[0101] Among them, the correlation coefficient R of the CG structure 2The highest is 0.99744, and the correlation coefficients of the other structures are also higher than 0.98 and close to 1. This shows that there is a strong correlation between the beta values predicted by the GBeta model and the beta values calculated by the MBIS method, proving that the GBeta model has excellent performance and its prediction results are comparable to those of the MBIS method based on the QM wave function file.

[0102] The beta value of the structure to be predicted is input into the optimal mapping function model PMapping to obtain the unique LJ parameters of each atom.

[0103] To demonstrate the quality of the LJ parameters assigned by the present invention, the present invention uses the simulation software GROMACS to call molecular dynamics simulation to calculate the interaction energy between molecules of the structure to be predicted. The structure to be predicted will undergo two molecular dynamics simulations, one using the LJ parameters assigned by the GAFF force field (the GAFF force field is derived from the widely used AMBER force field and is an upgraded version of the AMBER force field), and the other using the LJ parameters calculated by the PMapping model; we performed MD simulations on 343 nucleotide fragments, of which 158 were paired nucleotide fragments and the rest were stacked nucleotide fragments, and drew a linear regression graph of the intermolecular interaction energy calculated by the LJ parameters of GAFF and the intermolecular interaction energy at the high-precision QM level, and a linear regression graph of the intermolecular interaction energy calculated by the LJ parameters of PMapping simulation and the intermolecular interaction energy at the QM level, as shown in FIG. Figure 3 As shown:

[0104] The LJ parameters of the PMapping model are compared with those of GAFF, and the slope and R of the paired nucleotide segments are 2 The values changed from 0.84 and 0.87 to 0.99 and 0.94, and the slope and R 2 The values changed from 1.64 and 0.72 to 0.92 and 0.84, which shows that the LJ parameters of the PMapping model can better fit the interaction energy between molecules at the QM level. The intermolecular interaction energy calculated by the LJ parameters of the PMapping model is more strongly correlated with the interaction energy between molecules at the QM level, proving that the RNA LJ parameters dynamically assigned by the GBeta-PMapping model are a better LJ parameter and can improve the effect of MD simulation.

[0105] The above are merely preferred embodiments of the present application, but the scope of protection of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in this application should be included in the scope of protection of the present application. Therefore, the scope of protection of the present application should be based on the scope of protection of the claims.

Claims

1. A method for dynamically allocating RNA LJ parameters using a GBeta-PMapping model, characterized in that: include: The three-dimensional structure of all RNAs obtained from the PDB database was processed to form paired and stacked RNA nucleotide fragments; A high-precision RNA quantum chemical level dataset is produced, and the two RNA nucleotide fragments are quantitatively calculated to obtain the intermolecular interaction energy at the quantum mechanical (QM) level of the RNA nucleotide fragments, which is used as the target data for the PMapping model optimization parameters; Vectorizing the two RNA nucleotide fragment structures to obtain RNA nucleotide fragment features and prepare a data set for training the GBeta model; A deep learning model GBeta is constructed based on the graph neural network GraphSAGE model. The trained deep learning model GBeta is obtained by training the data set of the GBeta model. Based on particle swarm optimization algorithm, LJ mapping function theory and molecular dynamics simulation, the optimal mapping function model PMapping is obtained; The RNA structure to which the LJ parameters are to be assigned is obtained, and the RNA structure to which the parameters are to be assigned is converted into a matrix to be predicted. The matrix to be predicted is predicted by a trained deep learning model to obtain the electron cloud density attenuation coefficient (beta) of each atom in the RNA molecule. Beta is passed into the PMapping model to calculate the LJ parameters of each atom.

2. The method according to claim 1, characterized in that The processing of the RNA three-dimensional structure includes two processes: cutting and saturation. The cutting process uses the Python script program written in this experiment to cut, cutting from top to bottom in sequence without skipping some nucleotides in the middle. The paired nucleotide fragment refers to the two base-paired dimers after the double-stranded RNA is cut, and the stacked nucleotide fragment refers to the two base stacked dimers after the single-stranded RNA is cut (or the double-stranded RNA is first split into two single strands); the saturation process uses Discovery Studio software for automatic saturation. For stacked nucleotide fragments, because the distance is too close, the phosphate group of the nucleotide below is removed and then hydrogen atoms are added for saturation. For paired nucleotide fragments, hydrogen atoms are added directly at the bond break point of the program for saturation.

3. The method according to claim 1, characterized in that The high-precision RNA quantum chemical level data set was prepared by performing quantum chemical calculations using the ORCA software. The stacked nucleotide fragments were directly calculated using the calculation program provided by the software to obtain the molecular interaction energy at the QM level of the stacked fragments. To ensure high data accuracy, the paired fragments were corrected using BSSE. ORCA was used to first calculate the single point energies of three systems: the single point energy of the entire paired nucleotide fragment, the single point energy of the left nucleotide shielded, and the single point energy of the right nucleotide shielded. The molecular interaction energy of the paired fragments was calculated using the BSSE formula: AND QM =And AB -AND A(AB) -AND B(AB) Among them, E QM is the intermolecular interaction energy of the paired fragments, E AB is the single-site energy of paired nucleotide segments, E A(AB) It is the single point energy that shields the A nucleotide, E B(AB) It is the single point energy that shields the B nucleotide; The target data for PMapping model optimization refers to the QM-level intermolecular interaction energy of stacked and paired nucleotide fragments, which will be compared with the molecular mechanics (MM)-level intermolecular interaction energy calculated by the molecular dynamics simulation in each iteration of the particle swarm algorithm as the actual value.

4. The method according to claim 1, wherein The two RNA nucleotide fragment structure vectorizations include atomic feature vectorizations, bond feature vectorizations between atoms, and construction of edge adjacency matrices; The atomic feature vector quantization is performed by using the RDKit toolkit to convert atomic feature vectors for atoms in the nucleotide fragment, wherein the atomic feature vectors include features such as atomic coordinates, charge, mass, element category, hybridization type, and whether it is an aromatic atom; The bond feature vectors between the atoms are quantified by using the RDKit toolkit to convert the structural connection relationship between the atoms in the nucleotide fragment into a bond feature vector. The bond feature vector includes three features: bond type, whether it is conjugated, and whether it is in a ring. The feature of a nucleotide is the sum of all its atomic features and bond features; The edge adjacency matrix construction process includes: Construct the adjacency matrix based on the structural connection relationship between atoms: Among them, the structural connection relationship between atoms, if there is a connection relationship between the i-th atom and the j-th atom in the same nucleotide fragment, the matrix element A(i,j) in the corresponding correlation matrix in the adjacency matrix is set to 1, otherwise it is 0, where i and j represent the labels of different atoms; The preparation of the dataset for training the GBeta model includes two processes: vectorization of the structures of two RNA nucleotide fragments and calculation of the electron cloud density attenuation coefficient beta, which is the label of each nucleotide fragment; The calculation of the electron cloud density attenuation coefficient beta is based on the Minimal Basis Iterative Stockholder (MBIS) method of the Atoms in Molecules (AIM) theory to obtain the beta value of each atom in the nucleotide fragment.

5. The method according to claim 1, wherein The process of constructing the GBeta model includes: The GBeta model is constructed based on the graph neural network GraphSAGE model, which includes two GCNConv layers and an output layer. The number of nodes in the GCNConv layer and the activation function used in each layer are modified. The process of training the GBeta model includes: The dataset prepared for training the GBeta model was shuffled using PyTorch's built-in shuffle() function, and the training set and test set were divided into two sets in a ratio of 8:

2. Then, the GBeta model constructed above was trained.

6. The method according to claim 1, characterized in that The LJ mapping function theory is: Among them, σ i and ε i represents the LJ parameter of the i-th atom, β i represents the electron cloud density attenuation coefficient beta of the i-th atom, and The mapping function parameter represents the mapping function theory, which is related to the type of element, e i Represents the element type of the i-th atom (such as oxygen). The optimal mapping function model PMapping refers to a mapping function theory that determines the mapping function parameters of all RNA elements, and the mapping function parameters are determined by a particle swarm optimization algorithm calling a molecular dynamics simulation.

7. The method according to claim 1, characterized in that The matrix to be predicted includes the sum of all atomic eigenvectors, all bond eigenvectors, and edge adjacency matrices described in point 4 of the RNA.

Images