Free energy difference estimation method and simulation apparatus

In a free energy difference estimation method, the partial free energy difference indicates a bound state between a target compound and a first candidate compound bindable to the target compound, and is stored in a storage part for each value of a binding constant. A change region, in which a partial free energy difference is equal to or greater than a predetermined change value, is specified in a region of the binding constant. The partial free energy difference of a second candidate compound in the change region is interpolated based on the partial free energy difference acquired in each of a previous region and a following region by using an approximation function.

CROSS-REFERENCE TO RELATED APPLICATIONS

This patent application is based upon and claims the benefit of priority of the prior Japanese Patent Application No. 2011-055828 filed on Mar. 14, 2011, the entire contents of which are incorporated herein by reference.

FIELD

The embodiments discussed herein are related to a free energy difference estimation method and a simulation apparatus for calculating a free energy difference related to a binding between a macromolecular compound and a candidate compound in a solvent at higher speed.

BACKGROUND

Conventionally, molecular design is performed by using a molecular simulation. A problem of the molecular simulation is mainly the calculation amount, and a means has been provided to effectively process the calculation amount at higher speed.

For example, an electric charge determination method or the like is provided in a basis function, which is continuous for at least more than second order derivative, is applied. In the electric charge determination method, when either one grid point (or a particle) in a pair of grid points (or a pair of particles) is overlapped with a grid point in a roughly level grid which is at one level higher than the current level, an electric charge is determined for a grid point of a higher level so that an energy value of a pair in the higher level always corresponds to a proper energy value of the pair in the current level.

In drug development using a computer (computational drug development), calculation time and calculation resources cost less than a laboratory experiment. Therefore, computational drug development is simply used to comprehend a microscopic phenomenon which is not acquired by the laboratory experiment. If attempting to acquire data which are in no way inferior to the experiment, five days may be spent with approximately 320 Central. Processing Units (CPUs) for one reagent.Patent Document 1 JP2007-80044A

SUMMARY

According to one aspect of the embodiment, there is provided a free energy difference estimation method performed in a computer to estimate a free energy difference between compounds, the method including specifying, by the computer, a change region, in which a partial free energy difference is equal to or greater than a predetermined change value, in a region of a binding constant, by referring to the partial free energy difference for each value of the binding constant indicating a bound state between a target compound and a first candidate compound bindable to the target compound, the partial free energy difference stored in a storage part; and interpolating, by the computer, the partial free energy difference of a second candidate compound in the change region, based on the partial free energy difference acquired in each of a previous region and a following region by using an approximation function.

According to other aspects of the embodiment, there may be provided a simulation apparatus for performing the free energy difference estimation method, and a non-transitory computer-readable recording medium recorded with a program which, when executed by a computer, causes the computer to perform a free energy difference estimation process.

DESCRIPTION OF EMBODIMENTS

In the following, embodiments of the present invention will be described with reference to the accompanying drawings. First, a binding strength of a drug candidate compound for a target protein with noncovalent binding will be described with reference toFIG. 1andFIG. 2. The target protein is regarded as a biological functional agent concerning disease. A candidate compound is regarded as a compound which strongly binds to the target protein, and may be a compound which is considered as a drug for binding to the target protein and for blocking a function of the target protein.

FIG. 1is a diagram for explaining the binding strength. InFIG. 1, an unbound state2aand a fully bound state2bin water5(in a solvent) are illustrated. In the unbound state2a, a target protein3is separated from a candidate compound4in the water5. In the fully bound state2b, the target protein3is bound to the candidate compound4in the water5. A binding strength ΔG between the target protein3and the candidate compound4is represented by a relative free energy difference (hereinafter, simply called energy difference) between the unbound state2aand the fully bound state2b.

The binding strength ΔG is acquired by adding a first energy difference ΔG1and a second energy difference ΔG2. The first energy difference ΔG1indicates a partial energy difference between the unbound state2aand a state2cin which the candidate compound4in the water5is extracted into vacuum6. The second energy difference ΔG2indicates another partial energy difference between the fully bound state2band the state2c.

