protein structure elastic network model considering joint atom modeling strategy and corresponding physical force field

Abstract

This paper considers the Protein Joint Atom Modeling (SPICA) strategy and the corresponding elastic network model of physical interactions, belonging to the field of protein dynamics characterization (computational) technology. First, test sets were constructed, including 135 and 6 protein experimental crystal structures respectively. Then, the following five steps were performed: 1. Construction of the protein joint atom model; 2. Initial parameter setting and force field parameter pairing; 3. Calculation of spring force constants; 4. Construction of the Hessian matrix; 5. Parameter optimization and dynamic analysis. This method is the first to introduce the protein joint atom modeling strategy into the traditional elastic network model computational framework, while integrating the corresponding physical interaction forces, significantly improving computational accuracy while maintaining computational efficiency.

Description

Technical Field

[0001] This patent belongs to the field of protein dynamics characterization (computation) technology, and is a protein structural elastic network model developed based on a protein joint atomic modeling strategy and corresponding physical force field. Background Technology

[0002] The unique spatial conformation of proteins determines their distinctive biological functions, and they play a crucial role in maintaining the homeostasis of complex biological systems. [1, 2] The kinetic properties of proteins are widely considered a key factor in understanding their conformational-functional relationships. [3, 4] Common experimental techniques for characterizing protein kinetics include cryo-electron microscopy (cryo-EM). [5] Nuclear Magnetic Resonance (NMR) spectroscopy [6] And X-ray crystallography [7] The most representative theoretical method is molecular dynamics simulation (MD). [8] While the methods described above can provide atomic-resolution protein dynamics data, they are costly and struggle to fully cover the dynamic characteristics of proteins across their entire conformational space. Therefore, developing high-performance and cost-effective methods for protein dynamics research is of significant scientific and practical value.

[0003] To improve the computational efficiency of protein dynamics, Normal Mode Analysis (NMA) was first proposed in the 1980s. This method analyzes the protein crystal structure based on the resonant potential function, and the obtained normal vibrational modes can be used to characterize the dynamic behavior of proteins. [9] Compared to MD simulations, NMA significantly reduces computational costs, which is beneficial for high-throughput acquisition of protein dynamics data. However, it still requires the use of a full-atom model, and the input structure needs to undergo energy minimization preprocessing.

[10] With the continuous development of coarse-graining theory calculation methods [11, 12] The Elastic Network Model (ENM), a further simplification of the NMA framework, further improves computational efficiency.

[13] ENM employs a modeling strategy that uses Cα atoms to represent amino acid residues, simplifying the all-atom high-resolution protein model into a residue-scale model, thereby significantly reducing the number of atoms. Furthermore, a simple single-parameter resonant potential function has been shown to effectively capture the low-frequency dynamics of proteins, and this simplified potential function has enabled ENM to achieve a breakthrough in efficiency. Currently, the main traditional ENM variants include the Gaussian Network Model (GNM).

[14] And the Anisotropic Network Model (ANM).

[15] .

[0004] Traditional Energetic Network Model (ENM), a computational method relying solely on Cα atoms to construct models, suffers from limitations in accuracy when characterizing protein dynamics. Recent advancements in computing hardware have enabled the practical application of more computationally expensive algorithms. These factors have collectively driven researchers to integrate more refined protein feature information into the traditional ENM framework, propelling the development of novel algorithms. For example, the rigid domain elastic network model (domain-ENM) proposed by Song's team significantly improves the prediction accuracy of interdomain motion patterns within proteins by introducing the rigidity of structural domains.

[16] Similarly, Khade's team further combined the flexible characteristics of the hinge region to establish the Hinge-domain Anisotropic Network Mode (hdANM).

[17] This model can simultaneously describe the translational / rotational motion of rigid domains and the high-degree-of-freedom flexible motion of hinge regions, more realistically reproducing the conformational dynamics during allosteric processes and providing support for the study of protein function-allosteric relationships. Zhou's team introduced surface hydration energy as an additional contribution term into classical canonical mode analysis, thus incorporating solvent environment effects into protein dynamics characterization. The high-frequency vibrational modes obtained by this method are similar to the results of solvation systems obtained from molecular dynamics simulations, broadening the perspective for protein dynamics research.

[18] Kaynak et al. were the first to introduce residue specificity into the ENM computational framework and used cumulative distance-dependent force constants to describe residue-residue and residue-ligand interactions, achieving efficient characterization of protein-ligand complexes.

[19] Zhang et al. developed an improved ENM method for efficient dynamic characterization of membrane proteins by combining vibrational subsystem analysis (VSA) with recursive formulas. This method significantly reduces computational memory and time overhead by decomposing the lipid bilayer into concentric cylindrical domains and recursively solving the reduced Hessian matrix, thus efficiently characterizing the dynamic properties of large-scale membrane protein systems.

[20] Jia's team proposed a protein dynamics prediction strategy based on sequence evolution coupling: reconstructing residue contact maps through direct coupling analysis (DCA) and then deriving dynamic behavior using GNM.

[21] Guarnera and Berezovsky constructed an elastic network model integrating Miyazawa-Jernigan residue-residue potentials, established a sequence-dependent allosteric communication mechanism, and revealed two key determinants of allosteric regulation: the free energy contribution of regulatory sites and the overall structural dynamics background.

[22] .

[0005] Existing examples of improved ENM models fully demonstrate the advantages of their simplified and easily scalable computational framework. The core of these improved ENM models lies in establishing a more complete and flexible topology, enabling a more accurate description of the interactions between proteins and their surrounding environment. Therefore, introducing complex interatomic physical interactions within proteins into the traditional ENM computational framework can also be an effective way to comprehensively improve the model's computational performance. Traditional NMA successfully obtains high-precision kinetic properties by considering complex interatomic physical interactions within proteins; however, this method requires a complete protein structure. The protein structure data used in traditional NMA methods mainly comes from X-ray crystallography experiments. However, X-ray crystallography cannot determine the spatial coordinates of hydrogen atoms. Therefore, before using all-atom computational models such as NMA to analyze the experimentally measured crystal structure kinetic properties, hydrogen atoms must be artificially added to form a complete topological structure. Artificially added hydrogen atoms often place proteins in a non-equilibrium conformation, which may lead to energy saddle point problems, reducing the reliability of kinetic analysis and making the method highly sensitive to the quality of the input protein structure.

[23] Furthermore, to ensure that the protein structure input to NMA is in equilibrium, a rigorous energy-minimization preprocessing is required. This invention aims to balance the accuracy advantage of NMA due to its consideration of interatomic physical interactions with the efficiency advantage of ENM due to its simplified computational framework. To this end, it innovatively incorporates a protein joint atomic model (Surface Property fItting Coarse Grained, SPICA).

