A personalized parameter generation method for protein force field simulation

Abstract

The present invention relates to a personalized parameter generation method for protein force field simulation, which uses structural data in the PDB database for preprocessing and calculates sequence features. These features are used as input to train a bidirectional long short-term memory recurrent neural network to obtain the predicted values ​​of ψ and the distribution area of ​​the residue dihedral angle on the Laplace diagram, and the CMAP parameters of each residue are optimized by a reweighted method, so as to finally achieve the purpose of generating personalized force field parameters. Compared with the prior art, the present invention has the advantages of improving the accuracy of force field simulation and the like.

Description

Technical Field

[0001] The present invention relates to the technical field of bioinformation computer processing, and in particular to a personalized parameter generation method for protein force field simulation. Background Art

[0002] The structure of a protein determines its function. Natural random proteins have diverse flexible conformations. It is very challenging to experimentally describe the combination of conformations. Molecular dynamics (MD) simulations can sample continuous conformations, which can provide a valuable supplement to experimental data. However, MD simulations still face two major limitations: 1. Current computational efficiency limits the time scale of MD simulations to within microseconds, which is still far below many important biological time scales (sampling problem); 2. The accuracy of the current potential energy function may lead to biased conformations and energies (force field problem).

[0003] In order to improve the accuracy of the force field, adding grid-based energy correction map (Connectivity map, CMAP) parameters to distinguish amino acids has become one of the feasible methods for force field optimization and has been applied in many studies on force field development. Eight amino acid CMAP parameters that promote disorder were added to the ff99SB and ff14SB force fields, and then expanded to 20 amino acids, resulting in a set of IDP-specific force fields for structural proteins and natural disordered proteins. The environment-specific force field ESFF1 was developed based on ff14SB, adding nearly 80 CMAP parameters to describe the effect of the local environment on protein residues. Since the conformation of amino acids is determined not only by the local environment, but also by remote homologous proteins and homologous amino acids. Therefore, how to provide more accurate personalized parameters for each amino acid in the simulated protein in response to the problem of insufficient simulation accuracy of existing force fields has become a problem that needs to be solved urgently. Summary of the invention

[0004] The purpose of the present invention is to provide a personalized parameter generation method for protein force field simulation in order to overcome the defects of the above-mentioned prior art.

[0005] The purpose of the present invention can be achieved by the following technical solutions:

[0006] A personalized parameter generation method for protein force field simulation, the method comprising:

[0007] S1: Collect protein structure data and construct a sequence data set for subsequent processing; in this step, protein structure data is collected through the PDB database. The specific contents of constructing a sequence data set for subsequent processing are:

[0008] 11) For protein structure data collected from the PDB database, select those with X-ray resolution less than or equal to The data with R factor index less than or equal to 0.25 were selected, and protein structure data with sequence homology less than 30% were screened using cd-hit;

[0009] 12) Process PDB files and extract sequence information for each PDB structure. The angles, the areas corresponding to each pair of dihedral angles, and the prediction results of DSSP are used to remove small molecules and sequences other than the 20 amino acids that make up human proteins to obtain the remaining sequences. According to the length of the remaining sequences, a stratified sampling method is used to extract 10% of the data as the test group, and the remaining 90% of the data as the training group to complete the construction of the sequence data set.

[0010] S2: Perform feature extraction on the constructed sequence data set to obtain sequence features.

[0011] The extracted sequence features include a 20×sequence length vector consisting of the one-hot encoded vector of each amino acid in the remaining sequence, a position-specific scoring matrix obtained by searching the sequence library, a probability matrix for predicting residue secondary structures based on the remaining sequence, and an information matrix for predicting residue secondary structure tendencies based on the remaining sequence.

[0012] S3: A bidirectional long short-term memory recurrent neural network is used to build a deep learning model, train the model based on sequence features, and output residue dihedral angles and ψ, as well as the area over which the residue dihedral angles are distributed on the Ramachandran plot.

[0013] The deep learning model includes a classification deep learning model for predicting dihedral angle areas and a regression deep learning model for predicting dihedral angles.

[0014] Furthermore, the loss function adopted by the classification deep learning model is a cross entropy loss function, and the loss function adopted by the regression deep learning model is a sum of squares loss function.