Approximately 100 candidate compounds, which are narrowed down beforehand, are used as drug candidates in which the binding strength is simulated by a drug development simulator100applying Molecular Dynamics (MD). For one drug candidate, a simulation time of the first energy difference ΔG1is approximately one day. On the contrary, the simulation time is approximately two weeks until a value of the second energy difference ΔG2is converged.

FIG. 2is a diagram for explaining the second energy difference ΔG2. InFIG. 2, λ indicates a binding constant. When λ=0 (zero), the target protein3and the candidate compound4are in the fully bound state2b. When λ=1, the target protein3and the candidate compound4are in the state2c. In other words, λ is used as an index for representing strength of interaction between the candidate compound4and atoms surrounding the candidate compound4. When λ=0, an inherent interaction is indicated. When λ=1, a state in which there is no interaction is indicated. At intervals λ from λ0to λnin a range of 0≦λ≦1, an energy difference Δg(λi,i+1) kcal/mol (i=0, 1, . . . n−1) is calculated by using the following calculation. A total value becomes the second energy difference ΔG2.

FIG. 3is a graph illustrating a relationship between the binding strength and a calculation amount at the intervals λ. InFIG. 3, a horizontal axis indicates λ, a left vertical axis indicates the energy difference Δg(λi,i+1) at the intervals λ, and a right vertical axis indicates the calculation amount. The calculation amount may be either days or time spent for the calculation. Also, points8aindicate the energy differences Δg respective to the intervals λ, and a dashed line8bindicates the calculation amount.

InFIG. 3, when a value of λ is smaller than approximately 0.6, in a case of being closer to the fully bound state2b, energy between the target protein3and the candidate compound4changes less and therefore is easily converted. Thus, approximately only one day is needed for the calculation amount. On the other hand, in a region7where the value of λ is approximately 0.7 to 0.9 which is closer to the unbound state2a, since the energy difference Δg greatly changes, the interval may be given by λ precisely in order to acquire an accurate calculation result. Also, in the region7, more time is spent for convergence. Thus, a calculation time becomes significantly longer.

Accordingly, the inventor has found the fact that it is possible to acquire plots depicting a smooth line, without a specific point, by representing the binding strength with a differential value Δg/Δλ, regarding the binding strength at the intervals λ illustrated inFIG. 3, and the plot precisely indicates 0 where λ=1 (FIG. 4andFIG. 5).

FIG. 4is a diagram illustrating a relationship between the binding strength and the calculation amount by using the differential value. InFIG. 4, a horizontal axis indicates λ, a left vertical axis indicates a differential value Δg(λi,i+1)/Δλi,i+1of the energy difference Δg at the intervals λ, and a right vertical axis indicates the calculation amount.

Referring toFIG. 4, points8′ a indicate differential values Δg/Δλ of the energy difference Δg at the intervals λ, and form a smooth line without a specific point. Also, the plot precisely indicates 0 where λ=1. In a region exceeding the differential value 0 (zero), the calculation amount rapidly increases. The region where the calculation amount rapidly increases is regarded as a region where a calculation value is difficult to converge.

Accordingly, in the embodiment, a region of λ is divided into three regions by using the differential value, values in a region7′, where the calculation value is difficult to converge, are estimated from values in a previous region and a following region.

FIG. 5is a diagram for explaining a calculation method for each of the divided regions according to the embodiment. InFIG. 5, the horizontal axis and vertical axes are the same as those inFIG. 4, and the explanation thereof will be omitted. InFIG. 5, λ1indicates a value of λ where the differential value becomes 0 (zero), and λ2indicates another value of λ where the differential value becomes less than a predetermined value after a peak value. In the embodiment, the region from λ=0 to λ=1 is divided into a region A (0≦λ≦λ1), and a region B (λ1<λ<λ2), and a region C (λ2≦λ≦1). The region A is from λ=0 to the differential value of the energy difference Δg=0 (zero). The region B is where the differential value of the energy difference Δg=0 (zero), and continues until λ=λ2where the differential value becomes less than the predetermined value after the peak. The region C is from λ=λ2where the differential value becomes less than the predetermined value after the peak of the differential value of the energy difference Δg to λ=1.

In the region A, the second energy difference ΔG2between the fully bound state2band the state2cis calculated. In the region C, the first energy difference ΔG1between the unbound state2aand the state2c, which is calculated beforehand, is substituted, since the first energy difference ΔG1indicates the same value as the second energy difference ΔG2. In the region B, the energy difference Δg is approached and interpolated by using an approximation function such as a three dimensional spline or the like based on values in the region A and region C.

