CN115662496A - Motion mode simulation method of protein molecule machine - Google Patents

Abstract

The invention provides a motion mode simulation method of a protein molecule machine, which comprises the following steps: s10, acquiring a protein molecule machine structure data file; s20, simulating the molecular dynamics of a protein molecular machine to obtain a topological file and a coordinate file of the movement of the compound, a receptor thereof and a ligand thereof, and calculating the binding free energy of the compound; s30, predicting the motion track of the protein molecular machine based on the built NRI model and the coordinate file, wherein the NRI model comprises an encoder for predicting the given dynamic system track interaction and a decoder for predicting the given dynamic system track in the interaction diagram; s40, observing and recording the three-dimensional structure and the motion trail of the protein molecule machine based on VR equipment. The method combines molecular dynamics simulation and a deep neuron network to realize accurate simulation and prediction of the motion mode of a protein molecular machine.

Description

Motion mode simulation method of protein molecule machine

Technical Field

The invention relates to the technical field of vehicle control, in particular to an online estimation method and device for a tire-road adhesion coefficient.

Background

The protein molecular machine is a machine which is composed of substances with molecular scale and can perform certain processing function, and the components of the machine are mainly biomolecules such as protein and the like, and participate in a plurality of important life processes. The research on the structure, the function and the motion mode of a protein machine is the basis for understanding the high-dimensional complex molecular machine and has important significance for revealing the essence of life phenomena. The latest severe acute respiratory syndrome coronavirus variant 2 (SARS-CoV-2) Omicron caused panic worldwide due to its infectivity and vaccine escape mutations. The infectivity and antibody resistance of the SARS-CoV-2 variant depends on its mutation in the spike S protein Receptor Binding Domain (RBD). However, a large number of mutations have been made in the RBD of spike protein by Omicron, which has raised a high level of interest to the scientific community and the public. It is also important to explore the binding between the protein of the Omicron variant S and the human receptor ACE 2.

Disclosure of Invention

Aiming at the problems, the invention provides a motion mode simulation method of a protein molecular machine, which combines molecular dynamics simulation and a deep neuron network to realize accurate simulation and prediction of the motion mode of the protein molecular machine.

The technical scheme provided by the invention is as follows:

a method of motion pattern simulation of a protein molecular machine, comprising:

s10, acquiring a protein molecular machine structure data file, wherein the protein molecular machine comprises a complex of a ligand wild type and an Omicron Ba.2 RBD region S protein and a receptor ACE 2;

s20, simulating the molecular dynamics of the protein molecular machine to obtain a topological file and a coordinate file of the movement of the compound and the receptor and the ligand thereof, and calculating the binding free energy of the compound;

s30, predicting the motion track of the protein molecule machine based on the built NRI model and the coordinate file, wherein the NRI model comprises an encoder used for predicting the given dynamic system track interaction and a decoder used for predicting the given dynamic system track in the interaction diagram; in the NRI model, the motion coordinate of the compound is input, and the predicted motion coordinate of the compound is output;

s40, observing and recording the three-dimensional structure and the motion trail of the protein molecule machine based on a VR device.

The invention provides a motion mode simulation method of a protein molecular machine, which explores the dynamic interaction between the RBD region S protein of a ligand Omicron Ba.2 and a receptor ACE2 by adopting a molecular dynamics simulation method, and compares the dynamic interaction with a wild type and ACE2 system. The simulation finds that Omicron shows stronger binding capacity with human cells, and many important mutations occur in the RBD contacted with ACE2, which indicates that the mutation on the RBD may cause the infectivity to be enhanced, and provides a partial basis for subsequent medical diagnosis and treatment.

Drawings

The foregoing features, technical features, advantages and embodiments are further described in the following detailed description of the preferred embodiments, which is to be read in connection with the accompanying drawings.

FIG. 1 is a schematic flow chart of a method for simulating a motion pattern of a protein molecular machine according to the present invention.