[0015] Each deep learning model includes two layers of LSTM-BRNN networks, a fully connected hidden layer of 1024 nodes and a fully connected hidden layer of 512 nodes, and the LSTM-BRNN network is an LSTM-BRNN network with 256 nodes in each direction. Furthermore, the fully connected hidden layer adopts a ReLU activation function.

[0016] S4: The CMAP parameters are specifically corrected for the dihedral angles of the residues and the regions where they are located using the reweighting method to obtain the final personalized force field parameters, that is, to generate the optimized parameters of each residue in the protein sequence except the first and last residues. Energy parameters of 576 grid points of ψ dihedral angle λ′ >:

[0017]

[0018] In the formula, <> represents the average value, A λ′ and A λ denote the properties simulated using force field parameters λ′ and λ, respectively, and E λ′ and E λ denote the potential energy of each simulated conformation calculated using the force field parameters λ′ and λ, respectively, and β denotes the inverse of the product of the thermodynamic temperature and the Boltzmann constant.

[0019] The personalized parameter generation method for protein force field simulation provided by the present invention has at least the following beneficial effects compared with the prior art:

[0020] 1) The method of the present invention generates personalized force field parameters based on a bidirectional long short-term memory recurrent neural network and a reweighting method, and obtains more accurate force field parameters by training a bidirectional long short-term memory recurrent neural network. The area where the ψ values ​​and the residue dihedral angles are distributed on the Ramachandran plot is used to generate different CMAP parameters for each residue using the reweighting method, and enhanced sampling is performed through replica exchange of solute tempering to further improve the accuracy of the force field simulation;

[0021] 2) The method of the present invention simultaneously considers the local environment and remote homology information of each residue, can well reproduce the experimental results of the test protein, can further improve the accuracy of force field simulation, and has guiding significance for the development of more accurate force fields. The generation of personalized force field parameters will help to use molecular dynamics simulation to study the structure and function of more complex proteins in the future. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 A network diagram of a deep learning model constructed using a bidirectional long short-term memory recurrent neural network in an embodiment.

[0023] Figure 2 A schematic diagram of the process of constructing a data set in the embodiment.

[0024] Figure 3 3 is the distribution of observed Rg values ​​of the natural random protein AB40 simulated by ESFF1, ff03CMAP, PPFF1 and PPFF1_af2 force fields in the example, wherein the vertical lines represent the average Rg values ​​of the force field simulations.

[0025] Figure 4 ​Figure 1 is a diagram of the secondary structure and simulated stability results of Trp-cage in the examples, wherein Part A shows the secondary structure and PDB structure of Trp-cage predicted by PPFF1 (PDB ID: 1L2Y). In the PDB structure, the residues in the helix, turn and coil are marked with subscripts H, T and C, respectively, and the terminal residues are unmarked. In the prediction results of PPFF1, the correctly predicted residues are marked with subscripts corresponding to the secondary structure, and the rest are unmarked. Parts B and C respectively represent the comparison of the folding conformations of ESFF1 and PPFF1 with the PDB structure of Trp-cage; Part D shows the experimental data and the folding scores of Trp-cage simulated using ESFF1 and PPFF1. DETAILED DESCRIPTION

[0026] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work should belong to the scope of protection of the present invention.

[0027] Example

[0028] The embodiment of the present invention provides a method for generating personalized parameters for protein force field simulation, and the implementation process of the method includes the following steps:

[0029] 1. Dataset preprocessing

