Variable-temperature infrared spectrum calculation method for protein biotoxins

Abstract

The invention provides a variable-temperature infrared spectrum calculation method for protein biotoxins, which comprises the following steps of: putting a macromolecular bioprotein structure into a box with a boundary, which is established according to an actual chemical environment, continuously and randomly vibrating molecules in the box by using a molecular dynamics method, and continuously and randomly vibrating the molecules in the box during the period, and the system records the change of the total dipole moment in the box at each moment, and then converts a dipole moment signal into infrared frequency and intensity through Fourier transform, so that the infrared spectrum of the macromolecular bioprotein is obtained, and accurate analysis and rapid detection of the macromolecular biotoxin are realized. The difficulty that spectrum calculation of a large molecular weight system cannot be carried out due to the serious scale problem of a first principle and a density functional theory method is overcome, and a complete route of infrared spectrum dynamic properties of the system is obtained by adopting a molecular dynamics calculation method; therefore, infrared spectrum prediction of large-scale biological protein molecules is realized through a calculation means.

Description

Technical Field

[0001] The present invention belongs to the technical field of biochemistry, and particularly relates to a variable-temperature infrared spectrum calculation method for protein biotoxins based on molecular dynamics simulation. Background Art

[0002] The groups of biotoxin organic molecules have infrared fingerprint characteristics, so infrared spectroscopy is often used to analyze the composition and changes of biotoxins. In the field of biosafety, the analysis and monitoring of biotoxins are of great significance. However, due to their large molecular weight and complex structure, the infrared spectral signals of macromolecular biotoxins are weak and easily affected by sample preparation and environmental factors, increasing the difficulty of spectral simulation. In addition, macromolecular biotoxins may contain multiple different chemical groups, and these groups may produce overlapping vibration absorption peaks in the infrared spectrum, making it complicated to accurately extract the information of each group from the spectrum. Therefore, in order to improve the application effect of infrared spectroscopy in the analysis of macromolecular biotoxins, it is necessary to improve the spectral analysis algorithm, and there is an urgent need to develop an effective method for calculating the infrared spectrum of biotoxins in an ultra-large protein system to achieve the precise analysis and rapid detection of macromolecular biotoxins.

[0003] At present, the infrared spectral simulation research on biotoxins is limited to small molecule types. For more complex macromolecular biotoxins, the existing technical means are still insufficient. The methods for simulating the infrared spectra of biomolecules mainly include the Density Functional Theory (DFT) and the Molecular Models (MM) method. DFT is based on quantum mechanics and can provide accurate molecular structures and vibration information. However, the computational cost is relatively high. For large biomolecules, the complexity and scale of the simulation require DFT to have a large amount of computing resources and time, which cannot be satisfied currently. In addition, the DFT method has certain challenges in dealing with non-local interactions and weak interactions, which may lead to errors or deficiencies in simulating macromolecular systems. The MM method is faster and suitable for large systems, but the description of flexible molecules is relatively rough. For ultra-large systems (thousands of atoms) of protein substances, first-principles calculations cannot be used, and this calculation method has three main drawbacks: First, the harmonic approximation will bring relatively large errors; second, there are no intermolecular interactions in gas-phase single molecules, ignoring the molecular interactions between substances; third, only relatively small molecular systems, that is, only containing dozens to hundreds of atoms, can be considered. Therefore, effective algorithms are needed to solve the corresponding equations under various theoretical models. Molecular Dynamics (MD) simulation relying on molecular force fields is the most commonly used tool for studying the dynamic properties of biomacromolecules at present. Compared with strict quantum chemical calculations, the computational efficiency of molecular force fields is much higher. Therefore, molecular dynamics (MD) simulation is the most important tool for studying the structures and dynamic properties of biomacromolecules at present. Using molecular dynamics simulation to calculate the infrared spectra of condensed-phase molecules not only includes intermolecular interactions but also considers anharmonic effects, and can also be applied to large systems with thousands of atoms, which is a relatively reasonable calculation method for the infrared spectra of biomacromolecules. Therefore, it is extremely necessary to design a variable-temperature infrared spectral calculation method for protein biotoxins based on molecular dynamics simulation to calculate the infrared spectra of ultra-large protein system biotoxins. Summary of the Invention

