CN103699816A - Monte Carlo simulation-based protein thermomechanical-analysis method - Google Patents

Abstract

The invention relates to the technical fields of biological computing and Monte Carlo simulation, and particularly discloses a Monte Carlo simulation-based protein thermomechanical-analysis method, which comprises the following steps of A, determining a protein energy model; B, simulating and calculating the state density of a protein system. According to the method, the whole thermomechanical process of protein folding can be efficiently analyzed and researched to further explore and research a protein folding process.

Description

Protein thermodynamics analytical approach based on Monte Carlo simulation

Technical field

The present invention relates to biological computation and Monte Carlo simulation field, specifically disclose a kind of protein thermodynamics analytical approach based on Monte Carlo simulation.

Background technology

Protein engineering is the Disciplinary Frontiers of modern biotechnology development, and its basic goal is naturally occurring protein will be transformed according to people's imagination, or designs as required non-natural novel protein with certain specific function.And one of important foundation of this transformation and design is the prediction of protein folding structure.The form of protein folding structure has determined its biological function that may have to a great extent.In other words, between the structure of protein and function, there is consistance.Therefore to the prediction of protein folding structure and research in protein engineering, field of medicaments etc. have extremely important meaning.

Monte Carlo simulation is a kind of random search algorithm being widely used, by the particle position removal search conformational space of randomly changing protein system, the density of states of acquisition system, thereby various macroscopic properties that can computing system, and then the folding whole thermodynamic process of Study on Protein.But because protein energy space is complicated, between the local lowest region of energy value and have energy barrier between overall lowest region, the simulation at low temperatures of traditional Monte carlo algorithm is easily absorbed between the local lowest region of energy, is difficult to jump out and finds globally optimal solution.A lot of new Monte Carlo simulation algorithms have been produced thus, as many canonical ensembles (multicanonical) Monte carlo algorithm, simulated tempering (simulated tempering) algorithm, copy exchange (Replica Exchange Method) algorithm and Wang-Landau algorithm.

The basic thought of Wang-Landau algorithm is to use non-Boltzmann distribution function to realize the analog simulation of freely walking between energy range, thereby searches for more widely configuration space; Also can automatically obtain the density of states of protein system simultaneously, thereby can calculate the many cannoncial system thermodynamic quantities in broad temperature range, and then the folding whole thermodynamic process of Study on Protein.Wang-Landau is due to directly perceived succinct, and its range of application has expanded in magnetic system, liquid crystal, liquid, cluster, Spin glass model and protein folding.But Wang-Landau algorithm need further improvement in computational accuracy and speed.

Summary of the invention

The present invention is intended to overcome the defect of prior art computational accuracy and underspeed in protein thermodynamics is analyzed, and a kind of protein thermodynamics analytical approach based on Monte Carlo simulation is provided.

Technical scheme of the present invention is 1, a kind of protein thermodynamics analytical approach based on Monte Carlo simulation, it is characterized in that, comprises steps A: determine protein energy model; And B: simulation and the calculating protein system density of states.

Further, steps A further comprises: adopt ECEPP energy force field model and angle coordinate system, the expression-form in the ECEPP energy field of force is:

E _ECEPP=E _C+E _LJ+E _HB+E _Tor

Wherein

the coulomb acting force between two electric charges, r _ijrepresent the distance between atom i and j,

the Lan Na-Jones acting force between two atoms,

hyarogen-bonding, E _tor=∑ _lu _l(1 ± cos (n _lξ _l)) be dihedral turning effort power, ξ _ll dihedral.

Further, steps A further comprises: to carrying out discretize processing between used protein energy range, if get k energy bin interval value, need [E _min, E _max] on average divide k bin interval, with each interval middle energy value, represent this energy interval value.

Further, step B further comprises: by the MPI concurrent program algorithm of principal and subordinate's process mode, and the density of states of simulation and calculating protein system.

Further, in N minute process of described principal and subordinate's process mode, a minute process 1 is host process, and within all the other minute, process is subprocess.

Further, described host process comprises the following steps:

S11: logarithm S (the E)=lng (E)=0 of initialization protein system Density function, histogram H (E)=0(E _min≤ E≤E _max), s=1, modifying factor df| _e=(κ Θ (E ₀-E)+1) lnf, wherein Θ (E ₀-E) be Heaviside piecewise function, κ, E ₀, f is the parameter that is relevant to model;

S12：t=1；

S13: in host process, original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration, t=t+1, s=s+1, modifying factor is df| _e=(κ Θ (E ₀-E)+1) lnf;