[0030] The data set comes from all protein structure data in the PDB database in May 2019, with a total of 151,310 structure files. ( The R factor index is less than or equal to 0.25, and cd-hit is used to filter out protein structure data with sequence homology less than 30%. The MDTraj python package is used to process the pdb file and extract the sequence information of each pdb structure. Angles, the area corresponding to each pair of dihedral angles, and the prediction results of DSSP, small molecules were removed, and sequences other than the 20 amino acids that make up human proteins were removed. Finally, 10,651 sequences were obtained. According to the length of the sequence, 10% of the data was extracted as the test group using the stratified sampling method, and the remaining 90% of the data was used as the training group. Figure 2 The process of dataset construction is shown.

[0031] 2. Feature Calculation

[0032] The features include a 20×L (L is the sequence length) vector (SEQ) consisting of a one-hot encoding vector of each amino acid in the sequence, a position-specific scoring matrix (PSSM) obtained by searching the sequence library using PSI-BLAST, a probability matrix of residue secondary structures predicted by the PSIPRED software based on sequence prediction, and an information matrix of residue secondary structure tendencies predicted by the SPIDER3-Single software based on sequence prediction.

[0033] 3. Model Architecture

[0034] A deep learning model is constructed using a bidirectional long short-term memory recurrent neural network. First, the model takes the features of each sequence as input and finally outputs The area where the ψ values ​​and residue dihedral angles are distributed on the Ramachandran plot. The model architecture is as follows Figure 1 As shown in the figure, it includes two deep learning models, classification and regression. The first is a two-layer LSTM-BRNN network with 256 nodes in each direction, followed by two fully connected hidden layers with 1024 nodes and 512 nodes. The fully connected layer uses the ReLU activation function and sets the dropout to 0.8. Furthermore, Figure 1 In the above, the classification network is used to predict the dihedral angle area, and the regression network is used to predict the two sets of dihedral angles, and then the prediction results of the two categories are output respectively, that is, the The area of ​​the Ramachandran plot where the dihedral angles of each residue are distributed can be the area setting standard in the prior art and in the paper by Huang et al. (Nature Methods, 2017, 14: 71-73), where the α' and αR regions included in the α+ region are removed when the α+ region is considered.

[0035] IV. Model Training and Evaluation

[0036] The entire network is implemented using FAIR's PyTorch library, and the CUDA version of the NVIDIA Tesla V100 GPU is used to accelerate training. In model training, the cross entropy loss function is used for classification problems, and the square sum loss function is used for regression problems. 20,000 epochs are trained. Considering the training speed and computer memory, the batch size (mini-batchsize) is set to 50. The initial learning rate of the Adam optimizer is set to 0.001. The performance of the model is evaluated by the MAE (mean absolute error) of the ψ angle and the accuracy of the region.

[0037] 5. Reweighted Optimization of CMAP Parameters

[0038] Normally, optimized force field parameters will make the proportion of various conformations in the simulated conformation set more reasonable, and different types of simulation values ​​will be closer to the experimental values. Reweighting can link force field parameters with simulation properties that have no explicit relationship with the parameters, so that the properties simulated by the current parameters can be used to estimate the simulation properties under new parameters. The present invention uses a reweighting method to perform CMAP parameter-specific correction on the dihedral angles of the residues and the region in which they are located to obtain the final personalized force field parameters (PPFF1). Specifically, the present invention uses a reweighting method based on the output of the LSTM-BRNN model. The optimized corners and the regions where they are located are generated for each residue in the protein sequence (except the first and last residues). The energy (CMAP) parameters of the 576 grid points of the dihedral angle are obtained as follows:

[0039] For a set of conformations that converge under specific force field parameters, the value of its average property A under the new force field parameters can be expressed by the following formula (1):

[0040]

[0041] Among them, <> represents the average value, A λ′ and A λ denote the properties simulated using force field parameters λ′ and λ, respectively, and E λ′ and E λ denote the potential energy of each simulated conformation calculated using the force field parameters λ′ and λ, respectively, and β denotes the inverse of the product of the thermodynamic temperature and the Boltzmann constant.

[0042] 6. Molecular dynamics simulation

[0043] In order to solve the problem of efficient sampling of the test system, the solute tempered replica exchange (REST2) method is proposed to obtain as comprehensive conformational coverage as possible under limited computing resources. GROMACS 2018.8 and plumed 2.6.2 software are used to run REST2. The preprocessing of the simulation system is as follows:

[0044] First, water molecules were added to the initial structure of the natural random protein to simulate the solvent environment, and counterions were added to maintain the electrical neutrality of the system. Then, the system was energy optimized for 4000 steps using the steepest descent method, and then the system was heated to 298K under the NVT ensemble. Finally, energy balance was performed for 500ps under the NPT ensemble. Throughout the process, all covalent bonds connected to hydrogen atoms were constrained using the LINCS algorithm, and electrostatic interactions were constrained using the PME method. The cutoff values ​​for the Lennard-Jones potential energy and electrostatic interaction energy were set to The grid size is A fourth-order interpolation function is used. The temperature of the system is coupled by velocity scaling and random terms, and the pressure of the system is adjusted by the Parrinello-Rahman algorithm. In the simulation phase after the pretreatment, 16 parallel simulations are used to exchange conformations at a certain ratio (exchange replicas), and the "effective temperature" ranges from 298K to 500K. According to the principle of REST2, the temperature of different replicas i can be set using formula (2):

[0045]

[0046] Among them, T eff,i represents the temperature of different replicas i, T eff,0 It is the lowest temperature and the temperature of focus. eff,max represents the maximum temperature, and n represents the number of copies.

[0047] According to the above formula, the temperatures of the 16 replicas can be calculated to be 298.00 K, 308.46 K, 319.29 K, 330.50 K, 342.10 K, 354.11 K, 366.54 K, 379.40 K, 392.72 K, 406.51 K, 420.78 K, 435.55 K, 450.84 K, 466.66 K, 483.04 K, and 500.00 K. The simulation time for each replica is 100 ns, the integration step is 2 fs, a frame of conformation is output every 50 ps, ​​and a replica exchange attempt is made every 1000 steps (2 ps).

[0048] VII. Molecular dynamics simulation analysis

[0049] The root mean square deviation (RMSD) and radius of gyration (Rg) were calculated using the cpptraj module of AMBER. SPARTA+ software was used to calculate chemical shifts. The J coupling constant was calculated using the Karplus equation. In order to compare the accuracy of the force field, the present invention uses a normalized force field score. For folded proteins, equation (3) is used to calculate the average normalized RMSD of all experimental chemical shifts.

[0050]

[0051] Among them, Folded Proteins FF score is the average normalized value of all experimental chemical shifts, N is the number of skeleton chemical shift types, and FF is rmsd is the RMSD of the simulated value corresponding to the experimental value of type j, rmsd Norm is the lowest RMSD of the jth class among all force fields.

[0052] For natural disordered proteins, the force field score Disordered Protein FF score is calculated as in equation (4), which is the backbone NMR chemical shift score (CS score, CS score ), scores based on other NMR measurements (NMR scores, NMR score ) and Rg deviation penalty (Rg Penalty ) where the CS score is calculated in the same way as the folded protein score. For all other categories of NMR measurements, the NMR score is similar to the CS score. If the simulated Rg mean is within the experimental estimate error, then Rg Penalty If the deviation is greater than the experimentally estimated error, use equation (5) to calculate Rg Penalty Therefore, the closer the combined force field score is to 1, the closer it is to the experimental observations.

[0053]

[0054]

[0055] In the formula, Rg Exp Represents the experimental mean square radius of gyration, Rg sim Represents the simulated mean square radius of gyration, Rg Exp error Representative experimental mean square radius of gyration deviation.

[0056] Based on the above-designed method steps, in order to prove the effectiveness and accuracy of the method of the present invention, the example verification content of this embodiment is as follows:

[0057] (1) According to the above data set preprocessing steps, that is, combining Figure 2 The data set was constructed using the following process, and the obtained results include 10,651 sequences. The training set and test set have 9,586 sequences and 1,065 sequences respectively.

[0058] (2) Benchmark model: For comparison, this embodiment establishes a benchmark model, namely a Transformer Encoder model with a self-attention mechanism and compares it with the AlphaFold2 model. The Transformer Encoder model includes two deep learning models for classification and regression. First, there are two layers of Transformer Encoder networks, followed by two fully connected hidden layers of 1024 nodes and 512 nodes. The fully connected layers use the ReLU activation function and set the dropout to 0.8. The AlphaFold2 model is directly used to calculate the MAE of the dihedral angles and the accuracy of the dihedral angle area for one-third of the test set.

[0059] (3) Model performance evaluation

[0060] This example first trained the LSTM-BRNN model. After one round of iteration, the performance of the model improved, but the improvement was very small, so the model that had undergone two rounds of iterations was finally selected, and the subsequent Transformer Encoder model was not iterated. Then the LSTM-BRNN model was compared with the Transformer Encoder model and the AlphaFold2 model. Table 1 summarizes the performance of the three models on the test set.

[0061] Table 1: Performance of LSTM-BRNN model, Transformer Encoder model and AlphaFold2 model for predicting amino acid dihedral angles and dihedral angle distribution areas

[0062]

[0063] As shown in Table 1, among the three models, the AlphaFold2 model performs best in predicting dihedral angles and dihedral angle distribution areas. The MAE of the angle is 10.84, the MAE of the ψ angle is 12.62, and the prediction accuracy of the dihedral angle distribution area is 85.83%, followed by the LSTM-BRNN model, The MAE of the angle is 17.79, the MAE of the ψ angle is 26.36, and the prediction accuracy of the dihedral angle distribution area is 74.8%. However, because the prediction accuracy of AlphaFold2 in natural random proteins is not high, its model performance may be good because the prediction accuracy of structural proteins is high, thereby improving the overall model performance. Therefore, molecular dynamics of natural random proteins, structural proteins and rapidly folding proteins were performed after this example to test the accuracy of the PPFF1 force field.

[0064] (4) Protein force field accuracy test

[0065] (4.1) Natural random proteins

[0066] Aβ40 is a natural random protein, which is a 40-amino acid-long amyloid β-protein peptide. In order to describe the compactness of the conformation in MD simulations, the overall average Rg is an effective parameter. In this work, the TIP4PEW water model was used, and the environment-specific molecular force field (ESFF1), the molecular force field of structural proteins and natural random proteins in equilibrium (ff03CMAP), the personalized molecular force field (PPFF1), and the personalized molecular force field optimized by AlphaFold (PPFF1_af2) were used to calculate the average Rg value of Aβ40. Figure 3 As shown, the Rg values ​​of the four force fields cover a wide range, indicating that the simulation sampling is sufficient. and PPFF1 force field The Rg averages of the force field PPFF1_af2 are within the experimental error range. The Rg value of Aβ40 is overestimated. In addition, compared with ESFF1 (0.48), ff03CMAP (0.54), and PPFF1_af2 force fields (0.41), the residual polar coupling constant (RDC) of the PPFF1 force field is more consistent with the experimental data and has a lower Q factor (0.39). This result shows the superiority of PPFF1 in describing protein residues by adding personalized parameters.

[0067] Furthermore, this example analyzes the residual environment of the Aβ40 sequence according to the criteria proposed by the ESFF1 force field, and residues with the same environment are marked with the same color. Only Val and Gly are considered. Because these two residues account for a relatively high proportion in the sequence. In ESFF1, the same residual environment has the same CMAP parameters, while the CMAP parameters of these residues in PPFF1 are different. During the simulation, it was found that the CMAP parameters of the PPFF1 force field α The RMS error of the chemical shift of α The RMS error of the chemical shifts is 0.95 ppm. It is clear that the major differences between ESFF1 and PPFF1 are caused by these residues. The same results are observed in the chemical shifts of the backbone N atoms. However, the PPFF1 side chain C β The chemical shifts of the atoms in PPFF1 are slightly worse than those in ESFF1, indicating that the side chain dihedral angle parameters need to be modified in future force field optimization. At the same time, as shown in Table 2, this example calculates the CSscore, NMRscore and FFscore of the four force fields. The CSscore and NMRscore of the PPFF1 force field are the lowest among the four force fields, and the force field score is closest to 1. These results show that personalized CMAP parameters are more accurate in describing complex proteins with precise residue environments.

[0068] Table 2: RMSD of secondary chemical shifts and J coupling constants of Aβ40 simulated by four molecular force fields and force field scores

[0069]

[0070]

[0071] (4.2) Structural proteins

[0072] Ubiquitin is a folded protein with a large number of experimental measurements and has been extensively tested by many force fields in previous work. In this work, this example uses this protein to evaluate the PPFF1 force field and compares the simulation results with the ESFF1 and ff03CMAP force fields. This example calculates the force field scores of the three force fields, and it can be seen that the PPFF1 force field score is 1.37, which is higher than the ESFF1 force field score (1.00) and the ff03CMAP force field score (1.17). The simulation values ​​are then compared with the experimental values, including secondary chemical shifts and J coupling constants. It can be seen that the ESFF1 force field secondary chemical shift simulation value is closer to the experimental value, the ff03CMAP force field J coupling constant simulation value is closer to the experimental value, and the PPFF1 force field performs poorly on the ubiquitin protein, indicating that the PPFF1 force field is more suitable for simulating natural random proteins, as shown in Table 3.

[0073] Table 3: RMSD of secondary chemical shifts and J coupling constants of ubiquitin simulated by three force fields and force field scores

[0074]

[0075]

[0076] (4.3) Rapidly folding proteins

[0077] Trp-cage is one of the typical fast-folding proteins and has been simulated in many force field development and improvement works. In previous works, Trp-cage was simulated by a combination of several force fields and water models, and the results showed that only a few force fields can capture the temperature-dependent folding dynamics well. Here, this example uses REMD simulations of Trp-cage to evaluate the PPFF1 force field. Figure 4 As shown in part A of Figure 1, the secondary structure prediction of PPFF1 is in good agreement with the NMR structure (PDB: 1L2Y39). In addition, the REMD simulation with the PPFF1 force field can well capture the folded conformation of Trp-cage, with the lowest RMSD of 1.04% compared with the experimental NMR structure. Compared with the previous results of ESFF1 force field, the melting curves of the new force field are closer to the experimental data.

[0078] In summary, the present invention proposes a personalized parameter generation method for protein force field simulation based on a bidirectional long short-term memory recurrent neural network and a reweighting method. This method takes into account both the local environment and remote homology information of each residue and can well reproduce the experimental results of the test protein. This method has guiding significance for the development of more accurate force fields, and the idea of ​​personalized force fields will help to study the structure and function of more complex proteins using molecular dynamics simulations in the future.

[0079] The above is only a specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any person familiar with the technical field can easily think of various equivalent modifications or substitutions within the technical scope disclosed by the present invention, and these modifications or substitutions should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be based on the protection scope of the claims.

Claims

1. A personalized parameter generation method for protein force field simulation, characterized in that: include: Collect protein structure data and construct sequence data sets for subsequent processing; Perform feature extraction on the constructed sequence data set to obtain sequence features; A bidirectional long short-term memory recurrent neural network was used to construct a deep learning model. The model was trained based on sequence features and the residue dihedral angles φ and ψ were output, as well as the distribution area of ​​the residue dihedral angles on the Ramachandran plot. The dihedral angles of the residues and the regions in which they are located are subjected to CMAP parameter-specific correction using a reweighting method to obtain the final personalized force field parameters; protein structure data are collected through the PDB CNC library; The specific contents of constructing the sequence data set for subsequent processing are: 11) For protein structure data collected from the PDB database, screen out data with an X-ray resolution less than or equal to 2 Å and an R factor index less than or equal to 0.25, and use cd-hit to screen out protein structure data with sequence homology less than 30%; 12) Process the PDB files, extract the sequence information, φ / ψ angle, the area corresponding to each pair of dihedral angles, and the DSSP prediction results of each PDB structure, remove small molecules, and remove sequences other than the 20 amino acids that make up human proteins to obtain the remaining sequences. Based on the length of the remaining sequences, use the stratified sampling method to extract 10% of the data as the test group and the remaining 90% of the data as the training group to complete the construction of the sequence data set; The extracted sequence features include a vector of 20×sequence length consisting of a unique hot encoding vector of each amino acid in the remaining sequence, a position-specific scoring matrix obtained by searching the sequence library, a probability matrix for predicting residue secondary structures based on the remaining sequence, and an information matrix for predicting residue secondary structure tendencies based on the remaining sequence; the deep learning model includes a classification deep learning model for predicting dihedral angle regions and a regression deep learning model for predicting dihedral angles; The loss function used by the classification deep learning model is the cross entropy loss function, and the loss function used by the regression deep learning model is the square sum loss function; each deep learning model includes a two-layer LSTM-BRNN network, a fully connected hidden layer of 1024 nodes, and a fully connected hidden layer of 512 nodes, and the LSTM-BRNN network is an LSTM-BRNN network with 256 nodes in each direction; The fully connected hidden layer adopts the ReLU activation function; The final personalized force field parameters are energy parameters of 576 grid points of the optimized φ and ψ dihedral angles of each residue except the first and last residues in the generated protein sequence.

2. The personalized parameter generation method for protein force field simulation according to claim 1, characterized in that: The following formula is used to obtain the energy parameters of 576 grid points of the optimized φ and ψ dihedral angles of each residue in the protein sequence except the first and last residues: : ， In the formula, represents the average value, and Respectively represent the force field parameters and The simulated properties, and Respectively represent the force field parameters used and The calculated potential energy for each simulated conformation, Represents the reciprocal of the product of thermodynamic temperature and the Boltzmann constant.

Images