[0004] (1) Technical Problems to be Solved

[0005] The present invention provides a variable-temperature infrared spectral calculation method for protein biotoxins to solve the technical problem of how to overcome the limitations in the research on the infrared spectral calculation of biomacromolecules.

[0006] (2) Technical Solutions

[0007] To solve the above technical problems, the present invention provides a variable-temperature infrared spectral calculation method for protein biotoxins, and the variable-temperature infrared spectral calculation method includes the following steps:

[0008] S1. Establish the complete structure of the biotoxin protein

[0009] Obtain the initial structure of the biotoxin protein from the protein data bank website, and use DiscoveryStudio software to complete the hydrogen atoms to obtain the complete protein structure;

[0010] S2. Place the complete structure of the biotoxin protein into a box with periodic boundary conditions established according to the actual chemical environment. This box is a space that accommodates the protein and the solvent environment constructed according to the structural characteristics of the macromolecule when simulating biological macromolecules, enabling the protein to always be in the set chemical environment;

[0011] S3. Use the steepest descent method to optimize the structure of the constructed bioprotein system, and obtain the stable configuration of the system through energy minimization calculations;

[0012] S4. Perform restricted molecular dynamics simulations on the stable configurations of the bioprotein system with minimized energy in the NVT and NPT ensembles respectively, so that the system relaxes to a stable state to prepare for the formal long-time molecular dynamics simulation;

[0013] S5. Perform a formal molecular dynamics simulation on the bioprotein system in a stable state to raise the temperature of the bioprotein system to the target temperature;

[0014] S6. Perform an unconstrained molecular dynamics simulation on the bioprotein system to obtain more comprehensive vibration information of the protein molecules at different temperatures;

[0015] S7. According to the vibration information of the protein molecules at different temperatures obtained in step S6, use the program documentation of Gromacs, and use the velacc post-processing program to process the conformational data for Fourier transform, calculate the spectrum of the autocorrelation of the current vector vq, convert it into infrared spectrum information, and obtain the infrared spectra at different temperatures as unbroadened spectrum data;

[0016] S8. Use the Matlab program to write a script to broaden the infrared spectrum obtained in step S7, perform curve fitting, and obtain the final infrared spectrum diagram.

[0017] Furthermore, in step S2, establish a periodic cubic box on the Charmm-gui website, place the complete structure of the biotoxin protein into the constructed box, add the SPC / E three-point water model in the box to provide the solvent environment, and thus construct a protein system model. According to the charged situation of the system model, add chloride ions and sodium ions to make the system solution condition neutral and maintain the same concentration as the physiological saline in the organism for the environment where the protein is located. Thus, a bioprotein system for infrared spectrum calculation is constructed.

[0018] Further, in step S3, the key parameters of the submitted script for structural optimization calculation are set by the following three means: generating a neighbor list using the Verlet method to calculate the movement trajectories of particles; processing the long-range electrostatic interactions by the PME method; and constraining the bonds connected to hydrogen by the LINCS method.

[0019] Further, in step S4, based on the image of the protein system energy changing with time, it is judged whether the entire protein structure has reached a stable structure. When the energy gradually tends to be stable and there is basically no change in the subsequent optimization process, it indicates that the biological protein system has been relaxed to the state with the lowest and stable energy.

[0020] Further, in step S5, during the long-time molecular dynamics simulation, the root mean square deviation and the radius of gyration of the protein are plotted according to the protein molecule coordinates, energy, and velocity information to judge whether the simulation process has reached equilibrium. After reaching equilibrium, it indicates that the system has risen to the target temperature.

[0021] Further, in step S6, a molecular dynamics simulation without setting constraint conditions is performed to obtain the complete coordinate, energy, and velocity information.

[0022] Further, in step S8, a script is written using the Matlab program to perform curve fitting on the data. During the fitting process, the Lorentz function is used for broadening to obtain the infrared spectra of the protein at different temperatures.