[24] This invention incorporates the traditional ENM computational process, mapping the interactions between joint atoms to an elastic topological network. The protein joint atom model is constructed based on experimentally determined heavy atom coordinates, with a resolution between that of the all-atom model and the Cα coarse-grained model. Therefore, employing a joint atom modeling strategy within the ENM framework completely avoids the problems caused by artificially adding hydrogen atoms, maintains high computational efficiency, and outperforms the traditional Cα atom modeling strategy in preserving structural details. Another key improvement of this invention lies in the specialized optimization of the joint atom interaction potential function tailored to the characteristics of ENM, thereby further improving the accuracy of the model while maintaining its efficiency advantage.

[0006] In summary, considering the combined optimization requirements of computational model accuracy and efficiency, this invention innovatively proposes a highly efficient dynamic characterization method for proteins using a joint atom modeling strategy. By introducing the joint atom modeling strategy and a matching optimized interaction potential function into the traditional ENM framework, a united atom elastic network model (uaENM) is constructed to rapidly and accurately characterize the dynamic properties of proteins. This model achieves an optimal balance between reducing model complexity (i.e., improving computational efficiency) and maintaining accuracy; its simplified elastic network topology retains key features of protein structure, comprehensively and accurately mapping protein dynamic properties, and showing significant improvements in computational performance compared to traditional ENM and all-atom models. This method can provide high-throughput dynamic information for protein functional mechanism research at a lower computational cost, while also providing a feasible technical path for the further development of elastic network models. References

[0007] [1]Fraser JS, Clarkson MW, Degnan SC, et al. Hidden alternative structures of proline isomerase essential for catalysis[J]. Nature, 2009, 462(7273): 669-673.

[0008] [2]Echenique P. Introduction to protein folding for physicists[J].Contemporary Physics, 2007, 48(2): 81-108.

[0009] [3]Liang Z, Verkhivker G M, Hu G. Integration of network models andevolutionary analysis into high-throughput modeling of protein dynamics andallosteric regulation: theory, tools and applications[J]. Briefings inbioinformatics, 2020, 21(3): 815-835.

[0010] [4]Motlagh H N, Wrabl J O, Li J, et al. The ensemble nature ofallostery[J]. Nature, 2014, 508(7496): 331-339.

[0011] [5]Fernandez-Leiro R, Scheres S H W. Unravelling biologicalmacromolecules with cryo-electron microscopy[J]. Nature, 2016, 537(7620):339-346.

[0012] [6]Fenwick R B, van den Bedem H, Fraser J S, et al. Integrateddescription of protein dynamics from room-temperature X-ray crystallographyand NMR[J]. Proceedings of the National Academy of Sciences, 2014, 111(4):E445-E454.

[0013] [7]Kohn J E, Afonine P V, Ruscio J Z, et al. Evidence of functionalprotein dynamics from X-ray crystallographic ensembles[J]. PLoS computationalbiology, 2010, 6(8): e1000911.

[0014] [8]McCammon J A, Gelin B R, Karplus M. Dynamics of folded proteins[J]. nature, 1977, 267(5612): 585-590.

[0015] [9]Go N, Noguti T, Nishikawa T. Dynamics of a small globular proteinin terms of low-frequency vibrational modes[J]. Proceedings of the NationalAcademy of Sciences, 1983, 80(12): 3696-3700.

[0016]

[10] Ma J. Usefulness and limitations of normal mode analysis inmodeling dynamics of biomolecular complexes[J]. Structure, 2005, 13(3): 373-380.

[0017]

[11] Tozzini V. Coarse-grained models for proteins[J]. Current opinionin structural biology, 2005, 15(2): 144-150.

[0018]

[12] Thorpe I F, Zhou J, Voth G A. Peptide folding using multiscalecoarse-grained models[J]. The Journal of Physical Chemistry B, 2008, 112(41):13079-13090.

[0019]

[13] Tirion M M. Large amplitude elastic motions in proteins from asingle-parameter, atomic analysis[J]. Physical review letters, 1996, 77(9):1905.

[0020]

[14] Haliloglu T, Bahar I, Erman B. Gaussian dynamics of foldedproteins[J]. Physical review letters, 1997, 79(16): 3090.

[0021]

[15] Atilgan A R, Durell S R, Jernigan R L, et al. Anisotropy offluctuation dynamics of proteins with an elastic network model[J].Biophysical journal, 2001, 80(1): 505-515.

[0022]

[16] Song G, Jernigan R L. An enhanced elastic network model torepresent the motions of domain-swapped proteins[J]. Proteins: Structure,Function, and Bioinformatics, 2006, 63(1): 197-209.

[0023]

[17] Khade P M, Scaramozzino D, Kumar A, et al. hdANM: a newcomprehensive dynamics model for protein hinges[J]. Biophysical Journal,2021, 120(22): 4955-4965.

[0024]

[18] Zhou L, Siegelbaum S A. Effects of surface water on proteindynamics studied by a novel coarse-grained normal mode approach[J].Biophysical Journal, 2008, 94(9): 3461-3474.

[0025]

[19] Kaynak B T, Findik D, Doruker P. RESPEC incorporates residuespecificity and the ligand effect into the elastic network model[J]. TheJournal of Physical Chemistry B, 2017, 122(21): 5347-5355.

[0026]

[20] Zhang Y, Zhang S, Xing J, et al. Normal mode analysis of membraneprotein dynamics using the vibrational subsystem analysis[J]. The Journal ofChemical Physics, 2021, 154(19): 195102.

[0027]

[21] Jia K, Kilinc M, Jernigan R L. Functional protein dynamicsdirectly from sequences[J]. The Journal of Physical Chemistry B, 2023, 127(9): 1914-1921.

[0028]

[22] Guarnera E, Berezovsky I N. Sequence-dependent model ofallosteric communication[J]. The Journal of Chemical Physics, 2023, 158(13):135101.

[0029]

[23] Na H, Song G. The performance of fine-grained and coarse-grainedelastic network models and its dependence on various factors[J]. Proteins: Structure, Function, and Bioinformatics, 2015, 83(7): 1273-1283.

[0030]

[24] Kawamoto S, Liu H, Miyazaki Y, et al. SPICA force field for proteins and peptides[J]. Journal of Chemical Theory and Computation, 2022,18(5): 3204-3217. Summary of the Invention