In the embodiment, a simulation is preformed at the intervals λ in all regions A, B, and C for a first candidate compound, and the differential values Δg/Δλ of the energy difference Δg at the intervals λ. Points8′aindicate results (the second energy differences ΔG2) from the simulation at the intervals λ in all regions A, B, and C.

Calculations depending on the regions A, B, and C are performed from a next candidate compound. Shaded points9aindicate the second energy difference ΔG2in the region A. Shaded points9cindicate the first energy difference ΔG1between the unbound state2aand the state2cof the candidate compound, with which the energy difference Δg is substituted. Shaded points9bindicate values interpolated by using the three dimensional spline or the like from values of the region A and the region B.

As described above, the drug development simulator100for estimating a binding strength of the candidate compound for the target protein includes a hardware configuration as illustrated inFIG. 6.FIG. 6is a block diagram illustrating the hardware configuration of the drug development simulator100.

As illustrated inFIG. 6, the drug development simulator100is regarded as a terminal controlled by a computer, and includes a processor such as a Central Processing Unit (CPU)11, a memory device12, a display device13, an output device14, an input device15, a communication device16, a storage device17, and a drive device18, which are mutually connected via a bus B.

The CPU11controls the drug development simulator100in accordance with a program stored in the memory device12. The memory device12includes a Random Access Memory (RAM), a Read-Only Memory (ROM), and the like, and stores the program executed by the CPU11, data for the process by the CPU11, data acquired in the process by the CPU11, and the like. Also, a part of an area of the memory device12is assigned as a working area used in processing by the CPU11.

The display device13is used to display various information items for control of the CPU11. The output device14includes a printer or the like, and is used to output various information items in response to an instruction of a user. The input device15includes a mouse, a keyboard, and the like, and is used by the user to input various information items for the process conducted by the drug development simulator100. For example, the communication device16is connected to the Internet, a Local Area Network (LAN), and the like, and is used to control communications with an external device. For example, the storage device17includes a hard disk unit, and stores the program and data for conducting various processes.

For example, the program for realizing the process by the drug development simulator100are provided to the drug development simulator100by a recording medium19such as the Compact Disc Read Only Memory (CD-ROM) or the like. That is, when the recording medium19storing the program is set to the drive device18, the drive device18reads out the program from the recording medium19. The program read from the recording medium19is installed onto the storage device17through the bus B. Accordingly, when the program is activated, the CPU11initiates the process in accordance with the program installed onto the storage device17. It should be noted that a medium for storing the program is not limited to the CD-ROM, and any kind of a computer-readable recording medium may be used. As the computer-readable recording medium, a portable recording medium such as a Digital Versatile Disk (DVD), a Universal Serial Bus (USB) memory, a semiconductor memory such as a flash memory, or the like may be used as well as the CD-ROM.

In addition, the program for realizing the process of the drug development simulator100may be provided from the external device through the communication device16.

The drug development simulator100according to the embodiment includes a functional configuration as illustrated inFIG. 7.FIG. 7is a diagram illustrating the functional configuration of the drug development simulator100according to the embodiment. InFIG. 7, the drug development simulator100includes a molecular structure input part31, a first energy difference calculation part32, a second energy difference calculation part33, a region division part34, a region A calculation part35, a region C substitution part36, a region B interpolation part37, a binding strength calculation part38, and the like. The parts31through38are realized by the CPU11executing corresponding programs.

Also, in the drug development simulator100, target protein data40, a candidate compound list41, a Δg1list42, a Δg1differential value list42-2for a region C, a Δg2list43, a Δg2differential value list43-2for a region A, a Δg2reference value44, boundary λ data45for λ1and λ2, and binding strength data46are stored in a storage part50. The storage part50corresponds to a part of the memory device12and/or the storage device17.

The molecular structure input part31inputs, by the user, a molecular structure of the target protein and the candidate compound to the drug development simulator100. The molecular structure input part31is regarded as a process part which stores the target protein data40indicating the molecular structure of the target protein in the storage part50, and also stores the candidate compound list41listing the molecular structure for identification information of each of the candidate compounds in the storage part50.