[0023] Further, in step S1, the initial structure of staphylococcal enterotoxin B is obtained from the proteindatabank website, and the hydrogen atoms are supplemented using the Discovery Studio software to obtain the complete structure of the staphylococcal enterotoxin B protein, which has a total of 3921 atoms; in step S2, a cube box with a size of 8.5 nm × 8.5 nm × 8.5 nm is established on the Charmm-gui website to simulate the boundary conditions and the internal environment. The complete structure of the staphylococcal enterotoxin B protein is placed in the box, and the distance between the box boundary and the protein is 0.8 nm. The SPC / E solvent model is used as the water model, and a total of 18790 water molecules are added; to maintain the electrical neutrality of the staphylococcal enterotoxin B protein system and to keep the concentration of the surrounding environment similar to that of the liquid salt solution in the organism, 56 chloride ions and 55 sodium ions are added to it.

[0024] Further, in step S3, the constructed biological protein system is minimized in energy using the steepest descent method, with the maximum number of steps set to 10,000 steps and the maximum step size to 0.001 nm; a neighbor list is generated using the verlet method to calculate the particle motion trajectory; long-range electrostatic interactions are treated by the PME method; bonds connected to hydrogen are constrained by the LINCS method; when the maximum force is less than 10.0 kJ / mol / nm, it is considered that the protein structure has converged, and according to the graph of the energy of the staphylococcal enterotoxin B protein system changing with time, it is judged whether the entire protein structure has reached a stable structure; in step S4, after finding the stable configuration of the system, an NVT simulation of 125 ps is first performed, and then an NPT simulation of 500 ps is performed; the pressure of the staphylococcal enterotoxin B protein system is stabilized at 1.05 bar using the Berendsen pressure bath, and the temperature is controlled at 254 K by the Velocity-rescale heat bath method.

[0025] Further, in step S6, an unconstrained simulation is performed on the staphylococcal enterotoxin B protein system, with a simulation duration of 400 ps, and data is recorded every 2 fs.

[0026] (III) Beneficial effects

[0027] The present invention proposes a variable-temperature infrared spectrum calculation method for protein biological toxins, which places the macromolecular biological protein structure into a box with boundaries established according to the actual chemical environment, and uses the molecular dynamics method to continuously cause random vibrations of the molecules in the box. During this period, the system records the change in the total dipole moment in the box at each moment, and then transforms the dipole moment signal into infrared frequency and intensity through Fourier transform, so as to obtain the infrared spectrum of the macromolecular biological protein, realizing the precise analysis and rapid detection of macromolecular biological toxins.

[0028] The present invention proposes a method for calculating the infrared spectrum of protein biological macromolecular toxins at different temperatures based on molecular dynamics, overcoming the difficulty that the first-principles and density functional theory methods cannot perform spectral calculations on large molecular weight systems due to their serious scaling problems. By establishing a large-scale system model of biological protein molecules in a similar actual chemical environment, a complete route of the infrared spectrum kinetic properties of this system is obtained using the molecular dynamics calculation method, so as to realize the infrared spectrum prediction of large biological protein molecules through computational means.

[0029] The beneficial effects of the present invention specifically include:

[0030] 1. The calculation method of the present invention can provide atomic-level resolution information about the protein structure and dynamic behavior, making the simulation of the infrared spectrum more accurate and reliable.

[0031] 2. By simulating the structure and dynamics of proteins at different temperatures through the present invention, the influence of temperature on infrared spectra can be studied, and the temperature dependence of protein structure can be revealed.

[0032] 3. Compared with laboratory experiments, the molecular dynamics simulation adopted by the present invention can provide more dynamic information, such as changes in vibration frequency and intensity, as well as conformational changes at different temperatures, so as to more comprehensively understand the structure and dynamics of proteins.

[0033] 4. Compared with laboratory experiments, the molecular dynamics simulation adopted by the present invention can save a large amount of cost and time. Especially when studying the infrared spectra of proteins at different temperatures, the preparation and repetition of experimental conditions can be avoided. Description of the Drawings

[0034] Figure 1 It is a schematic diagram of the overall structure of SEB protein;

[0035] In the figure: The red particles are water molecules, and the dark green ribbon-like substances with helical structures are the protein structures of SEB;