Detailed Description

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following description will be made with reference to the accompanying drawings. It is obvious that the drawings in the following description are only some examples of the invention, and that for a person skilled in the art, other drawings and embodiments can be derived from them without inventive effort.

Molecular dynamics is a computer simulation experiment method based on molecular mechanics, and is mainly used for simulating and observing motion changes of molecules in different environments. The molecular dynamics simulation is a simulation method which is closest to experimental conditions in the molecular simulation, can give out a microscopic evolution process of a system from an atomic level, intuitively shows the mechanism and rule of the occurrence of experimental phenomena, and promotes the development of research towards a direction with higher efficiency, more economy and more predictability, so the molecular dynamics simulation plays an increasingly important role in the research of biology, pharmacy, chemistry and material science.

In recent years, deep learning has been remarkably advanced, and remarkable results are obtained from models, algorithms to large-scale application. The appearance of deep learning is an important revolution of machine learning and is a great thrust for the development of artificial intelligence. The traditional neural network is a shallow machine learning, and the deep learning is a new generation neural network under the development of the traditional neural network. Deep learning is to realize low-level to high-level feature extraction of externally input data by creating and simulating an information processing neural structure of a human brain, so that a machine can understand learning data to obtain information.

VR relies on computer platform and advanced electronic technology, through creating a simulation environment space of high fidelity, through using multiclass sensing equipment, further improves user's use experience, makes experience user carry out actual operation with more direct-viewing and tangible state to the object in the virtual space, has realized the optimization of simulation experience. The experiment utilizes VR equipment to look over the three-dimensional structure of the protein machine more intuitively, and the motion mode of the protein machine is better understood.

Based on this, the present invention provides a method for simulating a movement pattern of a protein molecular machine, as shown in fig. 1, comprising:

s10, acquiring a protein molecular machine structure data file, wherein the protein molecular machine comprises a complex of a ligand wild type and an Omicron Ba.2 RBD region S protein and a receptor ACE 2;

s20, simulating the molecular dynamics of a protein molecular machine to obtain a topological file and a coordinate file of the movement of the compound and a receptor and a ligand thereof, and calculating the binding free energy of the compound;

s30, predicting the motion track of the protein molecular machine based on the built NRI model and the coordinate file, wherein the NRI model comprises an encoder for predicting the given dynamic system track interaction and a decoder for predicting the given dynamic system track in the interaction diagram; in the NRI model, the motion coordinate of the compound is input, and the predicted motion coordinate of the compound is output;

s40, observing and recording the three-dimensional structure and the motion trail of the protein molecule machine based on the VR equipment.

Specifically, in step S10, wild-type (before mutation) and Omicron Ba.2 (Omicron variant) complexes of the RBD region S Protein and ACE2, which are all from the Protein Data Bank (PDB) website, were prepared, wherein the complex with Omicron Ba.2 was taken from pdbID 7ZF7 and the complex with wild-type was taken from pdbID 6LZG. In addition, the PDB file of the complex needs to be resolved into receptor and ligand using Pymol software to facilitate subsequent binding free energy calculations (in both complexes, wild-type and Omicron ba.2 RBD region S protein is the ligand and ACE2 is the receptor). The PDB file is also preprocessed by using UltraEdit before molecular dynamics simulation is carried out on the compound, the receptor and the ligand, only Atom field information is reserved, and the disulfide bond of the Atom field information is processed.