The first energy difference calculation part32is regarded as a process part which refers to the candidate compound list41, and calculates the first energy difference Δg1at the intervals λ based on a macromolecular structure of the target protein indicated by the target protein data40and a molecular structure of each of the candidate compounds. The first energy difference Δg1for each of the intervals λ corresponds to the identification information of the candidate compound and is stored in the Δg1list42in the storage part50.

The second energy difference calculation part33is regarded as a process part which calculates the second energy difference Δg2at the intervals λ in the entire regions A, B, and C, based on the molecular structure of a first candidate compound and the target protein data40. The second energy difference Δg2calculated at the intervals λ corresponds to the identification information of the first candidate compound and is stored in the Δg2list43. Also, the second energy difference calculation part33calculates the differential value of the second energy difference Δg2of the first candidate compound at the intervals λ, and stores the differential value as a Δg2reference value44in the storage part50.

The region division part34refers to the Δg2reference value44, specifies λ1where the differential value is 0 (zero) and λ2, and stores the boundary λ data45indicating λ1and λ2in the storage part50.

The region A calculation part35is regarded as a process part which calculates the second energy difference Δg2in the region A (0≦λ≦λ1) based on the molecular structure for each of the candidate compounds other than the first candidate compound, which is acquired by referring to the candidate compound list41, and the target protein data40. The second energy difference Δg2is calculated and additionally stored into the Δg2list43in the storage part50, for each of the intervals λ in the region A, so as to correspond to the identification information of the candidate compound. Also, the region A calculation part35calculates a differential value of the second energy difference Δg2at the intervals λ in the region A, and stores the differential value in the Δg2differential value list43-2.

The region C substitution part36is regarded as a process part which acquires the first energy difference Δg1between the candidate compound and the target protein data40, which is calculated by the first energy difference calculation part32, from the Δg1list42, and substitutes with the first energy difference Δg1for the region C. The first energy difference Δg1acquired from the Δg1list42corresponds to the identification information of the candidate compound and is additionally stored in the Δg2list43for each of the intervals λ in the region C. Also, the region C substitution part36calculates the differential value of the first energy difference Δg1at the intervals λ in the region C, and stores the differential value for each interval in the Δg1differential value list42-2.

The region B interpolation part37is regarded as a process part which refers to the differential value of the second energy difference Δg2for each of the intervals λ in the region A from the Δg2differential value list43-2and the differential value of the second energy difference Δg1for each of the intervals λ in the region C from the Δg1differential value list42-2, and interpolates the differential value of the second energy difference Δg2in the region B by using the three dimensional spline or the like. The region B interpolation part37estimates and interpolates the differential values in the region B based on the differential values at λ1and λ2of boundaries of the region B.

Also, the region B interpolation part37inversely calculates the second energy difference Δg2at the intervals λ in the region B from the differential value thereof being interpolated (multiplies the differential value by a respective value of λ), and additionally stores the calculated second energy difference Δg2, which corresponds to the identification information of the candidate compound, to the Δg2list43at the intervals λ in the region B.

The binding strength calculation part38is regarded as a process part which calculates the binding strength. The binding strength calculation part38acquires the first energy difference ΔG1by calculating a total value for the first energy difference Δg1at the intervals λ by referring to the Δg1list42for each of the candidate compounds. Moreover, the binding strength calculation part38acquires the second energy difference ΔG1by calculating a total value for the second energy difference Δg2at the intervals λ by referring to the Δg2list43for each of the candidate compounds. Furthermore, the binding strength calculation part38calculates a total of the first energy difference Δg1and the second energy difference Δg2, to acquire the binding strength. The binding strength acquired by adding the first energy difference Δg1and the second energy difference Δg2is stored as the binding strength data46corresponding to the identification information of the candidate compound in the storage part50. When the binding strengths are calculated for all candidate compounds listed in the candidate compound list41, the binding strengths in the binding strength data46may be sorted in a descending order of the binding strength or the like and may be displayed at the display device13in response to a request of the user.