Step S13 circulation tmax time;

S14: all process intercommunications, host process is collected all S from process _tmpand H (E) _tmp(E) and accumulation calculating go out overall S (E) and H (E), the overall situation S (E)=S (E)+all S from process _tmp(E), H (E)=H (E)+all H from process of the overall situation _tmp(E), being broadcast to of overall S (E) and H (E) is all from process, try to achieve H _real(E)=H (E) * [1+ κ Θ (E ₀-E)], the mild condition of judgement histogram:

$\frac{\max \left(H_{\mathrm{real}}\left(E\right)\right)-\min \left(H_{\mathrm{real}}\left(E\right)\right)}{\max \left(H_{\mathrm{real}}\left(E\right)\right)+\min \left(H_{\mathrm{real}}\left(E\right)\right)}<\mathrm{\&phi;}$ （0<φ<1）

If do not meet, turn back to S12 and continue iteration; If satisfied to S15;

S15: change modifying factor f, then return and carry out S12 continuation iteration, until satisfy condition

Δ E=E wherein _max-E _minenergy range for system;

S16: once satisfy condition

original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration, s=s+1, modifying factor is $\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{new}})+1);$

S17: continue S16 iteration, until meet procedure termination condition

by the S of all processes (E) is averaging

? and then obtain the relative density of states of protein system

Further, describedly from process, comprise the following steps:

S21: logarithm S (the E)=lng (E)=0 of the density of states of initialization system, S _tmp(E)=lng _tmp(E)=0, histogram H (E)=0, H _tmp(E)=0 (E _min≤ E≤E _max), s=1, modifying factor df| _e=(κ Θ (E ₀-E)+1) lnf, wherein Θ (E ₀-E) be Heaviside piecewise function, κ, E ₀, f is the parameter that is relevant to model;

S22：t=1；

S23: from process, original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration, t=t+1, s=s+1, modifying factor is df| _e=(κ Θ (E ₀-E)+1) lnf;

Step S23 circulation tmax time;

S24: all process intercommunications, from process by S _tmpand H (E) _tmp(E) send to host process, then receive the overall S (E) and the H (E) that through host process, calculate and upgrade original S (E) and H (E), by S _tmpand H (E) _tmp(E) be initialized as 0, try to achieve H _real(E)=H (E) * [1+ κ Θ (E ₀-E)], the mild condition of judgement histogram:

$\frac{\max \left(H_{\mathrm{real}}\left(E\right)\right)-\min \left(H_{\mathrm{real}}\left(E\right)\right)}{\max \left(H_{\mathrm{real}}\left(E\right)\right)+\min \left(H_{\mathrm{real}}\left(E\right)\right)}<\mathrm{\&phi;},$ 0< φ <1 wherein

If do not meet, return and carry out S22 continuation iteration; If satisfied perform step S25;

S25 changes modifying factor f, then returns and carry out S22 continuation iteration, until satisfy condition

Δ E=E wherein _max-E _minenergy range for system;

S26: once condition

meet, original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration;

S27: continue step S26 iteration, until meet procedure termination condition

wherein

Further, in step S13, according to Metropolis criterion, determine that the received probability of new configuration further comprises:

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration:

S(E _new)=S(E _new)+df| _E,H(E _new)=H(E _new)+1；

Otherwise:

S(E _old)=S(E _old)+df| _E,H(E _old)=H(E _old)+1。

Further, in step S23, according to Metropolis criterion, determine that the received probability of new configuration further comprises:

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration:

S(E _new)=S(E _new)+df| _E,H(E _new)=H(E _new)+1，

S _tmp(E _new)=S _tmp(E _new)+df| _E,H _tmp(E _new)=H _tmp(E _new)+1；

Otherwise:

S(E _old)=S(E _old)+df| _E,H(E _old)=H(E _old)+1，

S _tmp(E _old)=S _tmp(E _old)+df| _E,H _tmp(E _old)=H _tmp(E _old)+1。

Further, in step S16 and S26, according to Metropolis criterion, determine that the received probability of new configuration further comprises:

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration:

$S\left(E_{\mathrm{new}}\right)=S\left(E_{\mathrm{new}}\right)+\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{new}})+1),H\left(E_{\mathrm{new}}\right)=H\left(E_{\mathrm{new}}\right)+1;$

Otherwise:

$S\left(E_{\mathrm{old}}\right)=S\left(E_{\mathrm{old}}\right)+\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{old}})+1),H\left(E_{\mathrm{old}}\right)=H\left(E_{\mathrm{old}}\right)+1.$

Further, in step S15 and S25, change modifying factor f mode is:

First carry out continuously the f=f of N iteration ^α(0< α <1), then carry out 1 iteration

Repeatedly repeat aforesaid way.

Compare with classical Wang-Landau algorithm, the present invention uses the modifying factor stage by stage with overall regeneration characteristics, can accelerate search and analog rate; Utilize the mode of renewal modifying factor flexibly of annealing (Annealing) mechanism can improve simulation precision and speed; Adopt the parallel mode of this algorithm can greatly accelerate search and analog rate.Based on method provided by the invention, can efficient analysis and the folding whole thermodynamic process of Study on Protein, and then protein folding procedure is explored and studied.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of the protein thermodynamics analytical approach based on Monte Carlo simulation of the present invention.

Fig. 2 is the parallel algorithm process flow diagram of the density of states of simulation and calculating protein system.

Embodiment

In order to make object of the present invention, technical scheme and advantage clearer, below in conjunction with drawings and Examples, the present invention is further elaborated.Should be appreciated that specific embodiment described herein, only in order to explain the present invention, is not intended to limit the present invention.

In addition,, in each embodiment of described the present invention, involved technical characterictic just can not combine mutually as long as do not form each other conflict.

The invention provides a kind of protein thermodynamics analytical approach based on Monte Carlo simulation, with reference to figure 1, the method is shown and mainly comprises two basic steps: steps A is determined protein energy model; And step B simulation and the calculating protein system density of states.

In steps A, can adopt ECEPP energy force field model and angle coordinate system.Wherein the expression-form in the ECEPP energy field of force is:

E _ECEPP=E _C+E _LJ+E _HB+E _Tor

Wherein,

the coulomb acting force between two electric charges, r _ijrepresent the distance between atom i and j,

the Lan Na-Jones acting force between two atoms,

hyarogen-bonding, E _tor=∑ _lu _l(1 ± cos (n _lξ _l)) be dihedral turning effort power, ξ _ll dihedral.Because the ECEPP energy field of force adopts angle coordinate system, so counting yield is higher than the energy field of force based on cartesian coordinate system.

For the ease of computer simulation emulation, can also, to carrying out discretize processing between used protein energy range, if get k energy bin interval value, need [E _min, E _max] on average divide k bin interval, with each interval middle energy value, represent this energy interval value.

In step B, by the MPI concurrent program algorithm of principal and subordinate's process mode, the density of states of simulation and calculating protein system.In N minute process of described principal and subordinate's process mode, a minute process 1 is host process, and within all the other minute, process is subprocess.

With reference to figure 2, the parallel algorithm process flow diagram of the density of states of simulation and calculating protein system is shown.As can be seen from Fig. 2, host process comprises the following steps:

S11: logarithm S (the E)=lng (E)=0 of the density of states of initialization system, histogram H (E)=0(E _min≤ E≤E _max), s=1, modifying factor df| _e=(κ Θ (E ₀-E)+1) lnf, wherein Θ (E ₀-E) be Heaviside piecewise function, κ, E ₀, f is the parameter that is relevant to model, as desirable κ=5, E ₀=-2, f=e.

S12：t=1。

S13: in host process, original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration: (referred to as MCS-f step, t=t+1, s=s+1, modifying factor is df| _e=(κ Θ (E ₀-E)+1) lnf)

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration: S (E _new)=S (E _new)+df| _e, H (E _new)=H (E _new)+1;

Otherwise: S (E _old)=S (E _old)+df| _e, H (E _old)=H (E _old)+1.

Tmax (as 100 times) MCS-f step of step S13 circulation.

S14: all process intercommunications, host process is collected all S from process _tmpand H (E) _tmp(E) and accumulation calculating go out overall S (E) and H (E), the overall situation S (E)=S (E)+all S from process _tmp(E), H (E)=H (E)+all H from process of the overall situation _tmp(E).Then being broadcast to of overall S (E) and H (E) is all from process, then try to achieve H _real(E)=H (E) * [1+ κ Θ (E ₀-E)], finally judge the mild condition of histogram:

$\frac{\max \left(H_{\mathrm{real}}\left(E\right)\right)-\min \left(H_{\mathrm{real}}\left(E\right)\right)}{\max \left(H_{\mathrm{real}}\left(E\right)\right)+\min \left(H_{\mathrm{real}}\left(E\right)\right)}<\mathrm{\&phi;}$ (0< φ <1 as desirable in φ 0.2)