[0031] The purpose of this invention is to develop a high-performance improved Elastic Network Model (ENM) for efficiently and accurately characterizing protein dynamics, providing high-throughput reference information for studying protein conformation-function relationships. This method maintains (or even improves) computational accuracy comparable to all-atom ENM while retaining the efficiency of traditional ENM models, offering a high-throughput, high-precision, and low-cost technical path for protein conformation-function research. The key features of this invention are: ① The first-ever introduction of the Protein Joint Atom Modeling (SPICA) strategy into the traditional ENM computational framework. This modeling strategy utilizes experimentally obtained protein heavy atom coordinates for modeling, eliminating the need to add hydrogen atom coordinates, significantly reducing model complexity (number of atoms). ② Mapping corresponding physical interactions into the computational flow, matching SPICA force field parameters based on atom pair connection types (bonded, non-bonded, etc.) and atom types, accurately calculating elastic force constants, and constructing an elastic network topology reflecting reasonable physical interactions. Compared to the Cα modeling strategy used in traditional ENM, the joint atom modeling strategy has significant advantages in representing structural details; compared to the all-atom modeling strategy, this strategy significantly reduces the complexity of the model (total number of atoms) while maintaining detail and accuracy. ③ In this computational model, an optimized potential energy function adapted to ENM was used. Based on the characteristics of ENM, the interatomic interaction potential energy in the SPICA force field was specifically optimized to ensure the stability and accuracy of the computational model.

[0032] The united atom elastic network model (uaENM) described in this invention includes the following five core steps (as shown in the appendix). Figure 1 The process involves four steps: 1) Constructing a joint atomic model of the protein: The input protein's full atomic structure (from experimental or theoretical sources, requiring only the coordinates of heavy atoms in the model) is converted into a joint atomic model using tools provided in SPICA, and a topology is constructed. 2) Setting initial parameters and pairing force field parameters: Initial parameters such as the cutoff radius and non-bonded interaction range are set, the elastic network connecting atomic pairs are determined (spring connections are set), and the corresponding SPICA force field parameters are matched according to the connection type and atom type. 3) Calculating spring force constants: Using the matched SPICA force field parameters (bonding potential, van der Waals potential, etc.), the elastic force constants of the springs connecting each atomic pair are calculated by solving the second derivative of the corresponding potential energy function. 4) Constructing a Hessian matrix: Based on the established connection topology and the calculated force constants, a Hessian matrix (including contributions from two-body, three-body, and four-body potentials) mapping the protein's dynamic properties is constructed. Fifth, kinetic analysis and parameter optimization: Residue fluctuations are calculated based on the pseudo-inverse operation of the Hessian matrix, and their similarity is compared with experimental B-factors. The Pearson Correlation Coefficient (PCC) is then used to quantify these fluctuations. By maximizing this coefficient, the optimal parameters that best match the theoretical calculations with the experimental data are determined. Subsequently, key kinetic properties, including but not limited to residue fluctuations, inter-residue motion cross-correlation, and transfer entropy, are calculated based on these parameters.

[0033] First, the protein test set used to test the computational method (uaENM) proposed in this invention will be described:

[0034] This test uses the same protein test set (containing 135 protein structures) used to evaluate the performance of computational models employing the protein whole-atom modeling strategy (sbNMA). This dataset is available from [source missing]. https: / / github.com / gsongISU / sbNMAe To further comprehensively evaluate the performance of uaENM, this paper also constructed an additional test set containing six proteins with different topological features. All structures were obtained from the Protein Data Bank database. https: / / www.rcsb.orgSpecifically, these include: apolipoprotein kringle type IV 7 (PDB ID: 1I71), rich in cyclic secondary structures; the sixth PDZ domain of GRIP1 (PDB ID: 1N7E), with a complex and loose topology containing a mixture of α-helices, β-sheets, and cyclic regions; bleomycin resistance determinant (PDB ID: 1QTO), also with a complex topology, its secondary structure being a combination of α-helices, β-sheets, and loops; PedB (PDB ID: 2IP6), with a topology dominated by α-helices; the PKCiota kinase domain (PDB ID: 3A8X), with a relatively large number of residues; and linear diubiquitin (PDB ID: 3AXC), exhibiting a dimer conformation. The structures of the six proteins in the above test set are attached. Figure 2 As shown.

[0035] The uaENM calculation process consists of the following five steps:

[0036] Step 1: Construction of the Protein Joint Atom Model. Following the modeling strategy of SPICA, the protein joint atom input structure used in uaENM is generated from the heavy atom coordinates obtained through experimental analysis. The modeling details and topology construction rules are as follows:

[0037] 1. Construction of the joint atomic model

[0038] SPICA's joint atom modeling scheme achieves model simplification and residue property preservation through the following methods: ① Backbone (BB) simplification: The main chain atoms (Cα, C, O, N, and H) are simplified to an equivalent joint atom located at the original Cα position. ② Side chain (SC) hierarchical simplification: Based on chemical complexity and spatial configuration characteristics, the side chains are classified into four mapping structures (number of joint atoms 0-4). Examples: Glycine is represented by only BB atoms; alanine / valine / leucine adopts a BB + hydrophobic SC (Cβ / Cγ centroid) double joint atom configuration; charged and polar residues such as serine / threonine / aspartic acid are distinguished by charged and polar joint atoms, while retaining the polar characteristics of hydroxyl / carboxyl groups; aromatic residues (phenylalanine / tyrosine / tryptophan) are characterized by a three-joint atom ring topology to represent the benzene ring configuration, while retaining the joint atom representation of substituents. ③ Ionization site deployment: The charge of ionized residues such as arginine, lysine, and aspartic acid is located at designated side chain bonding atoms, with the N-terminal amino group and C-terminal carboxyl group counted as independent charge sites. In addition, the different residue types are also distinguished by the differences in the nonbonding potential parameters (Lennard-Jones (9-6)) of the different types of bonding atoms they contain.

[0039] 2. Rules for constructing topological structures (bonding interactions)

[0040] To maintain the stability of the joint atom model structure, the SPICA modeling strategy introduces the following interatomic bonding interactions: ① Bond constraint: acting between adjacent main chain / side chain joint atoms. ② Bond angle constraint: configured for three linearly adjacent atom groups (including BB and SC joint atoms). ③ Dihedral angle constraint: applied only to the phenylalanine / threonine side chains, fixing the spatial configuration of the benzene ring. ④ Secondary structure constraint: applying an elastic constraint (coefficient 1.195 kcal / mol / Ų) to the BB joint atomic spacing < 9.0 Å, ensuring the stability of the α-helix / β-sheet.