FIG. 8is a flowchart for explaining the simulation for acquiring the binding strength by the drug development simulator100. InFIG. 8, in the drug development simulator100, the first energy difference calculation part32reads out a macromolecular structure of the target protein from the target protein data40stored in the storage part50, and reads out the molecular structure of the candidate compound from the candidate compound list41(step S11), and calculates the first energy difference Δg1at the intervals λ (step S12). The first energy difference calculation part32stores the first energy difference Δg1corresponding to the identification information of the candidate compound at the intervals λ in the Δg1list42in the storage part50.

After that, it is determined whether or not a current candidate compound is the first candidate compound (step S13). That is, it may be checked whether the current candidate compound is a candidate compound stored in a beginning of the candidate compound list41. In a case of the first candidate compound, the second energy difference calculation part33calculates the second energy difference Δg2at the intervals λ in the entirety of the regions A, B, and C, based on the molecular structure of the first candidate compound and the target protein data40(step S14). The second energy difference calculation part33outputs and stores the second energy difference Δg2at the intervals λ in the Δg2list43in which the second energy difference Δg2corresponds to the identification information of the first candidate compound. Also, the second energy difference calculation part33calculates the differential value of each second energy difference Δg2in the entirety for the first candidate compound in the entirety of the regions A, B, and C, and stores the differential value of the second energy difference Δg2calculated at the intervals λ to the Δg2reference value44in the storage part50.

Next, the region division part34refers to the Δg2reference value43stored in the storage part50, specifies λ1where the differential value of the second energy difference Δg2indicates 0 (zero) and λ2based on the differential values of the second energy difference Δg2at the intervals λ, and stores the boundary λ data45indicating λ2and λ2(step S15).

On the other hand, in the step S13, when it is determined that the current candidate compound is the first candidate compound, the region A calculation part35calculates the second energy difference Δg2at the intervals λ in the region A, based on the molecular structure of the candidate compound which the first energy difference calculation part32reads out from the candidate compound list41and the macromolecular structure of the target protein indicated in the target protein data40(step S18). The region A calculation part35acquires λ1from the boundary A data45and specifies the region A (0≦λ≦λ1). Then, the region A calculation part35calculates the second energy difference Δg2at the intervals λ in the specified region A, and additionally stores the second energy difference Δg2, which is corresponded to, at the intervals λ in the region A, in Δg2list43.

Also, the region A calculation part35calculates the differential value of the second energy difference Δg2at the intervals λ in the region A, and stores the differential value calculated at the intervals λ in the Δg2differential value list43-2.

Next, the region C substitution part36acquires the first energy difference Δg1between the candidate compound and the target protein data40, which is calculated by the first energy difference calculation part32for the region C, and substitutes with the acquired first energy difference Δg1for the region C (step S19). The region C substitution part36specifies the region C (λ2≦λ≦1) by acquiring λ2from the boundary λ data45, and acquires the first energy difference Δg1of each of the intervals λ for the specified region C from the Δg1list42. The acquired first energy difference Δg1is corresponded to the identification information of the candidate compound and is stored in the Δg2list43.

Also, the region C substitution part36calculates the differential value for the first energy difference Δg1of each of the intervals λg1acquired from the Δg1list42, and the differential value for each of the intervals λg1in the region C to the Δg1differential value list42-2.

After the differential values in the region A and the region C are acquired, the region B interpolation part37refers to the differential value of the second energy difference Δg2corresponding to λ1of the region A and the differential value of the second energy difference Δg2corresponding to λ2of the region C, and interpolates the differential value of the second energy difference Δg2in the region B by using the three dimensional spline or the like (step S20).

The region B interpolation part37inversely calculates the second energy difference Δg2at the intervals λ in the region B from the differential value, which is interpolated by the region B interpolation part37, of the second energy difference Δg2in the region B, and additionally stores the second energy difference Δg2, which is corresponded to the identification information of the candidate compound, at the intervals λ in the region B into the Δg2list43in the storage part50.