[0036] Figure 2 It is a schematic diagram of the protein structure of SEB protein;

[0037] In the figure: The light green particles are chloride ions, the blue particles are sodium ions, and the dark green ribbon-like substances with helical structures are the protein structures of SEB;

[0038] Figure 3 It is the change of SEB system energy over time;

[0039] In the figure: The vertical coordinate is the potential energy, with the unit of kilojoules per mole; the horizontal coordinate is the energy minimization time, with the unit of picoseconds;

[0040] Figure 4 It is a graph showing the change of the root mean square deviation (RMSD) of SEB over time at different temperatures;

[0041] In the figure: The vertical coordinate is the root mean square deviation, with the unit of nanometers; the horizontal coordinate is time, with the unit of picoseconds; among them, the black line represents the curve of the change of RMSD over time at 254K; the red line represents the curve of the change of RMSD over time at 264K; the blue line represents the curve of the change of RMSD over time at 274K; the green line represents the curve of the change of RMSD over time at 298K; the purple line represents the curve of the change of RMSD over time at 314K;

[0042] Figure 5Graph of the radius of gyration (Rg) of SEB as a function of time at different temperatures;

[0043] In the figure: The vertical axis is the mean square radius of gyration, with the unit of nanometer; the horizontal axis is time, with the unit of picosecond. Among them, the black line represents the curve of Rg changing with time at 254K; the red line represents the curve of Rg changing with time at 264K; the blue line represents the curve of Rg changing with time at 274K; the green line represents the curve of Rg changing with time at 298K; the purple line represents the curve of Rg changing with time at 314K;

[0044] Figure 6 Infrared spectra of SEB in the range of 2800 - 4000 cm -1 at different temperatures;

[0045] In the figure: The vertical axis is absorbance; the horizontal axis is wave number, with the unit of cm -1 ; among them, the black line represents the spectral absorption intensity at 254K; the red line represents the spectral absorption intensity at 264K; the blue line represents the spectral absorption intensity at 274K; the green line represents the spectral absorption intensity at 298K; the purple line represents the spectral absorption intensity at 314K;

[0046] Figure 7 Infrared spectra of SEB in the range of 3410 - 3500 cm -1 at different temperatures;

[0047] In the figure: Vertical axis: absorbance; Horizontal axis: wave number, with the unit of cm -1 ; among them, the black line represents the spectral absorption intensity at 254K; the red line represents the spectral absorption intensity at 264K; the blue line represents the spectral absorption intensity at 274K; the green line represents the spectral absorption intensity at 298K; the purple line represents the spectral absorption intensity at 314K;

[0048] Figure 8 Infrared spectra of SEB protein in the range of 500 - 2000 cm -1 at different temperatures;

[0049] In the figure: The vertical axis is absorbance; the horizontal axis is wave number, with the unit of cm -1 ; among them, the black line represents the spectral absorption intensity at 254K; the red line represents the spectral absorption intensity at 264K; the blue line represents the spectral absorption intensity at 274K; the green line represents the spectral absorption intensity at 298K; the purple line represents the spectral absorption intensity at 314K. Detailed implementation manners

[0050] To make the objectives, content and advantages of the present invention clearer, the following further describes in detail the specific implementation manners of the present invention with reference to the accompanying drawings and embodiments.

[0051] This embodiment provides a variable-temperature infrared spectroscopy calculation method for protein biotoxins. The method specifically includes the following steps:

[0052] S1. Establish the complete structure of the biotoxin protein

[0053] Obtain the initial structure of the biotoxin protein from the protein data bank website, and use DiscoveryStudio software to complete the hydrogen atoms to obtain the complete protein structure.

[0054] In this embodiment, the initial structure of staphylococcal enterotoxin B (SEB) (PDB: 3SEB) is obtained from the protein databank website, and the hydrogen atoms are completed using Discovery Studio software to obtain the complete SEB protein structure. The complete SEB protein structure has a total of 3921 atoms.

[0055] S2. Place the complete structure of the biotoxin protein into a box with periodic boundary conditions established according to the actual chemical environment. This box is a space for accommodating proteins and solvent environments constructed according to the structural characteristics of macromolecules when simulating biological macromolecules, so that the protein can always be in the set chemical environment.