[0041] Step 2: Initial parameter setting and force field parameter pairing. uaENM, while adhering to the traditional ENM calculation framework, further introduces corresponding physical force fields. Before starting the calculation, the initial cutoff radius (usually set to 5 Å) and the non-bonded interaction range parameter (n) need to be preset. rc(Usually set to 0.1), and corresponding parameters are matched for the interaction potential function between each pair of bonded atoms. According to the construction rules of the protein bonded atom model in step one, the interactions between bonded atoms can be divided into two categories: bonding interactions and non-bonding interactions. The model contains 28 different types of bonded atoms, and the non-bonding interaction potential parameters between them have been provided in the SPICA force field. The bonding interactions in the protein bonded atom model are further divided into three categories: bond stretching, bond angle bending, and dihedral angle torsion. Referring to the parameters provided by the SPICA force field, the potential constants for these three types of bonding interactions are set to 50 kcal / mol / Ų, 20 kcal / mol / Ų, and 50 kcal / mol / Ų, respectively, in this invention.

[0042] Step 3: Calculation of the spring force constant. The semi-empirical interaction potential energy between atoms cannot usually be directly applied to the uaENM calculation framework; it needs to be converted into the form of second-order resonant potential energy. This invention achieves this conversion by calculating the second derivative of the semi-empirical interaction potential energy to obtain the energy field parameters of the spring resonant potential. The detailed calculation process is as follows: The common form of the non-bonded van der Waals interaction potential energy between atoms is given by Equation 1: (1) Where ε represents the potential well depth, r0 is the equilibrium distance, both provided by the SPICA force field; r represents the actual distance between atoms. Its second derivative is shown in Equation 2: (2) The van der Waals potential nonbonded elastic constant k derived from Equation 2 l (vdW) may contain negative values. Similar to traditional ENM-type resonant potential energy calculation models, negative elastic constants can create energy saddle points in uaENM calculations, leading to protein structural instability. To remove k... l For negative values ​​in (vdW), this invention uses the reduced form of Formula 2 to calculate the spring force constant, as detailed in Formula 3: (3) in, (4) n rc This is a parameter used to control the range of the van der Waals potential energy. A typical form of the non-bonded electrostatic interaction potential energy between atoms is shown in Equation 5: (5) The second derivative of this potential energy is: (6) Where q i and q jThe magnitudes of the charges carried by the joint atoms i and j are provided in the SPICA force field, respectively.

[0043] Since the spring force constant obtained from Equation 6 still has a small number of negative values, this will lead to a decrease in the non-bonded spring force constant k. nonbonding (defined as k) l (elec) and k l Negative values ​​also appear in the sum of (vdW); to ensure the stability of the computational model, these negative values ​​are set to zero. Furthermore, although the atomic density of the protein joint atom model is relatively low, experimental errors may cause some joint atoms to be too close together, resulting in extremely large k values. nonbonding Numerical values. Theoretically, non-bonded interactions between atoms should be weaker than the upper limit of hydrogen bonding interactions (14.29 kcal / mol / Å). 2 Therefore, non-bonded interaction values ​​exceeding this threshold are adjusted to the upper limit of hydrogen bond interaction values. For bonded interactions between joint atoms (including bond stretching, bond angle bending, and dihedral torsion), the values ​​in step 2 (50 kcal / mol / Ų, 20 kcal / mol / Ų, and 50 kcal / mol / Ų) are directly used as elastic force constants.

[0044] Step 4: Constructing the Hessian matrix. Similar to ENM, the method for calculating the elements of the Hessian matrix constituting uaENM is to solve for the second derivative of the resonant potential energy. When i ≠ j: (7) When i = j: (8) The uaENM integrates interatomic interactions, which can be categorized into two-body, three-body, and four-body potentials based on the number of participating atoms. Therefore, the Hessian matrix of the uaENM needs to be constructed sequentially using the two-body, three-body, and four-body potentials. The specific calculation details are as follows:

[0045] 1. Two-body pose (V) l Hessian matrix construction