If do not meet, turn back to S12 and continue iteration; If satisfied to S15.

S15: change modifying factor f, then turn back to S12 continuation iteration, until satisfy condition

(Δ E=E wherein _max-E _minenergy range for system).Changing modifying factor f mode can be: first carry out continuously the f=f of N iteration ^α(0< α <1), then carry out 1 iteration

repeatedly repeat this mode.

S16: once condition

meet, original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration: (referred to as MCS-s step, s=s+1, modifying factor is

)

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration: $S\left(E_{\mathrm{new}}\right)=S\left(E_{\mathrm{new}}\right)+\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{new}})+1),H\left(E_{\mathrm{new}}\right)=H\left(E_{\mathrm{new}}\right)+1;$

Otherwise: $S\left(E_{\mathrm{old}}\right)=S\left(E_{\mathrm{old}}\right)+\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{old}})+1),H\left(E_{\mathrm{old}}\right)=H\left(E_{\mathrm{old}}\right)+1.$

S17: continue S16 and continue iteration, until meet procedure termination condition

( as

desirable 0.0001), so just can be by the S of all processes (E) be averaging

?

and then obtain the relative density of states of protein system $g\left(E\right)=e^{S_{\mathrm{real}}\left(E\right)}.$

from process, comprise the following steps:

S21: logarithm S (the E)=lng (E)=0 of the density of states of initialization system, S _tmp(E)=lng _tmp(E)=0, histogram H (E)=0, H _tmp(E)=0 (E _min≤ E≤E _max), s=1, modifying factor df| _e=(κ Θ (E ₀-E)+1) lnf, wherein Θ (E ₀-E) be Heaviside piecewise function, κ, E ₀, f is the parameter that is relevant to model, as desirable κ=5, E ₀=-2, f=e.

S22：t=1。

S23: from process, original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration: (referred to as MCS-f step, t=t+1, s=s+1, modifying factor is df| _e=(κ Θ (E ₀-E)+1) lnf)

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration:

S(E _new)=S(E _new)+df| _E,H(E _new)=H(E _new)+1，

S _tmp(E _new)=S _tmp(E _new)+df| _E,H _tmp(E _new)=H _tmp(E _new)+1；

Otherwise:

S(E _old)=S(E _old)+df| _E,H(E _old)=H(E _old)+1，

S _tmp(E _old)=S _tmp(E _old)+df| _E,H _tmp(E _old)=H _tmp(E _old)+1。

Tmax (as 100 times) MCS-f step of step S23 circulation.

S24: all process intercommunications, from process by S _tmpand H (E) _tmp(E) send to host process, then receive the overall S (E) and the H (E) that through host process, calculate and upgrade original S (E) and H (E), then by S _tmpand H (E) _tmp(E) be initialized as 0, then try to achieve H _real(E)=H (E) * [1+ κ Θ (E ₀-E)], finally judge the mild condition of histogram:

$\frac{\max \left(H_{\mathrm{real}}\left(E\right)\right)-\min \left(H_{\mathrm{real}}\left(E\right)\right)}{\max \left(H_{\mathrm{real}}\left(E\right)\right)+\min \left(H_{\mathrm{real}}\left(E\right)\right)}<\mathrm{\&phi;}$ (0< φ <1 as desirable in φ 0.2)

If do not meet, turn back to S22 and continue iteration; If satisfied to S25.

S25: change modifying factor f, then turn back to S22 continuation iteration, until satisfy condition

(Δ E=E wherein _max-E _minenergy range for system).Changing modifying factor f mode can be: first carry out continuously the f=f of N iteration ^α(0< α <1), then carry out 1 iteration

repeatedly repeat this mode.

S26: once satisfy condition

original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration: (referred to as MCS-s step, s=s+1, modifying factor is

)

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration: $S\left(E_{\mathrm{new}}\right)=S\left(E_{\mathrm{new}}\right)+\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{new}})+1),H\left(E_{\mathrm{new}}\right)=H\left(E_{\mathrm{new}}\right)+1;$

Otherwise: $S\left(E_{\mathrm{old}}\right)=S\left(E_{\mathrm{old}}\right)+\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{old}})+1),H\left(E_{\mathrm{old}}\right)=H\left(E_{\mathrm{old}}\right)+1.$

S27: continue S26 iteration, until meet procedure termination condition