[0056] Establish a periodic cubic box of appropriate size on the Charmm-gui website, place the complete structure of the biotoxin protein into the established box, add the SPC / E three-point water model in the box to provide a solvent environment, and thus construct a protein system model. According to the charged situation of the system model, add chloride ions and sodium ions to make the system solution condition neutral and maintain the environment where the protein is located similar to the concentration of physiological saline (0.15M) in the organism, thereby constructing a biological protein system for infrared spectroscopy calculation.

[0057] In this embodiment, a cubic box of 8.5nm × 8.5nm × 8.5nm is established on the Charmm-gui website to simulate the boundary conditions and internal environment. The complete SEB protein structure is placed in the box, and the box boundary is 0.8nm away from the protein. The explicit solvent model SPC / E is used as the water model, and a total of 18790 water molecules are added, as Figure 1As shown in the figure. To maintain the electrical neutrality of the SEB protein system and keep the environment similar to the concentration of the liquid salt solution in the organism, 56 chloride ions and 55 sodium ions are added to it. This biological protein system is used for subsequent calculations and simulations, such as Figure 2 as shown.

[0058] S3. The steepest descent method (steep) is used to optimize the structure of the constructed biological protein system, and the stable configuration of the system is obtained through energy minimization calculation.

[0059] The key parameters of the submitted script for the structure optimization calculation are set by the following three means:

[0060] (1) The Verlet method is used to generate the neighbor list and calculate the motion trajectory of the particles.

[0061] (2) The PME method is used to handle the long-range electrostatic interactions.

[0062] (3) The LINCS method is used to constrain the bonds connected to hydrogen.

[0063] In this embodiment, the steepest descent method (steep) is used to minimize the energy of the constructed biological protein system. The maximum number of steps is set to 10,000 steps, and the maximum step size is 0.001 nm. The Verlet method is used to generate the neighbor list and calculate the motion trajectory of the particles for efficient calculation of the interactions. The PME method is used to handle the long-range electrostatic interactions. To ensure that the structure and properties of the molecules during the simulation are maintained within an appropriate range, the LINCS method is used to constrain the bonds connected to hydrogen. When the maximum force is less than 10.0 kJ / mol / nm, it is considered that the protein structure has converged. According to the graph of the energy of the SEB protein system changing with time, it is judged whether the entire protein structure has reached a stable structure, such as Figure 3 as shown. It can be seen from the change trend of the energy that at the initial stage of the simulation, the energy of the SEB protein drops rapidly. Starting from 2000 steps, the energy gradually stabilizes, and there is basically no change in the subsequent optimization process, indicating that the SEB protein structure is relaxed to the state with the lowest energy and stability, and the stable configuration of the biological protein system is obtained.

[0064] S4. The restricted molecular dynamics simulations of the NVT and NPT ensembles are respectively carried out on the stable configuration of the biological protein system with minimized energy to relax the system to a stable state and prepare for the formal long-time molecular dynamics simulation.

[0065] After finding the stable configuration of the biological protein system, perform restricted molecular dynamics simulations in the NVT and NPT ensembles. The system after structural optimization is at the minimum energy, i.e., the most stable state, and this system is the protein and solvent environment system with periodic boundary conditions established above. Only after structural optimization can there be a stable configuration, and then performing restricted molecular dynamics simulations in the NVT and NPT ensembles can prevent the protein molecules from moving randomly before the water relaxes, so as to avoid unreasonable phenomena in the system.

[0066] According to the image of the protein system energy changing with time, judge whether the entire protein structure has reached a stable structure. When the energy gradually tends to be stable and basically does not change during subsequent optimization processes, it indicates that the biological protein system has been relaxed to the state with the lowest and stable energy.

[0067] In this embodiment, after finding the stable configuration of the system, first perform a 125 ps NVT simulation, and then perform a 500 ps NPT simulation. Use the Berendsen pressure bath to stabilize the pressure of the SEB protein system at 1.05 bar, and control the temperature at 254 K by the Velocity - rescale heat bath method.