[0046] The calculation method for the second-order partial derivatives of the elements in the Hessian matrix of a two-body potential is consistent across all directions. First, define the Cartesian coordinates of atoms i and j as (X... i , Y i Z i ) and (X j , Y j Z j This invention uses the sub-elements h that constitute the Hessian matrix.ij Taking the X and XY directions as examples, the calculation method is demonstrated through the following formula: (9) (10) r is the distance between atoms i and j, k l Let represent the elastic coefficient of the two-body potential interaction between atomic pairs i and j, including bonding (set as a constant, 50 kcal / mol / Ų) and non-bonding terms (Equations 3 and 6). Diagonal element h ii The calculation method is as follows: (11) (12)

[0047] 2. Three-body pose (V) θ Hessian matrix construction

[0048] The bond angles formed by the three atoms i, j, and k (see appendix) Figure 3 The process of constructing the Hessian matrix of a three-body potential is illustrated using an example. (See Appendix...) Figure 3 In the diagram, coordinate vectors r represent atoms i, j, and k, respectively. i r j and r k Then attach Figure 3 The vectors A, B, and G used for auxiliary calculations can be represented as: (13) (14) (15) The sine and cosine values ​​of the key angle can be expressed as: (16) (17) The first derivatives of the key angle θ with respect to vectors A and B, and the first derivatives of vectors A and B with respect to vector r i r j and r k The first derivatives can be expressed as: (18) (19) (20) (twenty one) (twenty two) (twenty three) The second derivative of the bond angle vibration potential is shown below: (twenty four) The protein structure input to uaENM is generally considered to be its native conformation; therefore, the internal torque of the protein represented by the second term in Equation 24 should be zero. (Equation 24) The term represents the bond angle vibrational potential energy parameter, which is set to a constant (20 kcal / mol / Ų) based on the parameters provided in the SPICA force field. In Equation 24... The term is calculated using the following formula: (25) in, (26) Substitute the results of Equations 18-23 into Equation 26, and then substitute the results of Equation 26 into Equation 24 to solve for the elements of the Hessian corresponding to the three-body potential.

[0049] 3. Four-body posture (V) φ Hessian matrix construction

[0050] The construction process is basically the same as that of the Hessian potential of the three-body potential. The bond angles formed by the four atoms i, j, k, and l (see appendix) Figure 4 The process of constructing the Hessian matrix of a four-body potential is illustrated using an example. (See Appendix...) Figure 4 In the diagram, coordinate vectors r represent atoms i, j, k, and l, respectively. i r j r k and r l Appendix Figure 4 The vectors F, G, and H used for auxiliary calculations can be represented as: (27) (28) (29) Appendix Figure 4 Plane vectors A and B in the diagram can be defined as follows: (30) (31) The dihedral angle φ formed by four atoms can be defined as the angle between two planes, and its sine and cosine values ​​can be written as: (32) (33) The first derivatives of the dihedral angle φ with respect to vectors F, G, and H, and the first derivatives of vectors F, G, and H with respect to vector r i r j r k , and r l The first derivatives can be expressed as: (34) (35) (36) (37) (38) (39) (40) (41) (42) The second derivative of the dihedral torsional potential is shown below: (43) Similar to Equation 24, the second term in Equation 43 is zero. In Equation 43... The term represents the dihedral torsional potential energy parameter, which is set to a constant (50 kcal / mol / Ų) based on the parameters provided in the SPICA force field. In Equation 43... The term is calculated using the following formula: (44) in, (45) Similarly, substituting the results of equations 34-42 into equation 45, and then substituting the result of equation 45 into equation 43, allows us to solve for the elements of the Hessian matrix corresponding to the four-body potential. Finally, by adding the Hessian matrices of the two-body, three-body, and four-body potentials, we obtain the Hessian matrix that fully describes protein dynamics for subsequent calculations.

[0051] Step 5: Kinetic Analysis and Parameter Optimization. Protein kinetic properties in uaENM all require calculation using the Hessian pseudo-inverse matrix. (46) Where λ k and u k These are the eigenvalues ​​and eigenvectors corresponding to the k-th motion pattern, respectively. Protein residue fluctuations and their cross-correlation are defined as follows: (47) (48) Where k B γ is the Boltzmann constant, T is the thermodynamic temperature, and γ is the elastic constant, which is already embedded in the Hessian matrix. The residue motion cross-correlation C ij_ENM Calculated using the following formula: (49) Among them, when C ij_ENM As the value approaches 1, the residues move in the same direction; as it approaches -1, the residues move in opposite directions; as it approaches 0, the residue movements are unrelated. (From C) ij_ENM The matrix constructs a dynamic cross-correlation map (DCCM) to quantitatively characterize the collective motion properties of residues, providing a basis for large-scale conformational change analysis.

[0052] To verify the computational accuracy of the model, this patent uses the Pearson correlation coefficient (PCC) to quantify the similarity between the theoretical value and the reference value calculated by uaENM. The calculation method for the PCC (B_PCC) between the theoretical and experimental B factors calculated by uaENM is as follows: (50) in, and The B factor and its mean value are calculated for ENMs, respectively. and These represent the experimentally obtained B factor and its average value. The method for calculating the theoretical B factor using uaENM is as follows: (51)

[0053] The calculation method for the PCC (DCCM_PCC) between the theoretical DCCM obtained through uaENM and the reference DCCM obtained through molecular dynamics simulation is as follows: (52) in, and The DCCM calculated for ENMs and its average value are shown below. and These are the DCCM calculated by MD simulation and its average value, respectively.

[0054] Typically, initial parameters cannot yield optimal B_PCC values. This invention addresses this by systematically scanning the cutoff radius (cutoff: 5.0-20.0 Å, step size: 1.0 Å) and the non-bonded interaction range parameter (n). rcThe parameter combination is optimized using a range of 0.1-4.0 (step size: 0.1) to ensure the accuracy of subsequent dynamic property calculations.

[0055] In this invention, we propose a joint atomic elastic network model (uaENM) that characterizes the dynamic properties of proteins by employing a joint atomic modeling strategy and corresponding physical force fields. This computational model aims to improve computational performance by using a refined protein model with matching physical potentials, while reducing sensitivity to the quality of the input protein structure. These advantages are essentially attributed to the following: uaENM uses a simplified joint atomic model, incorporates corresponding coarse-grained physical interactions, and significantly improves the stability and integrity of the constructed topology by optimizing the potential energy function of the spring constant. These improvements allow uaENM to effectively overcome the limitations of high-resolution ENMs, such as high sensitivity to input structures and insufficient computational efficiency, while preserving the refined protein topology. Attached Figure Description

[0056] Appendix Figure 1 Procedure for calculating protein dynamics for uaENM;

[0057] Appendix Figure 2 The structure of six proteins was tested;

[0058] Appendix Figure 3 A schematic diagram of the bond angles formed by three atoms;

[0059] Appendix Figure 4 A schematic diagram of a dihedral angle formed by four atoms;

[0060] Appendix Figure 5 Fluctuations in protein 1I71 residues calculated and experimentally measured using uaENM;

[0061] Appendix Figure 6 Comparison of B_PCC for uaENM, sbNMA and ANM;

[0062] Appendix Figure 7 A comparison of DCCM_PCC for uaENM, sbNMA, and ANM;

[0063] Appendix Figure 8 Comparison of computational efficiency between uaENM, sbNMA, and ANM.

[0064] The present invention will be further described below with reference to the embodiments, but the present invention is not limited to the following embodiments.

[0065] Example 1

[0066] The following section uses a protein with PDB ID 1I71 as an example to introduce model parameter optimization and residue fluctuation calculation, and provides a comparison between the calculation results and the residue fluctuations (B factor) obtained experimentally.

[0067] 1. Preparation

[0068] Before executing the uaENM calculation program, a Python 3.7 environment must be configured on the Linux system, and MATLAB 2024, the spica-tools tool (provided by SPICA Force Field), and the DSSP tool must be correctly installed. Download instructions and detailed installation guides for these three software programs are available on their respective official websites. It is important to note that the spica-tools tool has ambiguity regarding atom labeling when processing isoleucine (ILE) residues; the CD1 atom in the protein coordinate file must be renamed to CD manually or by writing a script.

[0069] 2. Model Transformation and Topology Construction

[0070] First, input the protein experimental crystal structure file (without adding hydrogen atoms) into the spica-tools tool to generate a joint atomic model file of the protein; then, input the joint atomic model file and the original experimental crystal structure file into the spica-tools tool again to generate a topology file of the joint atomic model.

[0071] 3. Characterization of protein kinetic properties

[0072] (1) Read the protein joint atomic model and match the force field parameters

[0073] The program "res_type_spica.m" is used to provide SPICA force field parameters, the program "read_top.m" reads the topology file data generated by spica-tools, and the program "pdbread_cg_allatoms.m" integrates the outputs of "read_top.m" and "res_type_spica.m" to match the corresponding force field parameters for each joint atom in the model.

[0074] (2) Construction of Hessian matrix

[0075] The programs “hessian_2.m”, “hessian_3.m”, and “hessian_4.m” are used to generate Hessian matrices for two-body, three-body, and four-body potentials, respectively.

[0076] (3) Solving the Hessian pseudo-inverse matrix and calculating residue fluctuations

[0077] The program "calculate_uaENM" is responsible for integrating the two-body, three-body, and four-body Hessian matrices calculated by "hessian_2.m", "hessian_3.m", and "hessian_4.m" into a single complete Hessian matrix. Subsequently, this Hessian matrix is ​​orthogonally decomposed to obtain eigenvalues ​​and their corresponding eigenvectors arranged in ascending order. Finally, using these eigenvectors and eigenvalues, residue fluctuations (considering three directions) containing k motion modes are calculated.

[0078] (4) Finding the optimal parameters

[0079] The program “select_parameters.m” systematically inputs the cutoff radius (range 5.0-20.0 Å, step size 1.0 Å) and non-bonded interaction range parameter (range 0.1-4.0, step size 0.1) into the program “calculate_uaENM” to calculate the Pearson correlation coefficient between the theoretical and experimental values ​​of residue fluctuations. The optimal parameters for the protein are determined by selecting the parameter combination that maximizes the Pearson correlation coefficient.

[0080] The final visualization results are shown in the appendix. Figure 5 middle.

[0081] Example 2

[0082] This invention compares the computational performance of the constructed uaENM (using a joint atom modeling strategy), traditional ANM (using a coarse-grained modeling strategy based on Cα atoms), and sbNMA (using a full atom modeling strategy) on two test sets (containing 135 and 6 protein structures, respectively). Some computational procedures are the same as in Example 1; only the computational steps not detailed in Example 1 and the results of this example are described here.

[0083] To comprehensively evaluate the computational performance of the proposed method, uaENM, ANM, and sbNMA were systematically tested on two test sets in terms of residue fluctuation (B factor) prediction, motion cross-correlation, protein collective motion capture ability, and computational efficiency.

[0084] On a test set containing 135 protein structures, the performance of three models in predicting residue fluctuations was evaluated. The B-factor prediction correlation coefficient (B_PCC, Equation 50) was used as the quantitative indicator. The results are shown in the appendix. Figure 6 And Table 1.

[0085] Table 1. Comparison of ENMs' predicted residue fluctuation properties at 135 protein structures.

[0086] Table 1 shows the average (standard deviation) of the B_PCC values ​​for 135 proteins calculated using the three methods, as well as the range of optimal cutoff values ​​and the total number of springs defined by each method. (From Table 1 and Appendix...) Figure 6 It can be seen that uaENM outperforms sbNMA and ANM in predicting residue fluctuations. Specifically, compared with sbNMA and ANM, the average B_PCC calculated by uaENM is improved by 18.18% and 12.07%, respectively. Figure 8 Furthermore, uaENM achieved higher B_PCC in a wider range of test samples, with improvements of 85.93% and 74.81%, respectively. Additionally, uaENM reduced the total number of springs defined during calculation by up to 94.95% compared to sbNMA. These test results demonstrate that uaENM has achieved its goal of maintaining high computational performance with fewer springs.

[0087] The ability of three models to capture residue motion cross-correlation was evaluated on a test set containing six protein structures. Performance was quantified using the Pearson correlation coefficient (DCCM_PCC) by comparing the dynamic cross-correlation matrix (DCCM, Equation 49) calculated by each model with the DCCM obtained from MD simulations. The results are shown in the appendix. Figure 7 As shown. (Attached) Figure 7 The results show that uaENM's DCCM_PCC is higher than that of traditional ANM, and it maintains a high correlation even after including mid-to-high frequency motion modes. Compared with sbNMA, which uses a full-atom protein model, uaENM achieves comparable DCCM computational accuracy with a more simplified atomic model.

[0088] The ability of the three models to capture functional slow motion patterns of proteins was further evaluated on a test set containing six protein structures. The evaluation was measured by comparing the overlap between the low-frequency motion patterns (low-frequency feature vectors) calculated by ENMs and the principal components (PCs) obtained by MD simulations. Specifically, the maximum overlap between the top 10% of low-frequency motion patterns in ENMs and PC1, PC4, PC7, and PC10 (denoted as ) was used. , , and ), calculate according to formulas 53 and 54 respectively: (53) (54) Where, p i and u iLet p represent the principal component eigenvectors of the MD simulation and the motion mode eigenvectors of the ENMs, respectively. i ‖and‖u i ‖ represents its modulus length.