In step S20, amber is used to simulate the molecular dynamics of the protein molecular machine to obtain a trajectory file of the complex. Specifically, wild-type and composite systems of Omicron ba.2 and ACE2 were simulated using software Amber16 and leaprc.ff03ua force fields, respectively: firstly, molecular construction is carried out under a linux system, and a topological file and a coordinate file (used as references for subsequent predicted tracks) of the compound, the receptor and the ligand are obtained. In this procedure, the complex is centered in a cubic box and appropriate ions are added to the system to electrostatically neutralize the molecular system and dissolve with a tip3pbox water model to obtain the topology and coordinate file of the complex in solvent. Thereafter, the solvent complex is energy minimized and after a 50ps temperature rise, the complex is subjected to a 50ps density balance under weakly limiting conditions and an isostatic balance of 500ps at a temperature of 300K. In the whole simulation process, SHAKE constraint is applied to hydrogen atoms during operation, and the time step length of 2fs and Langevin dynamics are used for controlling the temperature. The Root Mean Square Deviation (RMSD) was calculated relative to the protein backbone of the energy minimized structure to determine if the conformation had stabilized during equilibrium. Here, balance is madeWhether the conformation during the period is stable is determined by whether the calculated RMSD fluctuates within a certain interval, e.g., a certain time period after the value starts from 0

To

Indicating conformational stability during equilibrium. If a sudden increase or decrease occurs, it is not stable.

After assessing whether the molecular dynamics simulation balances the solvent complex state based on RMSD properties, the energy term for S protein-ACE 2 interaction was further calculated using a script, which is based on a molecular mechanics with poisson-boltzmann surface area (MM-PBSA) method.

Ideally, the binding free energy is calculated directly using receptor and ligand binding in solvent, but in simulations of these solvation states, most of the energy contribution will come from solvent-solvent interactions and the total energy fluctuations will be an order of magnitude greater than the binding energy, so the calculation of the binding energy will take a very long time to converge. It is therefore more efficient to calculate the parts of the thermodynamic cycle separately, simplifying the procedure by first combining the receptor and ligand in the gas phase and then placing the combination in a solvent environment. Since the path states are the same between the beginning and the end, the binding free energy Δ G is considered from the viewpoint of energy conservation _Solv,bind Is represented by the formula (1):

ΔG _Solv,bind ＝ΔG _Gas,bind +ΔG _Solv,comp -(ΔG _Solv,Lig +ΔG _Solv,Rec )(1)

wherein, Δ G _Gas,bind Indicates the energy of ligand-receptor binding in solvent, Δ G _Solv,comp Representing the energy converted by the complex in the gaseous and solvent states, AG _Solv,Lig Representing the energy of conversion of the ligand in the gaseous and solvent states, AG _Solv,Rec Represents the energy of the acceptor converted in the gaseous and solvent states;

bonding free from the viewpoint of bonding energyEnergy Δ G _bind Further can be represented by formula (2): :

ΔG _bind ＝ΔH-TΔS(2)

wherein Δ H represents enthalpy change and can be decomposed into gas phase energy Δ E _MM And solvation energy AG _sol Neglecting entropy (-T Δ S) versus energy Δ G _bind The enthalpy change Δ H is expressed by the formulas (3) to (5):

ΔH＝ΔE _MM +ΔG _sol (3)

ΔE _MM ＝ΔE _ele +ΔE _vdW +ΔE _int (4)

ΔG _sol ＝ΔG _pb +ΔG _np (5)

wherein, delta E _ele Denotes the electrostatic term,. DELTA.E _vdW Denotes the Van der Waals term,. DELTA.E _int Represents the internal energy term (Δ E) _int Internal energy including bond energy, bond angle energy, and dihedral angle energy), Δ G _pb Denotes polar solvation energy,. DELTA.G _np Representing a non-polar solvation energy.

In addition to calculating the binding free energy, after obtaining a trajectory file based on the above dynamic simulation process, performing calculation analysis on Root Mean Square Fluctuation (RMSF), solvent reachable surface area (sasa) and minimum residue distance (minist) by amber software, wherein RMSF shows the comparison of the moving distance between a residue and the original position, and the smaller the movement is, the better the binding is shown, the stronger the binding free energy is, and the higher the solvent reachable surface area sasa is; the smaller the minimum residue distance minist, the stronger the binding free energy, the distance being less than