[0068] S5. Perform a formal molecular dynamics simulation on the stable biological protein system to raise the temperature of the biological protein system to the target temperature.

[0069] During the long - time molecular dynamics simulation, there will be output files recording information such as the coordinates, energy, and velocity of protein molecules. By using this information to plot the Root Mean Square Deviation (RMSD) and the Radius of Gyration (Rg) of the protein, to judge whether the simulation process has reached equilibrium. After reaching equilibrium, it means that the system has been raised to the target temperature.

[0070] In this embodiment, after completing the restricted molecular dynamics simulation, start a 10 ns formal molecular dynamics simulation and record information such as coordinates, energy, and velocity. Judge whether the simulation process of the SEB protein structure has reached equilibrium through RMSD and Rg. As Figure 4 shown, the RMSD value of the SEB protein backbone increases significantly in the initial stage of the simulation and then tends to be stable. The RMSD values at 264 K are generally lower than those at other temperatures. There is a sudden drop in the RMSD value at 6900 ps, dropping to 0.112 nm, and then returning to the stable fluctuation range at 7600 ps. The molecular dynamics simulation process of the SEB protein at 274 K has the largest vibration amplitude. The RMSD value at 2000 ps is 0.140 nm, and it still shows fluctuations around 0.158 nm during the subsequent process of the protein system tending to be stable. As Figure 5As shown, at 264K, a relatively large Rg value indicates that the arrangement of SEB protein is relatively loose, while the Rg value at 274K is relatively small, indicating that the protein has a higher density at this temperature. In the later stage of the kinetic simulation, the Rg values at these temperatures are within a stable numerical range, approximately 1.865 ± 0.01 nm. The final result shows that the SEB protein at different temperatures has reached an ideal state and can be used for subsequent infrared spectroscopy simulation.

[0071] S6. Conduct an unconstrained molecular dynamics simulation on the biological protein system to obtain more comprehensive vibration information of the protein molecules at different temperatures.

[0072] Conduct a molecular dynamics simulation without setting constraints to ensure obtaining complete information such as coordinates, energy, and velocity, so as to obtain complete infrared spectroscopy information subsequently.

[0073] In this embodiment, an unconstrained simulation is conducted on the SEB protein system. The simulation duration is 400 ps, and data is recorded every 2 fs.

[0074] S7. According to the vibration information of the protein molecules at different temperatures obtained in step S6, using the program documentation of Gromacs, use the velacc post-processing program to process the conformational data for Fourier transform, calculate the spectrum of the autocorrelation of the current vector vq, and convert it into infrared spectroscopy information. At this time, the infrared spectra at different temperatures are unbroadened spectral data.

[0075] In this embodiment, using the program documentation of Gromacs, the velacc post-processing program calculates the spectrum of the autocorrelation of the current vector vq through 200,000 conformational data obtained to perform infrared spectroscopy simulation of SEB at different temperatures.

[0076] S8. Use a Matlab program to write a script to broaden the infrared spectrum obtained in step S7 and perform curve fitting to obtain the final infrared spectrum diagram.

[0077] Use a Matlab program to write a script to perform curve fitting on the data. During the fitting process, use the Lorentzian function to broaden and obtain the infrared spectra of the protein at different temperatures.