[0089] The root mean square inner product (RMSIP) is used to evaluate the overall similarity between the top l principal component eigenvectors of the MD simulation and the top 10% low-frequency motion patterns of the ENMs. The calculation formula is as follows: (55) like A value close to 1 for RMSIP indicates that the motion patterns of ENMs are highly consistent with the principal component directions of MD simulations; a value close to 0 indicates no overlap. The results are shown in Table 2. Table 2 shows that uaENM... It significantly outperforms ANM and sbNMA in both RMSIP and other aspects, indicating its superior performance in capturing low-frequency protein motion patterns related to biological functions.

[0090] Table 2. Average overlap and RMSIP of motion patterns obtained by ENMs and MD simulations on six protein datasets.

[0091] The bolded values ​​represent the best performance values ​​for each metric.

[0092] Finally, the computational efficiency of uaENM was compared and tested. The time required for the three methods to process protein structures with tens to thousands of residues was statistically analyzed, and the results are attached. Figure 8 As shown. (Attached) Figure 8 The results show that as the number of residues increases, the computation time and Hessian matrix size of sbNMA based on the all-atom model grow exponentially, while uaENM using the joint atom model is slightly longer than the traditional ANM, and its computational efficiency is significantly better than sbNMA.

[0093] In summary, this example comprehensively tested uaENM's capabilities and computational efficiency in residue fluctuation prediction, motion correlation prediction, and capturing low-frequency motion patterns of protein functionality. The results show that uaENM's joint atom modeling strategy not only inherits the accuracy advantages brought by introducing refined protein models and corresponding physical interactions, but also significantly reduces the quality requirements for the input protein structure, completely avoids interference that may be introduced by artificially added hydrogen atoms, and maintains computational efficiency comparable to traditional ANM.

Claims

1. An improved method for calculating elastic network models based on a joint atom modeling strategy for proteins and a real physical force field, referred to as the united atom elastic network model (uaENM), is used for high-throughput and high-precision characterization of the dynamic properties of proteins, characterized by... This study employs a protein joint atom modeling strategy, incorporating corresponding interatomic interaction potentials. This approach improves computational accuracy while offering advantages such as low requirements on input structure quality and high computational efficiency. The strategy comprises five core steps:

1. Constructing the protein joint atom model: The input protein full-atom model (experimental or theoretical) is converted into a joint atom model using tools provided in the SPICA force field, and a topology is constructed.

2. Initial parameter setting and force field parameter pairing: Initial cutoff radius and non-bonded interaction range parameters are set, elastic network connecting atom pairs are determined (spring connections are established), and corresponding SPICA force field parameters are matched based on the atom pair connection type and atom type.

3. Calculating spring force constants: Using the matched SPICA force field parameters (bonding potential, van der Waals potential, etc.), the elastic force constants of the springs connecting each atom pair are calculated by solving the second derivative of the corresponding potential energy function.