(

as

desirable 0.0001).

The present invention uses the protein thermodynamics analytical approach based on Monte Carlo simulation, can effectively simulate and calculate the protein system density of states, thereby just can calculate the many canonical ensemble thermodynamic quantities in broad temperature range, the whole thermodynamic process that so just energy Study on Protein folds, and then protein folding procedure is explored and studied.

Compare with classical Wang-Landau algorithm, method of the present invention has the modifying factor stage by stage of overall regeneration characteristics can accelerate search and analog rate; Utilize the mode of renewal modifying factor flexibly of annealing (Annealing) mechanism can improve simulation precision and speed; Adopt the parallel mode of this algorithm can greatly accelerate search and analog rate.Can efficient analysis and the folding whole thermodynamic process of Study on Protein based on this, and then protein folding procedure is explored and studied, and can promote the use of the high efficiency technical in a lot of fields.

One of ordinary skill in the art will appreciate that all or part of step in the whole bag of tricks of embodiment is to come the hardware that instruction is relevant to complete by program, this program can be stored in a computer-readable recording medium, storage medium can comprise: ROM (read-only memory) (ROM, Read Only Memory), random access memory (RAM, Random Access Memory), disk or CD etc.

The above the specific embodiment of the present invention, does not form limiting the scope of the present invention.Various other corresponding changes and distortion that any technical conceive according to the present invention has been done, all should be included in the protection domain of the claims in the present invention.

Claims (10)

1. the protein thermodynamics analytical approach based on Monte Carlo simulation, is characterized in that, comprises step

A: determine protein energy model; And

B: simulation and the calculating protein system density of states.

2. protein thermodynamics analytical approach as claimed in claim 1, wherein, steps A further comprises:

Adopt ECEPP energy force field model and angle coordinate system, the expression-form in the ECEPP energy field of force is:

E _ECEPP=E _C+E _LJ+E _HB+E _Tor

Wherein

the coulomb acting force between two electric charges, r _ijrepresent the distance between atom i and j, the Lan Na-Jones acting force between two atoms,

hyarogen-bonding, E _tor=∑ _lu _l(1 ± cos (n _lξ _l)) be dihedral turning effort power, ξ _ll dihedral.

3. protein thermodynamics analytical approach as claimed in claim 1, wherein, steps A further comprises:

To carrying out discretize processing between used protein energy range, if get k energy bin interval value, need [E _min, E _max] on average divide k bin interval, with each interval middle energy value, represent this energy interval value.

4. protein thermodynamics analytical approach as claimed in claim 1, wherein, step B further comprises:

By the MPI concurrent program algorithm of principal and subordinate's process mode, the density of states of simulation and calculating protein system.

5. protein thermodynamics analytical approach as claimed in claim 4, wherein, in N minute process of described principal and subordinate's process mode, a minute process 1 is host process, within all the other minute, process is subprocess.

6. protein thermodynamics analytical approach as claimed in claim 5, wherein, described host process comprises step:

S11: logarithm S (the E)=lng (E)=0 of initialization protein system Density function, histogram H (E)=0(E _min≤ E≤E _max), s=1, modifying factor df| _e=(κ Θ (E ₀-E)+1) lnf, wherein Θ (E ₀-E) be Heaviside piecewise function, κ, E ₀, f is the parameter that is relevant to model;

S12：t=1；

S13: in host process, original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration, t=t+1, s=s+1, modifying factor is df| _e=(κ Θ (E ₀-E)+1) lnf;

Step S13 circulation tmax time;

S14: all process intercommunications, host process is collected all S from process _tmpand H (E) _tmp(E) and accumulation calculating go out overall S (E) and H (E), the overall situation S (E)=S (E)+all S from process _tmp(E), H (E)=H (E)+all H from process of the overall situation _tmp(E), being broadcast to of overall S (E) and H (E) is all from process, try to achieve H _real(E)=H (E) * [1+ κ Θ (E ₀-E)], the mild condition of judgement histogram:

$\frac{\max \left(H_{\mathrm{real}}\left(E\right)\right)-\min \left(H_{\mathrm{real}}\left(E\right)\right)}{\max \left(H_{\mathrm{real}}\left(E\right)\right)+\min \left(H_{\mathrm{real}}\left(E\right)\right)}<\mathrm{\&phi;}$ （0<φ<1）

If do not meet, turn back to S12 and continue iteration; If satisfied to S15;

S15: change modifying factor f, then return and carry out S12 continuation iteration, until satisfy condition