[0078] In this embodiment, the full width at half maximum (FWHM) value used is 20 cm -1 , and obtain the infrared spectra of the SEB protein at different temperatures. As Figure 6 shown, the infrared spectra of the SEB protein cover 2800 - 4000 cm -1Range, within which there are two strong vibration peaks, respectively caused by the stretching vibrations of O-H and N-H. Comparing the O-H stretching vibration peaks of SEB protein at different temperatures, the vibration frequency does not change significantly, and the shape of the vibration peak does not show obvious differences. However, there is a slight change in intensity. The intensity value of the O-H stretching vibration peak is the highest at 254K. In contrast, it is the lowest at 314K, with a 4% decrease in intensity. It is worth noting that the peak with a higher frequency shows a significant blue shift and peak broadening at 298K, as indicated by the arrow in the figure. Further observation Figure 7 of the infrared spectra at different temperatures (3410 - 3500 cm -1 ) shows that the peak in this range exhibits two bands (most significant at two temperatures of 264K and 274K). However, at 298K, this peak shows a single band. These changes in infrared spectral characteristics indicate that the structural dynamic changes of proteins at different temperatures can be revealed by analyzing the changes in characteristic peaks, thus providing a deeper understanding of the internal structure of SEB protein. As Figure 8 shown, the infrared spectra of SEB protein in the range of 500 - 2000 cm -1 at different temperatures are presented. This range has the vibration peaks of the amide I band and amide II band of the protein, which play an important role in studying the structural changes and kinetic processes of proteins. The amide I band is mainly the stretching vibration of the C=O bond on the protein backbone, and the peak position is usually at 1690 - 1630 cm -1 , and this vibration mode has strong sensitivity and specificity to the secondary structure of proteins. The amide II band is the stretching and bending vibration of N-H, which can be used to understand the changes and responses of proteins under different conditions. The frequency of the amide I band of SEB protein is the same at different temperatures, and the vibration peak is located at 1650 cm -1 . However, the presence of isomers can be clearly observed at 254K and 314K, which affects the shape of this vibration peak ( Figure 8 the red oval part in). At 264K, there are two relatively weak peaks beside the amide I vibration peak, and the peak positions are 1606 cm -1 and 1630 cm -1 . Perhaps at this temperature, the kinetic energy of the molecules is low, and the protein molecules are more likely to stay in the energy states corresponding to various vibration modes, so there are two obvious peaks. At other temperatures, the peak at 1630 cm -1 is very weak, especially at 298K, one peak disappears, and the peak position is 1613 cm -1 . And at 1520 cm -1The above situation also exists in the amide II band vibration of the protein, and the vibration peaks at 298K and 314K both disappear. This change indicates that temperature affects the average kinetic energy of the protein, thereby causing changes in the occupancy of the molecular vibration modes.

[0079] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the technical principle of the present invention, several improvements and modifications can be made, and these improvements and modifications should also be regarded as the protection scope of the present invention.

Claims

1. A method for calculating the temperature-variable infrared spectrum of proteinaceous biological toxins, characterized in that: The variable temperature infrared spectrum calculation method comprises the following steps: S1. Establish the complete structure of the biotoxin protein The initial structure of the biotoxin protein was obtained from the protein data bank website, and the hydrogen atoms were filled in using Discovery Studio software to obtain the complete protein structure; S2. Place the complete structure of the biotoxin protein into a box with periodic boundary conditions established according to the actual chemical environment. The box is a space that contains the protein and the solvent environment and is constructed according to the structural characteristics of the macromolecule when simulating biological macromolecules, so that the protein can always be in the set chemical environment; S3. Use the steepest descent method to optimize the structure of the constructed biological protein system and obtain the stable configuration of the system through energy minimization calculation; S4. Perform restricted molecular dynamics simulations of NVT and NPT ensembles on the stable configurations of the energy-minimized biological protein system, so that the system can relax to a stable state and prepare for formal long-term molecular dynamics simulations; S5. Perform a formal molecular dynamics simulation on the biological protein system in a stable state to raise the temperature of the biological protein system to a target temperature; S6. Perform an unconstrained molecular dynamics simulation on the biological protein system to obtain more comprehensive vibration information of protein molecules at different temperatures; S7. Based on the vibration information of the protein molecule at different temperatures obtained in step S6, the conformational data is processed using the velacc post-processing program using the Gromacs program document to perform Fourier transform, the spectrum of the current vector vq autocorrelation is calculated, and the spectrum is converted into infrared spectrum information, and the infrared spectra at different temperatures are obtained as unbroadened spectrum data; S8. Use Matlab program to write a script to broaden the infrared spectrum obtained in step S7, perform curve fitting, and obtain the final infrared spectrum.

2. The variable temperature infrared spectrum calculation method of protein biotoxins according to claim 1, characterized in that: In step S2, a periodic cube box is established on the Charmm-gui website, the complete structure of the biological toxin protein is placed in the established box, and the SPC / E three-point water model is added to the box to provide a solvent environment to construct a protein system model. According to the charge of the system model, chloride ions and sodium ions are added to make the system solution condition neutral and maintain the protein environment with the same concentration as the physiological saline in the organism, thereby constructing a biological protein system for infrared spectrum calculation.