The residue of (b) is considered a contact residue. In addition, the VMD tool is used to analyze hydrogen bonding conditions such that the distance between the donor atom and the acceptor atom is less than

And the angle formed by the hydrogen atoms connected with the donor-acceptor-donor is less than 45 degrees as a standard, and ligaplot software is used for analyzing the 2D structure of the protein-ligand so as to help explain the conclusion that the combination free energy is strengthened after mutation.

After the calculation of the binding free energy is completed, in step S30, the motion trajectory of the protein molecular machine is predicted based on the built NRI model and the coordinate file, where the NRI model includes an encoder for predicting the interaction of the given dynamic system trajectory and a decoder for predicting the given dynamic system trajectory in the interaction graph; in the NRI model, the input is the motion coordinates of the compound, and the output is the predicted motion coordinates of the compound. Specifically, the input to the NRI model consists of N nodes, with the feature vector (position and velocity in the X, y and z dimensions) of node i represented as X at time t _i ^t The feature set of all N nodes is denoted as X ^t ＝{X ₁ ^t ,...,X ^t _N }. The locus of node i is denoted X _i ＝{X _i ¹ ,...,X _i ^t Where t represents the number of time steps. Finally, the process is carried out in a closed loop, all track data recorded as X = { X ¹ ,...,X ^t }. The NRI model learns edge values simultaneously and reconstructs future trajectories of the dynamic system in an unsupervised manner based on the unknown map z. Interaction between nodes i and j with latent variable Zj ⁱ Is of the form {1, ·, K }, where K is the number of interaction types modeled. These interaction types do not have any predefined meaning, but the model learns to assign a meaning to each type.

After the predicted track file is obtained in step S30, the predicted track file is placed in Assets file folder, and resources are placed in these folders in a classified manner. Here, a model file in fbx format is used, and attributes such as animation, material, and cells can be derived and used simultaneously with the model. And optimizing the three-dimensional scene and the role model, and realizing real-time rendering of the three-dimensional grid in the interactive system. Specifically, an Xfrog Organic modeling tool is used for making the whole process of combining wild type and Omicron S protein with human ACE2, and the efficiency of making three-dimensional animation is improved.

The VR device includes: the physical space positioning module is used for positioning the relative position of the user in the real space through the head-mounted sensing equipment to acquire the spatial position data of the user; the VR environment interaction module is used for generating interaction information based on the space position acquired by the physical space positioning module and the interaction action captured by the handheld sensor; the message processing module is used for processing the interactive information generated by the VR environment interactive module; the VR virtual environment simulation module is used for updating the parameters of the scene object according to the interactive information generated by the message processing module and finishing the rendering output of a new frame; and feeding back the virtual environment interaction information to the user sensor and the physical space positioning module.

In work, the physical space positioning module positions the relative position of the user in the real space through the head-mounted sensing equipment, and spatial position data of the user are obtained. And after the spatial position and the interactive behavior data are acquired, the interactive behavior data are handed to an interactive system for interactive behavior interpretation. The position information in the real space can be mapped to different position nodes in the corresponding virtual environment through coordinate change and conversion. Interactive action data captured by the handheld sensor are uniformly converted into different types of message structural bodies and sent to the message processing modules, corresponding interactive response data and processing events are transmitted to the VR virtual environment simulation system through the interactive behavior feedback module to update parameters of scene objects, rendering output of a new frame is completed, and virtual environment interactive information is fed back to the user sensor and the spatial positioning device, so that the motion trajectory of the protein subjected to molecular dynamics simulation can be observed more clearly.

According to the invention, a molecular dynamics simulation method is adopted to explore the dynamic interaction between ACE2 and RBD, and the comparison with a wild type ACE2 system and an ACE2 system is carried out, and the result shows that Omicron shows stronger binding capacity with human cells, and a plurality of important mutations occur in RBD contacted with ACE2, which indicates that the mutation on the RBD may cause the infectivity enhancement. Meanwhile, the binding property of a receptor ACE2 and two RBD region S proteins of a wild type and an Omicron Ba.2 is researched by a molecular dynamics method, the aspects of hydrogen bonds, the number of contact residues, the accessible surface area of a solvent and the like are analyzed, and the binding free energy of the RBD region S proteins and the ACE2 is calculated.

In the process of predicting the motion trail, an NRI model is set up, and an optimal prediction model is selected according to the accuracy of a training set and a testing set. Specifically, anaconda software is applied, and a deep learning system is built by adopting a python language. python is an object-oriented interpreted computer programming language that uses features that are cross-platform. The BP network nerve has the biggest characteristic that signals are propagated in the forward direction, errors are propagated in the backward direction, and the method can be understood as adding input data into the signals which are propagated in the forward direction continuously, and observing and feeding back the obtained result, so that a model is optimized better, and the method is understood as the backward propagation of the errors. The learning rule is that the steepest descent method is used, and the weight and the threshold value of the network are continuously adjusted through back propagation, so that the error square sum of the network is minimum.

Protein allosterism is a biological process facilitated by spatially long distance intra-protein communication in which ligand binding or remote amino acid changes affect the active site remotely. Molecular Dynamics (MD) simulations provide a powerful computational method for exploring allosteric effects. However, current MD simulations do not achieve the time scale of the entire allosteric process. The advent of deep learning makes it possible to evaluate short-range and long-range communications over space to understand the allosterism. To this end, a neural relationship inference model based on graph neural networks is applied, which employs an encoder-decoder architecture while inferring potential interactions to probe protein allosteric processes as a dynamic network of interacting residues. From the MD trajectories, the model successfully learns long-range interactions and pathways, and furthermore, the model can discover allosterically related interactions earlier in the MD simulated trajectories and predict relative free energy changes at sudden changes more accurately than other methods.

Molecular dynamics simulation can directly detect biomolecular motion, but due to the limited time scale of the simulation and the high dimensionality and complexity of the 3D trajectory data, meaningful functional information may not be captured. Furthermore, many challenging MD analysis problems lack a suitable method to detect remote communications. Computational techniques for modeling protein allosteric communication rely on graph-theoretic metrics to identify long-range couplings between two remote active sites. In general, a protein can be mapped to a graph, where each node represents a residue and each weighted edge represents an interaction between two nodes. The shortest path between the allosteric site and the active site in a protein may be important for propagating signals in allosteric communication. Early graphical models used static crystal structures to calculate the shortest path between one residue and the other, which may not account for all potential contacts in dynamic proteins and the associated allosteric behavior. The latter well-known allosteric method, the disturbance response sweep (PRS), uses a Hessian-based Elastic Network Model (ENM) to obtain the relative dynamics of a location. The model studies how perturbations on a single residue trigger perturbation (signal) cascades to other nodes in the elastic network, thus achieving allosteric communication. To more accurately model the response upon ligand binding or mutation, the inverse of the Hessian model was replaced with a covariance matrix containing the system dynamics.

The NRI model is suitable for learning simulated motion trajectories of biomacromolecules in MD simulation, wherein the biomolecular molecules are formed by chemical bond-connected atoms, and the motion rule of the biomolecular molecules is described by Newton mechanics. The model learns the embedding of network dynamics by minimizing reconstruction errors between reconstructed trajectories and simulated trajectories using GNNs; the NRI model then infers the edges between the residuals represented by the latent variables. The chemical insertions essentially abstract the fundamental role of key residues in conformational transitions, which helps to decipher the mechanisms of protein allosterism.

Unity3D is a specialized cross-platform game development and virtual reality engine developed by Unity Technologies, supporting multiple scripting languages such as C #, javascript, and the like. The invention uses VR equipment HTC Vive which is well received by the industry, combines a Unity engine as a development platform, and is researched and applied to a biological control VR simulation platform. Finally VR is used to observe the cell structure and its motion trajectory.

It should be noted that the above embodiments can be freely combined as necessary. The foregoing is only a preferred embodiment of the present invention, and it should be noted that, for persons skilled in the art, numerous modifications and adaptations can be made without departing from the principle of the present invention, and such modifications and adaptations should be considered as within the scope of the present invention.

Claims (6)

1. A method of simulating a movement pattern of a protein molecular machine, comprising:

s10, acquiring a structure data file of a protein molecular machine, wherein the protein molecular machine comprises a complex of ligand wild type and Omicron Ba.2 RBD region S protein and a receptor ACE 2;

s20, simulating the molecular dynamics of the protein molecular machine to obtain a topological file and a coordinate file of the movement of the compound and the receptor and the ligand thereof, and calculating the binding free energy of the compound;

s30, predicting the motion track of the protein molecular machine based on the built NRI model and the coordinate file, wherein the NRI model comprises an encoder for predicting the given dynamic system track interaction and a decoder for predicting the given dynamic system track in the interaction diagram; in the NRI model, the motion coordinate of the compound is input, and the predicted motion coordinate of the compound is output;

s40, observing and recording the three-dimensional structure and the motion trail of the protein molecule machine based on a VR device.

2. The method for simulating a motion pattern of a protein molecular machine according to claim 1, wherein in step S20, amber is used to simulate the molecular dynamics of the protein molecular machine to obtain a trace file of the complex.

3. The method for modeling a movement pattern of a protein molecular machine according to claim 1, wherein in step S20, the binding free energy Δ G _Solv,bind Can be expressed as:

ΔG _Solv,bind ＝ΔG _Gas,bind +ΔG _Solv,comp -(ΔG _Solv,Lig +ΔG _Solv,Rec )

wherein, Δ G _Gas,bind Indicates the energy of ligand-receptor binding in solvent, Δ G _Solv,comp Representing the energy converted by the complex in the gaseous and solvent states, AG _Solv,Lig Representing the energy of conversion of the ligand in the gaseous and solvent states, AG _Solv,Rec Representing the energy of the acceptor conversion in the gaseous and solvent states.

4. The method for modeling a movement pattern of a protein molecular machine according to claim 1, wherein in step S20, the binding free energy Δ G _bind Can also be expressed as:

ΔG _bind ＝ΔH-TΔS

wherein Δ H represents enthalpy change and can be decomposed into gas phase energy Δ E _MM And solvation energy AG _sol Neglecting the entropy T Δ S vs. energy Δ G _bind The enthalpy change Δ H is expressed as:

ΔH＝ΔE _MM +ΔG _sol

ΔE _MM ＝ΔE _ele +ΔE _vdW +ΔE _int

ΔG _sol ＝ΔG _pb +ΔG _np

wherein, delta E _ele Denotes the electrostatic term,. DELTA.E _vdW Denotes the Van der Waals term,. DELTA.E _int Represents an internal energy term, Δ G _pb Denotes polar solvation energy,. DELTA.G _np Representing a non-polar solvation energy.

5. The method for simulating the movement pattern of a protein molecular machine according to claim 1, 2, 3 or 4, wherein in step S40, VR device generates the corresponding fbx file based on pdb trace file generated by NRI model, and obtains and displays the whole process of the wild type and Omicron Ba.2 RBD region S protein binding with the receptor ACE 2.

6. The method of claim 5, wherein the VR device comprises:

the physical space positioning module is used for positioning the relative position of the user in the real space through the head-mounted sensing equipment to acquire the spatial position data of the user;

the VR environment interaction module is used for generating interaction information based on the space position acquired by the physical space positioning module and the interaction action captured by the handheld sensor;

the message processing module is used for processing the interactive information generated by the VR environment interactive module;

the VR virtual environment simulation module is used for updating the parameters of the scene object according to the interactive information generated by the message processing module and finishing the rendering output of a new frame; and feeding back the virtual environment interaction information to the user sensor and the physical space positioning module.

Images