Δ E=E wherein _max-E _minenergy range for system;

S16: once satisfy condition

original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration, s=s+1, modifying factor is $\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{new}})+1);$

S17: continue S16 iteration, until meet procedure termination condition

by the S of all processes (E) is averaging

?

and then obtain the relative density of states of protein system

7. protein thermodynamics analytical approach as claimed in claim 6, wherein, describedly comprises step from process:

S21: logarithm S (the E)=lng (E)=0 of the density of states of initialization system, S _tmp(E)=lng _tmp(E)=0, histogram H (E)=0, H _tmp(E)=0 (E _min≤ E≤E _max), s=1, modifying factor df| _e=(κ Θ (E ₀-E)+1) lnf, wherein Θ (E ₀-E) be Heaviside piecewise function, κ, E ₀, f is the parameter that is relevant to model;

S22：t=1；

S23: from process, original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration, t=t+1, s=s+1, modifying factor is df| _e=(κ Θ (E ₀-E)+1) lnf;

Step S23 circulation tmax time;

S24: all process intercommunications, from process by S _tmpand H (E) _tmp(E) send to host process, then receive the overall S (E) and the H (E) that through host process, calculate and upgrade original S (E) and H (E), by S _tmpand H (E) _tmp(E) be initialized as 0, try to achieve H _real(E)=H (E) * [1+ κ Θ (E ₀-E)], the mild condition of judgement histogram:

$\frac{\max \left(H_{\mathrm{real}}\left(E\right)\right)-\min \left(H_{\mathrm{real}}\left(E\right)\right)}{\max \left(H_{\mathrm{real}}\left(E\right)\right)+\min \left(H_{\mathrm{real}}\left(E\right)\right)}<\mathrm{\&phi;},$ 0< φ <1 wherein

If do not meet, return and carry out S22 continuation iteration; If satisfied perform step S25;

S25 changes modifying factor f, then returns and carry out S22 continuation iteration, until satisfy condition

Δ E=E wherein _max-E _minenergy range for system;

S26: once condition meet, original configuration is carried out to random fluctuation, produce new configuration, calculating energy E _new, according to Metropolis criterion, determine the received probability of new configuration;

S27: continue step S26 iteration, until meet procedure termination condition

wherein

8. protein thermodynamics analytical approach as claimed in claim 7, wherein,

In step S13, according to Metropolis criterion, determine that the received probability of new configuration further comprises:

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration:

S(E _new)=S(E _new)+df| _E,H(E _new)=H(E _new)+1；

Otherwise:

S(E _old)=S(E _old)+df| _E,H(E _old)=H(E _old)+1；

In step S23, according to Metropolis criterion, determine that the received probability of new configuration further comprises:

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration:

S(E _new)=S(E _new)+df| _E,H(E _new)=H(E _new)+1，

S _tmp(E _new)=S _tmp(E _new)+df| _E,H _tmp(E _new)=H _tmp(E _new)+1；

Otherwise:

S(E _old)=S(E _old)+df| _E,H(E _old)=H(E _old)+1，

S _tmp(E _old)=S _tmp(E _old)+df| _E,H _tmp(E _old)=H _tmp(E _old)+1。

9. protein thermodynamics analytical approach as claimed in claim 8, wherein, in step S16 and S26, according to Metropolis criterion, determine that the received probability of new configuration further comprises:

$P(\mathrm{old}\&RightArrow;\mathrm{new})=\min (1,e^{-[S\left(E_{\mathrm{new}}\right)-S\left(E_{\mathrm{old}}\right)]})$

If accept new configuration:

$S\left(E_{\mathrm{new}}\right)=S\left(E_{\mathrm{new}}\right)+\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{new}})+1),H\left(E_{\mathrm{new}}\right)=H\left(E_{\mathrm{new}}\right)+1;$

Otherwise:

$S\left(E_{\mathrm{old}}\right)=S\left(E_{\mathrm{old}}\right)+\frac{\mathrm{\&Delta;E}}{s}(\mathrm{\&kappa;\&Theta;}(E_0-E_{\mathrm{old}})+1),H\left(E_{\mathrm{old}}\right)=H\left(E_{\mathrm{old}}\right)+1.$

10. protein thermodynamics analytical approach as claimed in claim 9, wherein, in step S15 and S25, changes modifying factor f mode and is:

First carry out continuously the f=f of N iteration ^α(0< α <1), then carry out 1 iteration

Repeatedly repeat aforesaid way.

Images