3. The variable temperature infrared spectrum calculation method of protein biotoxins according to claim 1, characterized in that: In step S3, the key parameters of the script submitted for the structure optimization calculation are set by the following three means: using the verlet method to generate a neighbor list and calculate the motion trajectory of the particles; using the PME method to deal with long-range electrostatic interactions; and using the LINCS method to constrain the bonds connected to hydrogen.

4. The variable temperature infrared spectrum calculation method of protein biotoxins according to claim 1, characterized in that: In step S4, based on the image of the change of protein system energy over time, it is judged whether the entire protein structure has reached a stable structure. When the energy gradually tends to be stable and the energy basically does not change during the subsequent optimization process, it indicates that the biological protein system is relaxed to the lowest energy and stable state.

5. The variable temperature infrared spectrum calculation method of protein biotoxins according to claim 1, characterized in that: In step S5, during the long-term molecular dynamics simulation, the root mean square deviation and protein gyration radius are plotted according to the protein molecular coordinates, energy and velocity information to determine whether the simulation process has reached equilibrium. If equilibrium is reached, it means that the system has risen to the target temperature.

6. The variable temperature infrared spectrum calculation method of protein biotoxins according to claim 1, characterized in that: In step S6, a molecular dynamics simulation is performed without setting constraints to obtain complete coordinate, energy and velocity information.

7. The variable temperature infrared spectrum calculation method of protein biotoxins according to claim 1, characterized in that: In step S8, a script is written using the Matlab program to perform curve fitting on the data. During the fitting process, the Lorentz function is used for broadening to obtain infrared spectra of the protein at different temperatures.

8. The variable temperature infrared spectrum calculation method of protein biotoxins according to claim 1, characterized in that: In step S1, the initial structure of Staphylococcal enterotoxin B was obtained from the protein databank website, and the hydrogen atoms were supplemented using Discovery Studio software to obtain the complete structure of Staphylococcal enterotoxin B protein. The complete structure of Staphylococcal enterotoxin B protein has a total of 3921 atoms; in step S2, a 8.5nm×8.5nm×8.5nm cubic box was established on the Charmm-gui website to simulate boundary conditions and internal environment, and the complete structure of Staphylococcal enterotoxin B protein was placed in the box. The box boundary was 0.8nm away from the protein, and the display solvent model SPC / E was used as a water model, and a total of 18790 water molecules were added; in order to maintain the electrical neutrality of the Staphylococcal enterotoxin B protein system and maintain the environment similar to the concentration of liquid salt solution in the organism, 56 chloride ions and 55 sodium ions were added thereto.

9. The variable temperature infrared spectrum calculation method of proteinaceous biotoxins according to claim 8, characterized in that: In step S3, the steepest descent method is used to minimize the energy of the constructed biological protein system, and the maximum number of steps is set to 10,000 steps and the maximum step length is 0.001 nm; the verlet method is used to generate a neighbor list and calculate the motion trajectory of the particles; the long-range electrostatic interaction is handled by the PME method; the bonds connected to hydrogen are constrained by the LINCS method; when the maximum force is less than 10.0 kJ / mol / nm, the protein structure is considered to have converged, and the energy change of the Staphylococcus enterotoxin B protein system over time is used to determine whether the entire protein structure has reached a stable structure; in step S4, after finding a stable configuration of the system, a 125 ps NVT simulation is performed first, and then a 500 ps NPT simulation is performed; the pressure of the Staphylococcus enterotoxin B protein system is stabilized at 1.05 bar using a Berendsen pressure bath, and the temperature is controlled at 254 K using the Velocity-rescale heat bath method.

10. The variable temperature infrared spectrum calculation method of protein biotoxins according to claim 9, characterized in that: In step S6, a simulation without setting constraints is performed on the Staphylococcus enterotoxin B protein system, the simulation duration is 400 ps, ​​and data is recorded every 2 fs.

Images