4. Constructing the Hessian matrix: Based on the established connection topology and the calculated force constants, a Hessian matrix (including contributions from two-body, three-body, and four-body potentials) mapping the dynamic properties of the protein system is constructed. Fifth, kinetic analysis and parameter optimization: Residue fluctuations are calculated based on the pseudo-inverse operation of the Hessian matrix, compared with the experimental B factor, and the Pearson correlation coefficient is calculated. By maximizing this coefficient, the optimal parameters that best match the theoretical calculations with the experimental data are determined. Subsequently, key kinetic properties, including but not limited to residue fluctuations, inter-residue motion cross-correlation, and transfer entropy, are calculated based on these parameters.

2. The joint atomic elastic network model based on protein joint atomic modeling strategy and real physical force field as described in claim 1 is characterized in that: the calculation steps enable uaENM to accurately and quickly obtain protein dynamic properties. Step 1: Construction of the Protein Joint Atom Model. Following the modeling strategy of SPICA, the protein joint atom input structure used in uaENM is generated from the heavy atom coordinates obtained through experimental analysis. The modeling details and topology construction rules are as follows: Construction of the Joint Atom Model: SPICA's joint atom modeling scheme achieves model simplification and residue property preservation through the following methods: Backbone (BB) simplification: Residue backbone atoms (Cα, C, O, N, H) are simplified to an equivalent joint atom located at the original Cα position. Side Chain (SC) hierarchical simplification: Based on chemical complexity and spatial configuration characteristics, residue side chains are classified into four mapping structures (number of joint atoms 0-4). Examples: Glycine is represented by only BB atoms; alanine / valine / leucine adopts a BB + hydrophobic SC (Cβ / Cγ centroid) double joint atom configuration; charged and polar residues such as serine / threonine / aspartic acid are distinguished by charged and polar joint atoms, preserving the polar characteristics of hydroxyl / carboxyl groups; Aromatic residues (phenylalanine / tyrosine / tryptophan) are characterized by a three-unit ring topology, retaining the representation of the unit atom of substituents. Ionization site deployment: the charge of ionized residues such as arginine / lysine / aspartic acid is located at designated side-chain unit atoms, with the N-terminal amino group and C-terminal carboxyl group counted as independent charge sites. Furthermore, the various residue types are distinguished by differences in the nonbonding potential parameters (Lennard-Jones (9-6)) of the different types of unit atoms they contain. Topological structure (bonding interaction) construction rules: To maintain the stability of the joint atom model structure, the SPICA modeling strategy introduces the following interatomic bonding interactions: Bond constraint: Acts between adjacent main chain / side chain joint atoms. Bond angle constraint: Configured for three linearly adjacent atom groups (including BB and SC joint atoms). Dihedral angle constraint: Applied only to the phenylalanine / threonine side chain, fixing the benzene ring spatial configuration. Secondary structure constraint: Applying elastic constraints (coefficient 1.195 kcal / mol / Ų) to the BB joint atomic spacing < 9 Å, ensuring the stability of α-helix / β-sheet. Step 2: Initial parameter setting and force field parameter pairing. uaENM, while using the traditional Elastic Network Model (ENM) computational framework, further introduces real interactions between atoms. Before starting the calculation, the initial cutoff radius (usually set to 5 Å) and the non-bonded interaction range (n) need to be preset. rc (usually set to 0.1), and corresponding parameters are matched for the interaction potential function between each pair of joint atoms. According to the construction rules of the protein-bonded atom model in Step 1, the interactions between bonded atoms can be divided into two categories: bonding interactions and non-bonding interactions. The model contains 28 different types of bonded atoms, and the potential energy parameters for their non-bonding interactions have been provided in the SPICA force field. The bonding interactions in the protein-bonded atom model are further divided into three categories: bond stretching, bond angle bending, and dihedral angle torsion. Referring to the parameters provided by the SPICA force field, the potential energy constants for these three types of bonding interactions are set to 50 kcal / mol / Ų, 20 kcal / mol / Ų, and 50 kcal / mol / Ų, respectively, in this invention. Step 3: Calculation of spring force constant. The semi-empirical interaction potential energy between atoms cannot usually be directly applied to the uaENM calculation framework and needs to be converted into the form of second-order resonant potential energy. This invention obtains the energy field parameters of the spring resonant potential by calculating the second derivative of the semi-empirical interaction potential energy, thereby achieving the transformation. The detailed calculation process is as follows: The common form of the non-bonded van der Waals interaction potential energy between atoms is given by Equation 1: (1) Where ε represents the potential well depth, r0 is the equilibrium distance, both provided by the SPICA force field; r represents the actual distance between atoms. Its second derivative is shown in Equation 2: (2) The van der Waals potential nonbonded elastic constant k derived from Equation 2 l (vdW) may contain negative values. Similar to traditional ENM-type resonant potential energy calculation models, negative elastic constants can create energy saddle points in uaENM calculations, leading to protein structural instability. To remove k... l For negative values ​​in (vdW), this invention uses the reduced form of Formula 2 to calculate the spring force constant, as detailed in Formula 3: (3) in, (4) n rc This is a parameter used to control the range of the van der Waals potential energy. A typical form of the non-bonded electrostatic interaction potential energy between atoms is shown in Equation 5: (5) The second derivative of this potential energy is: (6) Where q i and q j The magnitudes of the charges carried by the joint atoms i and j are provided in the SPICA force field, respectively. Since the spring force constant obtained from Equation 6 still has a small number of negative values, this will lead to a decrease in the non-bonded spring force constant k. nonbonding (defined as k) l (elec) and k l Negative values ​​also appear in the sum of (vdW); to ensure the stability of the computational model, these negative values ​​are set to zero. Furthermore, although the atomic density of the protein joint atom model is relatively low, experimental errors can still cause some joint atoms to be too close together, resulting in extremely large k values. nonbonding Numerical values. Theoretically, non-bonded interactions between atoms should be weaker than the upper limit of hydrogen bonding interactions (14.29 kcal / mol / Å). 2 Therefore, non-bonded interaction values ​​exceeding this threshold are adjusted to the upper limit of hydrogen bond interaction values. For bonded interactions between joint atoms (including bond stretching, bond angle bending, and dihedral torsion), the values ​​in step 2 (50 kcal / mol / Ų, 20 kcal / mol / Ų, and 50 kcal / mol / Ų) are directly used as elastic force constants. Step 4: Constructing the Hessian matrix. Similar to ENM, the method for calculating the elements of the Hessian matrix constituting uaENM is to solve for the second derivative of the resonant potential energy. When i ≠ j: (7) When i = j: (8) uaENM integrates real interatomic interactions, which can be categorized into two-body, three-body, and four-body potentials based on the number of participating atoms. Therefore, the Hessian matrix of uaENM needs to be constructed sequentially using the two-body, three-body, and four-body potentials, with specific calculation details as follows: Two-body posture (V) l Construction of the Hessian matrix: The calculation method for the second-order partial derivatives of the elements in the Hessian matrix of a two-body potential is consistent across all directions. First, define the Cartesian coordinates of atoms i and j as (X... i Y i Z i ) and (X j ,Y j Z j This invention uses the sub-elements h that constitute the Hessian matrix. ij Taking the X and XY directions as examples, the calculation method is demonstrated through the following formula: (9) (10) r is the distance between atoms i and j, k l Let represent the elastic coefficient of the two-body potential interaction between atomic pairs i and j, including bonding (set as a constant, 50 kcal / mol / Ų) and non-bonding terms (Equations 3 and 6). Diagonal element h ii The calculation method is as follows: (11) (12) Three-body pose (V) θ Hessian matrix construction: Taking the bond angle formed by three atoms i, j, and k (Figure 3) as an example, the process of constructing the Hessian matrix of the three-body potential is explained. In Figure 3, coordinate vectors r representing atoms i, j, and k are defined respectively. i r j and r k Then, the vectors A, B, and G used for auxiliary calculation in Figure 3 can be represented as: (13) (14) (15) The sine and cosine values ​​of the key angle can be expressed as: (16) (17) The first derivatives of the key angle θ with respect to vectors A and B, and the first derivatives of vectors A and B with respect to vector r i r j and r k The first derivatives can be expressed as: (18) (19) (20) (21) (22) (23) The second derivative of the bond angle vibration potential is shown below: (24) The protein structure input to uaENM is generally considered to be its native conformation; therefore, the internal torque of the protein represented by the second term in Equation 24 should be zero. (Equation 24) The term represents the bond angle vibrational potential energy parameter, which is set to a constant (20 kcal / mol / Ų) based on the parameters provided in the SPICA force field. In Equation 24... The term is calculated using the following formula: (25) in, (26) Substitute the results of Equations 18-23 into Equation 26, and then substitute the calculation results of Equation 26 into Equation 24 to solve for the elements of the Hessian corresponding to the three-body potential. Four-body posture (V) φ Hessian matrix construction: The process is basically the same as that for the three-body potential. Taking the bond angle formed by four atoms i, j, k, and l (Figure 4) as an example, the Hessian matrix construction process for the four-body potential is explained. In Figure 4, coordinate vectors r representing atoms i, j, k, and l are defined respectively. i r j r k and r l The vectors F, G, and H used for auxiliary calculations in Figure 4 can be represented as: (27) (28) (29) The planar vectors A and B in Figure 4 can be defined as follows: (30) (31) The dihedral angle φ formed by four atoms can be defined as the angle between two planes, and its sine and cosine values ​​can be written as: (32) (33) The first derivatives of the dihedral angle φ with respect to vectors F, G, and H, and the first derivatives of vectors F, G, and H with respect to vector r i r j r k , and r l The first derivatives can be expressed as: (34) (35) (36) (37) (38) (39) (40) (41) (42) The second derivative of the dihedral torsional potential is shown below: (43) Similar to Equation 24, the second term in Equation 43 is zero. In Equation 43... The term represents the dihedral torsional potential energy parameter, which is set to a constant (50 kcal / mol / Ų) based on the parameters provided in the SPICA force field. In Equation 43... The term is calculated using the following formula: (44) in, (45) Similarly, substituting the results of equations 34-42 into equation 45, and then substituting the result of equation 45 into equation 43, allows us to solve for the elements of the Hessian matrix corresponding to the four-body potential. Finally, by adding the Hessian matrices of the two-body, three-body, and four-body potentials, we obtain the Hessian matrix that fully describes protein dynamics for subsequent calculations. Step 5: Kinetic Analysis and Parameter Optimization. Protein kinetic properties in uaENM all require calculation using the Hessian pseudo-inverse matrix. (46) Where λ k and u k These are the eigenvalues ​​and eigenvectors corresponding to the k-th motion pattern, respectively. Protein residue fluctuations and their cross-correlation coefficients are defined as follows: (47) (48) Where k B γ is the Boltzmann constant, T is the thermodynamic temperature, and γ is the elastic constant, which is already embedded in the Hessian matrix. The residue motion cross-correlation C ij_ENM Calculated using the following formula: (49) Where C ij_ENM As the value approaches 1, the residues move in the same direction; as it approaches -1, the residues move in opposite directions; as it approaches 0, the residue movements are unrelated. (From C) ij_ENM The matrix constructs a dynamic cross-correlation map (DCCM) to quantitatively characterize the collective motion properties of residues, providing a basis for large-scale conformational change analysis. To verify the computational accuracy of the model, this patent uses Pearson's Correlation Coefficient (PCC) to quantify the similarity between the theoretical value calculated by uaENM and the reference value. The calculation method for the PCC (B_PCC) between the theoretical and experimental B factors calculated by uaENM is as follows: (50) in, and The B factor and its mean value are calculated for ENMs, respectively. and These represent the experimentally obtained B factor and its average value. The method for calculating the theoretical B factor using uaENM is as follows: (51) The calculation method for the PCC (DCCM_PCC) between the theoretical DCCM obtained through uaENM and the reference DCCM obtained through molecular dynamics simulation is as follows: (52) in, and The DCCM calculated for ENMs and its average value are shown below. and These are the DCCM calculated by MD simulation and its average value, respectively. Typically, initial parameters cannot yield optimal B_PCC values. This invention addresses this by systematically scanning the cutoff radius (cutoff: 5.0-20.0 Å, step size: 1.0 Å) and the non-bonded interaction range parameter (n). rc The parameter combination is optimized using a range of 0.1-4.0 (step size: 0.1) to ensure the accuracy of subsequent dynamic property calculations.

3. The joint atomic elastic network model based on protein joint atomic modeling strategy and real physical force field as described in claim 1 is characterized by: for the first time introducing the protein joint atomic modeling strategy (based on the SPICA force field) into the computational framework of the traditional elastic network model (ENM). The resolution of this protein modeling strategy is between the all-atomic modeling strategy used in traditional NMA and the coarse-grained modeling strategy based on Cα atoms in traditional ENM, thus balancing the computational accuracy brought by the all-atomic modeling strategy and the computational efficiency brought by the Cα atom coarse-grained modeling strategy.

4. The joint atomic elastic network model based on protein joint atomic modeling strategy and real physical force field as described in claim 1 is characterized in that: this study introduces complex joint atomic force field parameters (SPICA force field parameters) into the traditional elastic network model (ENM) for the first time, and performs targeted optimization on the potential energy function in the SPICA force field to improve its adaptability in the traditional ENM calculation framework, thereby significantly enhancing the stability of uaENM calculation.

Images