After the second energy difference Δg2is calculated in the entirety of the region of 0≦λ≦1 for the first candidate compound, and the region of 0≦λ≦1 is divided into the regions A, B, and C, or after the above described processes respective to the regions A, B, and C for each of the candidate compounds other than the first candidate compound, the binding strength calculation part38calculates the binding strength (step S16). The binding strength calculation part38refers to the Δg1list42and acquires the first energy difference ΔG1by calculating the total value of the first energy difference Δg1at the intervals λ. Moreover, the binding strength calculation part38refers to the Δg2list43and acquires the second energy difference ΔG2by calculating the total value of the first energy difference Δg2at the intervals λ. Furthermore, the binding strength calculation part38adds the first energy difference ΔG1and the second energy difference ΔG2to acquire the binding strength. The binding strength calculation part38additionally stores the acquired binding strength, which corresponds to the identification information of the candidate compound, into the binding strength data46in the storage part50.

Thus, it is determined whether the above described processes are completed for all candidate compounds listed in the candidate compound list41(step S17). When the processes are completed, the drug development simulator100terminates the simulation for acquiring the binding strength.

In step S17, when it is determined that the processes are not completed for all candidate compounds, the drug development simulator100goes back to the step S11to acquire the binding strength of a next candidate compound in the candidate compound list41.

As described above, in the drug development simulator100, the CPU11executes a process60including the steps S18, S19, and S20for the candidate compounds other than the first candidate compound. That is, the second energy difference Δg2, which is a heavier process workload in the simulation, is not calculated at the intervals λ for all candidate compounds other than the first candidate compound. Therefore, it is possible to greatly reduce time spent for a simulation process.

For example, if 100 candidate compounds are targets to be processed, the second energy difference Δg2is calculated at the intervals λ in the entirety of the region of 0≦λ≦1 in the step S14for the first candidate compound. However, for the other 99 candidate compounds, calculations respective to the regions A, B, and C are conducted in the process60. Accordingly, as illustrated with the dashed line9dinFIG. 5, it is possible to omit a simulation for calculating the second energy difference Δg2at the intervals λ for the region B and the region C.

In a case in which the embodiment is not applied (the step S15and the process60are not performed), the simulation for calculating the first energy difference Δg2in the step S12, which spends approximately one day for one candidate compound, is performed 100 times.

1 day×100 candidate compounds=100 days Therefore, 100 days are spent for the simulation. Also, if the simulation for calculating the second energy difference Δg2at the intervals λ in the step S14is conducted 100 times, in which approximately 15 days are spent for one candidate compound.

In a case in which the embodiment is applied and performed (including the step S15and the process60), days spent for the simulation in step S12are 100 days which is the same in days. However, the simulation in step S14is conducted only for the first candidate compounds.

15 days×1 candidate compound=15 days The process60spending approximately one day is performed for the other 99 candidate compounds.

1 day×99 candidate compounds=99 days Accordingly, 214 days (100 days+15 days+99 days) may be spent. Compared with the case in which the embodiment is not applied, it is possible to greatly reduce the number of days for the simulation in the embodiment.

FIG. 9is a diagram illustrating a comparison of calculation accuracy. InFIG. 9, the binding strength between the target protein and each of candidate compounds A, B, C, and D is represented by the relative free energy difference (kcal/mol). For each of the candidate compounds A, B, C, and D, an experimental value, a calculation value A in a case in which the embodiment is not applied, and a calculation value B in a case in which the embodiment is applied are illustrated inFIG. 9.

A correlation coefficient R between the experimental value and the calculation value A indicates 0.95. The correlation coefficient R between the experimental value and the calculation value B of applying the embodiment indicates 0.98. The calculation in the embodiment improves accuracy. In the simulation by the drug development simulator100according to the embodiment, it is possible to calculate a free energy difference between the target protein and a drug candidate compound in the solvent at higher speed. Moreover, it is possible to similarly calculate the free energy difference between other macromolecular compounds and the candidate compounds.

As described above, in the embodiment, two states (the unbound state2aand the fully bound state2b) are defined by setting the binding constant λ to be a parameter. When the free energy difference between the two states is acquired by aggregating ΔGiof interval lengths Δλi, a region (region B) where a convergence speed of the calculated value becomes slower is interpolated by using the three dimensional spline. Accordingly, it is possible to reduce days spent for the simulation, and to improve calculation accuracy.

In the embodiment, in a case in which a protein, a Deoxyribo Nucleic Acid (DNA), an Ribo Nucleic Acid (RNA), an organic compound, and the like are mutually and noncovalently bound in the solvent, it is possible to reduce computer resources for estimating the free energy difference by the simulation.