Method for calculating free energy

Abstract

To provide a method for calculating free energy for calculating a change in binding free energy and hydration free energy even about proteins containing a hydrophilic region and a hydrophobic region in a molecule or a molecular unit with high accuracy, a program to be used for the method, a recording medium, a calculation device and a system.SOLUTION: A method for calculating binding free energy between a protein and a ligand by using a calculation device includes the steps for: respectively assuming the protein, the ligand, and a protein and ligand complex as a solute, and calculating five specific parameters about each solute on the basis of a form index of each solute including a displacement volume, an exposure surface area, an integral value of average curvature of exposure surface and an integral value of Gaussian curvature of exposure surface; and for calculating the binding free energy between the protein and the ligand on the basis of structure energy of each solute and the five parameters of each solute.SELECTED DRAWING: Figure 2

Description

Technical Field

[0001] The present disclosure relates to a method for calculating free energy, a program used therefor, a recording medium, a computing device, and a system. More specifically, the present disclosure relates to a method for calculating binding free energy or hydration free energy, a program used therefor, a recording medium, a computing device, and a system.

Background Art

[0002] Conventionally, as drug discovery targets, target proteins (free proteins, membrane proteins, etc.) and ligands therefor have been widely studied. In drug discovery, it is important to predict the binding affinity between a target protein and a ligand. And it is known that the above binding affinity correlates with the change in free energy (binding free energy) before and after the binding of the target protein and the ligand.

[0003] On the other hand, it has long been regarded as a problem that the research on binding affinity and binding free energy requires enormous costs and labor. Therefore, from the viewpoint of reducing the enormous costs and labor required for drug discovery, in recent years, various methods for calculating predicted values of binding free energy using a computing device such as a computer have been studied.

[0004] Among the calculation methods using a computing device such as a computer, the "free energy perturbation method (or thermodynamic integration method) + molecular dynamics simulation" using molecular dynamics simulation is the mainstream theory that can accurately predict the binding free energy and the hydration free energy that is crucial for its calculation or the change in hydration free energy before and after ligand binding (Non-Patent Documents 1 and 2), and studies have been conducted from the viewpoints of improving it and shortening the calculation time. In addition, attempts have also been made to apply the energy representation method (Non-Patent Document 3) that combines molecular dynamics simulation and solution theory based on energy representation to binding free energy calculation.

Prior Art Documents

Non-Patent Documents

[0005]

Non-Patent Document 1

Non-Patent Document 2

Non-Patent Document 3

Summary of the Invention

[0006] However, even when applied to water-soluble proteins, not to mention when applied to proteins containing hydrophilic regions and hydrophobic regions in the molecule or molecular unit such as membrane proteins, methods such as the free energy perturbation method (or thermodynamic integration method) + molecular dynamics simulation and the energy representation method are not practical because of reasons such as high computational cost.

[0007] Therefore, one object of the present disclosure is to provide a method capable of calculating the binding free energy, the hydration free energy change, etc., which is also applicable to proteins containing hydrophilic regions and hydrophobic regions in the molecule or molecular unit.

[0008] In view of the above problems, the present inventors have intensively studied and as a result, newly found that it is possible to calculate the binding free energy, the hydration free energy change, etc. with high accuracy from five parameters calculated based on the morphological index of the solute.

[0009] According to one embodiment of the present disclosure, a method for calculating the binding free energy between a protein and a ligand using a computing device, comprising: assuming the protein, the ligand, and the protein-ligand complex of the protein and the ligand as solutes respectively, and based on the morphological indices of each solute including the excluded volume, the exposed surface area, the integrated value of the mean curvature of the exposed surface, and the integrated value of the Gaussian curvature of the exposed surface, for each solute, the following five parameters: (1) A first hydration energy (|ε1|), which is the first parameter; (2) A product of the absolute temperature and the first hydration entropy (|TS1|), which is the second parameter; (3) A van der Waals potential component (|ε 2,vdw |) in the second hydration energy, which is the third parameter; (4) An electrostatic potential component (|ε 2,ES |) in the second hydration energy, which is the fourth parameter; and (5) A product of the absolute temperature and the electrostatic potential component in the second hydration entropy (|TS 2,ES |), which is the fifth parameter; calculating the above parameters by the computing device; calculating the binding free energy between the protein and the ligand by the computing device based on the structural energy of each solute and the above five parameters of each solute; A method comprising the above is provided.

[0010] Further, according to one embodiment of the present disclosure, a method for calculating the change in hydration free energy accompanying the formation of a protein-ligand complex using a computing device, comprising: assuming the protein, the ligand, and the protein-ligand complex of the protein and the ligand as solutes respectively, and calculating, by the computing device, the morphological indices of each solute including the excluded volume, the exposed surface area, the integrated value of the mean curvature of the exposed surface, and the integrated value of the Gaussian curvature of the exposed surface based on the molecular structure parameters of each solute; Based on the above morphological indices of each solute, for each solute, the following five parameters: (1) The first parameter, which is the first hydration energy (|ε1|); (2) The second parameter, which is the product of the absolute temperature and the first hydration entropy (|TS1|); (3) The third parameter, which is the van der Waals potential component in the second hydration energy (|ε 2,vdw |); (4) The fourth parameter, which is the electrostatic potential component in the second hydration energy (|ε 2,ES |); and (5) The fifth parameter, which is the product of the absolute temperature and the electrostatic potential component in the second hydration entropy (|TS 2,ES |); Calculating the above by the computing device; Calculating the hydration free energy change by the computing device by subtracting the sum of the five parameters in the protein and the sum of the five parameters in the ligand from the sum of the five parameters in the protein-ligand complex; A method comprising the above is provided.

[0011] Also, according to an embodiment of the present disclosure, a program for causing a computing device to execute the above method is provided.

[0012] Also, according to an embodiment of the present disclosure, a computer-readable recording medium recording the above program is provided.

[0013] Also, according to an embodiment of the present disclosure, a computing device having the above program recorded in an internal storage unit is provided.

[0014] Also, according to an embodiment of the present disclosure, a system including the above computing device is provided.

[0015] According to the present disclosure, it becomes possible to calculate the binding free energy and the hydration free energy change with high accuracy for proteins containing hydrophilic and hydrophobic regions within a molecule or molecular unit.

Brief Description of the Drawings

[0016]

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

Figure 14

Figure 15

Mode for Carrying Out the Invention

[0017] [Definition] As used herein, the "protein" is not particularly limited as long as it is a protein to which the ligand described later can bind. Therefore, proteins include free proteins (e.g., water-soluble proteins), membrane proteins, etc. The protein may be one in which two or more identical or different ones bind, associate, etc. to form a unit. Further, the protein may contain one or more amino acid mutations.

[0018] As used herein, the "membrane protein" means a protein attached to a biological membrane. Membrane proteins include integral membrane proteins that penetrate the lipid bilayer or bind to the lipid bilayer by fatty acid chains, etc., and peripheral membrane proteins that bind to the hydrophilic region of the lipid bilayer or other membrane proteins by non-covalent bonds. The membrane protein may be a single-pass type or a multi-pass type.

[0019] Examples of proteins include, but are not limited to, free proteins such as albumin, globulin (e.g., α1 globulin, α2 globulin, β globulin, γ globulin (including various antibodies)), fibrinogen; G protein-coupled receptors (e.g., adenosine receptor, adrenergic receptor, angiotensin receptor, opioid receptor, orexin receptor, gastrin receptor, glucagon receptor, cholecystokinin receptor, secretin receptor, serotonin receptor (excluding 5-HT3 receptor), somatostatin receptor, dopamine receptor, histamine receptor, muscarinic acetylcholine receptor, GABAB receptors, such as P2Y receptors), ligand-gated ion channels (e.g., glycine receptors, 5-HT3 receptors, AMPA receptors, GABA A receptors, NMDA receptors, P2X receptors, etc.), voltage-gated ion channels (e.g., chloride channels, potassium channels, calcium channels, sodium channels, proton channels, etc.), and other membrane proteins; and the like. In one embodiment of the present disclosure, the protein is preferably a mammalian protein, more preferably a primate protein, and even more preferably a human protein.

[0020] In one embodiment of the present disclosure, the protein is a membrane protein (preferably a human membrane protein). In one embodiment of the present disclosure, the protein is a G protein-coupled receptor (preferably a human G protein-coupled receptor).

[0021] As used herein, "ligand" means something that can bind to the above protein to form a protein-ligand complex. The ligand is not particularly limited in terms of molecular weight, size, origin (e.g., natural origin or synthetic origin), etc. Ligands include, but are not limited to, for example, small organic molecules (low molecular weight compounds), nucleic acid-derived substances (e.g., DNA, miRNA, RNA such as aptamers), and amino acid-derived substances (e.g., antibodies or their antigen-binding fragments, interleukins (IL-2, IL-7, IL-12, IL-15, etc.), chemokines, interferons (IFN-γ, etc.), insulin, erythropoietin, transforming growth factor, lymphotoxin, adiponectin, leptin, cytokines, peptide hormones, their modified forms, their mutants, etc.). Therefore, when the ligand is amino acid-derived, the protein-ligand complex can also be a protein-protein complex.

[0022] As used herein, the term "binding free energy" refers to the amount of change in free energy when a protein and a ligand bind to form a protein-ligand complex. Therefore, the binding free energy (ΔG) is typically calculated by multiplying the free energy of the protein-ligand complex-water system (G P-L ) to obtain the free energy of the protein-water system (G P ) and the free energy of the ligand-water system (G L ) can be calculated by subtracting (ΔG = G P-L -(G P +G L )). Alternatively, the binding free energy (ΔG) is the structural energy change (ΔE c ) from the hydration free energy change (Δμ) before and after ligand binding, absolute temperature (T) and the change in structural entropy of the solute before and after ligand binding (ΔS C ) product (TΔS C ) and subtract (ΔG = ΔE C -Δμ-TΔS C In addition, when calculating the binding free energy and hydration free energy changes, the above structural entropy change (ΔS C ) is assumed to be small and can be considered to be substantially cancelled out (i.e., ΔS c can be considered to be substantially 0). In this case, the binding free energy (ΔG) is substantially the same as the above structural energy change (ΔE C ) by subtracting the above hydration free energy change (Δμ) from the above (ΔG = ΔE C -Δμ). On the other hand, the structural entropy change (ΔS C When it is not appropriate to ignore the binding free energy (ΔG), ΔE C -Δμ-TΔS C It may be calculated as:

[0023] In addition, in the method described in the Examples below, the free energy of the ligand complex in the binding pose of the crystal structure (G P-L,wild) and the free energy (G P-L,decoy ) of the ligand complex in the false binding pose, the binding pose binding free energy (ΔG wild ) of the crystal structure and the free energy (ΔG decoy ) of the ligand complex in the false binding pose, the difference between them can be calculated as "ΔG decoy - ΔG wild = G P-L,decoy - G P-L,wild ". This is because the binding pose binding free energy of the crystal structure can be described as "ΔG wild = ΔG P-L,wild - (G P,wild + G L,wild )", and the free energy of the ligand complex in the false binding pose can be described as "ΔG decoy = G P-L,decoy - (G P,decoy + G L,decoy )". When the receptor and the ligand are separated, "G P,wild = G P,decoy " and "G L,wild = G L,decoy ", so the difference is canceled, and thus the above "ΔG decoy - ΔG wild = G P-L,decoy - G P-L,wild " is considered to hold.

[0024] As used herein, the "morphological index" represents the solvation thermodynamic quantity (typically, the hydration thermodynamic quantity) of an arbitrary solute as the morphological index. The "solvation thermodynamic quantity" means the change in the thermodynamic quantity that occurs when a solute with a fixed three-dimensional structure is inserted into a solvent (e.g., water). The morphological index includes at least four indices: the excluded volume, the exposed surface area, the integral value of the mean curvature of the exposed surface, and the integral value of the Gaussian curvature of the exposed surface.

[0025] As used herein, the "excluded volume" means the volume of the space (also referred to as the "excluded space" herein) into which the centers of solvent molecules (e.g., water) cannot enter in an arbitrary solute.

[0026] As used herein, the "exposed surface area" means the surface area of the excluded volume for any solute.

[0027] As used herein, the "integrated value of the mean curvature of the exposed surface" means the integrated value of the product of the mean curvature 1 / R and the exposed surface area when the solute is regarded as a combination of spheres of various sizes and a certain sphere has a radius R (that is, the value obtained by calculating the product of the mean curvature 1 / R and the exposed surface area for all spheres and then calculating their sum).

[0028] As used herein, the "integrated value of the Gaussian curvature of the exposed surface" means that when the solute is regarded as a combination of spheres of various sizes and a certain sphere has a radius R, the integrated value of the product of the Gaussian curvature 1 / R 2 and the exposed surface area (that is, the value obtained by calculating the product of the Gaussian curvature 1 / R 2 and the exposed surface area for all spheres and then calculating their sum).

[0029] As used herein, the "constant volume condition" means that there is no volume change.

[0030] In this specification, a symbol or numerical value enclosed by "||" (for example, |ε1|) means an absolute value. Therefore, a symbol or numerical value enclosed by "||" includes both positive and negative values ((for example, |ε1|) includes "ε1" and "-ε1").

[0031] As used herein, the "first hydration energy" (also denoted as "ε1" in this specification) means the hydration energy in Process 1 when the process of solute hydration is divided into two processes (Process 1 and Process 2) described below. In this specification, the first hydration energy is also denoted as "|ε1|", "ε1". In the present disclosure, the first hydration energy (preferably, the first hydration energy under the constant volume condition (also denoted as "|ε VH,1 |", "ε VH,1 ") means the first parameter.

[0032] As used herein, the "first hydration entropy" means the hydration entropy in Process 1 when the process of solute hydration is divided into the following two processes (Process 1 and Process 2). In this specification, the first hydration entropy is also denoted as "S1". In the present disclosure, the product of the absolute temperature (T) and the first hydration entropy (also denoted as "|TS1|", "-TS1" herein) (preferably, the product of the absolute temperature (T) and the first hydration entropy under constant volume conditions (also denoted as "|TS VH,1 |", "-TS VH,1 " herein)) means the second parameter.

[0033] As used herein, the "van der Waals potential component in the second hydration energy" means the component contributing to the van der Waals force among the hydration energy (ε2) in Process 2 when the process of solute hydration is divided into the following two processes (Process 1 and Process 2). In this specification, the van der Waals potential component in the second hydration energy is also denoted as "|ε 2,vdw |", "ε 2,vdw " herein. In the present disclosure, the van der Waals potential component in the second hydration energy means the third parameter.

[0034] As used herein, the "electrostatic potential component in the second hydration energy" means the component contributing to the electrostatic force among the hydration energy in Process 2 when the process of solute hydration is divided into the following two processes (Process 1 and Process 2). In this specification, the electrostatic potential component in the second hydration energy is also denoted as "|ε 2,ES |", "ε 2,ES " herein. In the present disclosure, the electrostatic potential component in the second hydration energy means the fourth parameter.

[0035] As used herein, the "electrostatic potential component in the second hydration entropy" means the component contributing to the electrostatic force among the hydration entropy in Process 2 when the process of solute hydration is divided into the following two processes (Process 1 and Process 2). Note that in this specification, the electrostatic potential component in the second hydration entropy is also denoted as "S" 2,ES ". In the present disclosure, the product of the absolute temperature (T) and the electrostatic potential component in the second hydration entropy (denoted as "|TS|", "-TS" 2,ES " in this specification) means the fifth parameter. 2,ES ).

[0036] As used herein, "structural energy" means the sum of the torsional energy, dihedral angle energy, van der Waals interaction energy, and electrostatic potential of the bond length energy of a protein or ligand. Also, in this specification, "change in structural energy" (denoted as "ΔE" C " in this specification) means the changes in bond length energy, torsional energy, dihedral angle energy, and the changes in van der Waals interaction energy and electrostatic potential within and between solutes before and after ligand binding. The change in structural energy can typically be calculated by subtracting the structural energy of the protein (E P-L ) and the structural energy of the ligand (E P ) from the structural energy of the protein-ligand complex (E L ) (ΔE C = E P-L - (E P + E L ))).

[0037] As used herein, "change in structural entropy" (denoted as "ΔS" C(also referred to as ") means the change in the structural entropy of the solute (protein) associated with ligand binding. Also, the "structural entropy" used in this specification means the thermodynamic index of the number of microscopic structures that a protein can take. The structural entropy change can be calculated by subtracting the structural entropy of the protein (S P-L ) from the structural entropy of the protein-ligand complex (S P ) and the structural entropy of the ligand (S L ) (ΔS C = S P-L - (S P + S L ))). Also, the structural entropy can be calculated using the Boltzmann-quasi-harmonic (BQH) method (S. Hikiri, T. Yoshidome, M. Ikeguchi, Journal of Chemical Theory and Computation, 12, 5990 (2016)) or the simple method developed by Kinoshita et al. (T. Yamada, T. Hayashi, S. Hikiri, N. Kobayashi, H. Yanagawa, M. Ikeguchi, M. Katahira, T. Nagata, and M. Kinoshita, Journal of Chemical Information and Modeling, 59, 3533 (2019)).

[0038] The "molecular structure parameter" used in this specification represents the molecular structure of proteins, ligands, etc. as parameters. The molecular structure parameter can be the basis for calculating the above morphological indexes and structural energies. Examples of the molecular structure parameter include atomic coordinates (x, y, z), the diameter of the molecule, the charge of the molecule, the lattice constant, etc. In one embodiment of the present disclosure, the molecular structure parameter preferably includes at least one selected from the group consisting of atomic coordinates (x, y, z), diameter, and charge, and more preferably includes atomic coordinates (x, y, z), diameter, and charge.

[0039] As used herein, the "generalized Born energy" refers to the electrostatic component of the solvation energy of a polyatomic solute calculated by considering the solvent (water) as a continuum.

[0040] As used herein, the "change in hydration free energy" (also denoted as "Δμ" herein) refers to the amount of change in hydration free energy before and after the binding of a ligand. Therefore, the change in hydration free energy (Δμ) is typically calculated by subtracting the hydration free energy of the protein (μ P-L ) and the hydration free energy of the ligand (μ P ) from the hydration free energy of the protein-ligand complex (μ L ) (Δμ = μ P-L - (μ P + μ L ))).

[0041] [Method for calculating binding free energy] According to one embodiment of the present disclosure, a method for calculating the binding free energy between a protein and a ligand using a computing device is as follows: Assume the above protein, the above ligand, and the protein-ligand complex of the above protein and the above ligand as solutes respectively, and based on the morphological indices of each solute including the excluded volume, the exposed surface area, the integrated value of the mean curvature of the exposed surface, and the integrated value of the Gaussian curvature of the exposed surface, for each solute, the following five parameters: (1) The first hydration energy (|ε1|), the first parameter; (2) The product of the absolute temperature and the first hydration entropy (|TS1|), the second parameter; (3) The van der Waals potential component in the second hydration energy (|ε 2,vdw |), the third parameter; (4) The electrostatic potential component in the second hydration energy (|ε 2,ES |), the fourth parameter; and (5) The product of the absolute temperature and the electrostatic potential component in the second hydration entropy (|TS 2,ES |), the fifth parameter; The step of calculating with the above computing device, The step of calculating, with the above computing device, the binding free energy between the protein and the ligand based on the structural energy of each solute and the above five parameters of each solute, comprises.

[0042] In one embodiment of the present disclosure, as described above, assuming the protein, the above ligand, and the protein-ligand complex of the protein and the ligand as solutes respectively, based on the morphological indices of each solute including the excluded volume, the exposed surface area, the integrated value of the mean curvature of the exposed surface, and the integrated value of the Gaussian curvature of the exposed surface, for each solute, (1) the first parameter (|ε1|); (2) the second parameter (|TS1|); (3) the third parameter (|ε 2,vdw |); (4) the fourth parameter (|ε 2,ES |); and (5) the fifth parameter (|TS 2,ES |) are calculated with the above computing device.

[0043] Hereinafter, one embodiment of the present disclosure will be described including the theoretical background of the steps. As described above, the binding free energy (ΔG) can be calculated from the structural energy change (ΔE c ) and the hydration free energy change (Δμ). Further, the hydration free energy change (Δμ) can be calculated from the hydration free energy of each solute.

[0044] On the other hand, when calculating the binding free energy, the change in hydration free energy, etc. in a protein (e.g., a membrane protein) containing a hydrophilic region and a hydrophobic region within a molecule or a molecular unit, not only the hydration in the hydrophilic region before and after ligand binding (change in hydration energy and change in hydration free energy) but also the solvation in the hydrophobic region (nonpolar solvent) before and after ligand binding (change in nonpolar solvation energy and change in nonpolar solvation free energy) need to be considered. And since the calculation of the binding free energy, etc. considering solvation in a nonpolar solvent is extremely complicated, it has been difficult to apply conventional calculation methods to proteins containing a hydrophilic region and a hydrophobic region within a molecule or a molecular unit.

[0045] Based on our own findings, we have obtained a new finding that the structure of the intramembrane region of a protein hardly changes before and after ligand binding, and the solvation of that region in a nonpolar solvent (change in nonpolar solvation energy and change in nonpolar solvation free energy) is almost canceled out when taking the difference before and after ligand binding in the calculation of ΔG. And based on such findings, when calculating the binding free energy, the change in hydration free energy, etc. of a protein containing a hydrophilic region and a hydrophobic region within a molecule or a molecular unit, the process of solvation of the solute (protein, ligand, protein-ligand complex) in water can be divided into the following Process 1 and Process 2: Process 1: The process of generating a cavity having the same (preferably the same at the atomic level) geometric characteristics as the solute in the solvent (water); Process 2: The process of inserting the van der Waals potential between the cavity-solvent (water) and the electrostatic potential between the solute-solvent (water). We have obtained a new finding that it is possible to apply a method of dividing them (Figure 1).

[0046] Note that Process 2 can also be considered separately as "Process 2 vdw " for inserting the van der Waals potential between the cavity-solvent (water) and "Process 2 ES " for inserting the electrostatic potential between the solute-solvent (water).

[0047] When considered as described above, the hydration free energy (μ) of any solute can be expressed as the sum of the hydration free energy (μ1) in Process 1 and the hydration free energy (μ2) in Process 2 (μ = μ1 + μ2). And since Process 2 can be considered separately as Process 2 vdw and Process 2 ES the hydration free energy (μ) can also be expressed as "μ = μ1 + μ 2,vdw + μ 2,ES ".

[0048] Furthermore, the present inventors obtained the finding from their own findings that the binding occurs almost under constant pressure and constant volume. Therefore, both Process 1 and Process 2 may be considered under the constant volume condition, which is easier to handle statistically. Note that since the constant volume condition is not affected by the expansion or compression of bulk water, the physical interpretation of the hydration energy and the hydration entropy may be easy. Therefore, in one embodiment of the present disclosure, Process 1 is set as the constant volume condition. In this specification, "VH" used together with symbols and numerical values means hydration under constant volume conditions.

[0049] Here, the hydration free energy can be expressed as "μ = ε - TS" (where ε is the hydration energy, T is the absolute temperature, and S is the entropy) from the common general knowledge in the field of thermodynamics. Therefore, the hydration free energy (μ1) in Process 1 can be expressed as "μ1 = ε1 - TS1" (where ε1 is the hydration energy in Process 1 (i.e., the first hydration energy under constant volume conditions), T is the absolute temperature, and S1 is the entropy in Process 1 (i.e., the first hydration entropy)). Also, the hydration free energy (μ vdw ) in Process 2 2,vdw is "μ 2,vdw = ε 2,vdw - TS 2,vdw " (where ε 2,vdw is the van der Waals potential component in the hydration energy in Process 2 (i.e., the van der Waals potential component in the second hydration energy), T is the absolute temperature, and S 2,vdwcan be expressed as the van der Waals potential component of the entropy in Process 2. Similarly, for the hydration free energy (μ ES ) in Process 2, 2,ES it can be expressed as "μ 2,ES = ε 2,ES -TS 2,ES " (where ε 2,ES is the electrostatic potential component of the hydration energy in Process 2 (i.e., the electrostatic potential component in the second hydration energy), T is the absolute temperature, and S 2,ES is the electrostatic potential component of the entropy in Process 2).

[0050] Among the above-described components, the van der Waals potential component (S 2,vdw ) of the entropy in Process 2 is considered to be a very small value compared to the electrostatic potential component (S 2,ES ) of the hydration energy in Process 2. Therefore, the product of S 2,vdw and the absolute temperature (T), "TS 2,vdw ", may be approximated to 0. In this case, the hydration free energy (μ vdw ) in Process 2 may be approximated as "μ 2,vdw = ε 2,vdw ". 2,vdw

[0051] As described above, the hydration free energy can be divided into the hydration free energies in Process 1 and Process 2, and further, these can be divided into each component. Therefore, by calculating each divided component, it becomes possible to calculate the hydration free energy, the change in hydration free energy, the binding free energy, and the like.

[0052] Based on further investigations based on the above findings, the present inventors have found that each of the above-described components (ε1, -TS1, ε 2,vdw , ε 2,ES , -TS 2,ESIt has been found that it can be calculated by applying morphological expressions. That is, the three-dimensional structure of the solute is expressed based on morphological indices (including excluded volume, exposed surface area, integral value of the mean curvature of the exposed surface, integral value of the Gaussian curvature of the exposed surface), and it has been found that each component can be calculated using this. Based on such findings, in one embodiment of the present disclosure, five parameters (the first parameter (ε1), the second parameter (-TS1), the third parameter (ε 2,vdw ), the fourth parameter (ε 2,ES ), and the fifth parameter (-TS 2,ES )) are calculated.

[0053] According to one embodiment, in the method of the present disclosure, each of the above components can be expressed, for example, by the following linear combination formula:

Equation

[0054] For the sake of computational convenience, these may be calculated after being made dimensionless using desired numerical values or constants. For example, these equations may use those divided by k B T (where k B is the Boltzmann constant and T is the absolute temperature).

[0055] According to morphometrics, each of the above Cs is considered to be a constant that does not depend on the morphological index of the solute. Since the morphological index of the solute can be calculated by a known method (for example, it can also be calculated by AlphaMol (https: / / github.com / pkoehl / AlphaMol) (see also, for example, Patrice Koehl et al., J. Chem. Inf. Model. 2023, 63, 973 - 985)), once each C is determined, for any solute, it becomes possible to calculate each of the above components (that is, the first to fifth parameters).

[0056] Each C may be calculated, for example, by the following methods (1) to (3), by a method applying these, or by a method combining these. (1) For spherical solutes (for example, proteins) having a plurality of different sizes (diameters), calculate each of the above components (the first to fifth parameters) by the angular-dependent integral equation theory (see, for example, Hansen J-P, McDonald LR Theory of simple liquids, 3rd edn. Academic Press, London (2006)). For each of these spherical solutes, based on the values of each component (the first to fifth parameters) calculated by the angular-dependent integral equation theory and the morphological index of the solute, apply the least squares method to the equation "C1V ex +C2A+C3X+C4Y" (where each symbol is as defined above) to calculate each C. (2) For a plurality of different solutes (e.g., proteins), each of the above components (the first to fifth parameters) is calculated by the three-dimensional reference interaction site model theory (also referred to herein as the "3D-RISM theory") (for example, see F. Hirata ed., Molecular Theory of Solvation, Kluwer, Dordrecht (2003)). The charges of these solutes may be 0 or non-zero. The method of the present disclosure is advantageous in that it can be calculated even when the charge of the solute is non-zero. For each solute, based on the values of each component (the first to fifth parameters) calculated by the 3D-RISM theory and the morphological index of the solute, the least squares method is applied to the formula "C1V ex +C2A+C3X+C4Y" (where each symbol is as defined above) to calculate each C. (3) For a plurality of different solutes (e.g., proteins) having a desired charge (e.g., -2e), each of the above components (the first to fifth parameters) is calculated by the 3D-RISM theory. For each solute having these desired charges, based on the values of each component (the first to fifth parameters) calculated by 3D-RISM and the generalized Born energy, and the morphological index of the solute, the least squares method is applied to the formula "C1V ex +C2A+C3X+C4Y+C5E GB +C6" (where each symbol is as defined above) to calculate each C.

[0057] Note that each C used in the examples described later was calculated using the following method. However, since this is only an example, it should be understood that each C calculated by another method may also be used. · Constants used in the calculation of "ε1" and "-TS1" (i.e., C 1-1 、C 2-1 、C 3-1 、C 4-1 、C 1-2 、C 2-2 、C 3-2 、and C 4-2 ): Use those calculated by the method of (1) above; · "ε 2,vdw " and "-TS2,ES Constants used in the calculation of " (i.e., C 1-3 , C 2-3 , C 3-3 , C 4-3 , C 1-5 , C 2-5 , C 3-5 , and C 4-5 ): Use those calculated by the method in (2) above; · Constants used in the calculation of "ε 2,ES " (i.e., C 1-4 , C 2-4 , C 3-4 , C 4-4 , C 5-4 , C 6-4 : Use those calculated by the method in (3) above.

[0058] And, according to one embodiment, in the method of the present disclosure, based on the structural energy of each solute and the above five parameters of each solute, the step of calculating the binding free energy between the protein and the ligand is executed by the above computing device.

[0059] [Processing Flow of the Method for Calculating Binding Free Energy] Hereinafter, with reference to FIG. 2, the processing flow of the method for calculating the binding free energy in one embodiment of the present disclosure will be described in more detail.

[0060] First, for each solute of the protein, ligand, and protein-ligand complex, for example, input molecular structure parameters (e.g., atomic coordinates, charges, diameters) obtained from a Protein Data Bank file and an Amber force field file (S101).

[0061] Based on the input molecular structure parameters, for each solute, morphological indices such as the excluded volume, the solvent-accessible surface area, the integral value of the mean curvature of the solvent-accessible surface, and the integral value of the Gaussian curvature of the solvent-accessible surface are calculated by the above computing device using, for example, AlphaMol (https: / / github.com / pkoehl / AlphaMol) (S102). Note that if the calculated values of the morphological indices already exist for all or part of each solute, S101 - S102 may be omitted as necessary, and these morphological indices may be directly input.

[0062] Based on these morphological indices, for each solute, the generalized Born energy is calculated by the above computing device using, for example, the Amber program (https: / / ambermd.org) (S103). Note that if the calculated values of the generalized Born energy already exist for all or part of each solute, S103 may be omitted as necessary, and the calculated values of the generalized Born energy may be directly input.

[0063] Next, based on the morphological indices and the generalized Born energy, five parameters (the first parameter (ε1), the second parameter (-TS1), the third parameter (ε 2,vdw )), the fourth parameter (ε 2,ES ), and the fifth parameter (-TS 2,ES ))) are calculated by the above computing device using, for example, the following formula:

Equation

[0064] Next, based on the above molecular structure parameters, for each solute, for example, the structural energy may be calculated by the Amber program in the above computing device (S105). In this step, in FIG. 2, for convenience, it is described after S104, but it may be at any location as long as it is a step after the input of the molecular structure parameters. Note that if the value of the structural energy already exists for all or part of each solute, S105 may be omitted within the necessary range, and these structural energies may be directly input. Note that the structural energy is the sum of the bond length energy, the angle bending energy, and the dihedral angle energy of the solute, "E int ", and the van der Waals interaction energy of the solute, "E int,vdw ", and can also be expressed as the sum of the electrostatic potential of the solute, "E int,ES ".

[0065] Based on the above five parameters and the structural energy obtained for each solute, the binding free energy is calculated in the above computing device.

[0066] According to one embodiment, in FIG. 2, first, the free energy (G) of each solute is, for example, the following formula:

Equation

[0067] And based on the free energy of each solute, for example, the following formula:

Equation

[0068] Also, according to another embodiment, in FIG. 2, the following formula: [Number] (where ΔE C means the change in conformational energy before and after binding of the ligand, Δε1 means the result of subtracting the first parameter of the protein and the first parameter of the ligand from the first parameter of the protein-ligand complex (the difference in the first parameter), TΔS1 means the result of subtracting the second parameter of the protein and the second parameter of the ligand from the second parameter of the protein-ligand complex (the difference in the second parameter), Δε 2,vdw means the result of subtracting the third parameter of the protein and the third parameter of the ligand from the third parameter of the protein-ligand complex (the difference in the third parameter), Δε 2,ES means the result of subtracting the fourth parameter of the protein and the fourth parameter of the ligand from the fourth parameter of the protein-ligand complex (the difference in the fourth parameter), TΔS 2,ESmeans the result of subtracting the fifth parameter of the protein and the fifth parameter of the ligand from the fifth parameter of the protein-ligand complex (the difference in the fifth parameter), ΔS C means the result of subtracting the structural entropy of the protein and the structural entropy of the ligand from the structural entropy of the protein-ligand complex (change in structural entropy)). The binding free energy may be calculated by the above computing device (S107). As described above, when the ligand does not change before and after binding of the ligand (i.e., considering the same ligand), the change in structural entropy (ΔS C ) may be considered to be 0.

[0069] Alternatively, according to another embodiment, in FIG. 2, the following formula:

Number

Number

[0070] Alternatively, according to another embodiment, in FIG. 2, the following formula:

Number

[0071] [Calculation method of hydration free energy change] Also, according to an embodiment of the present disclosure, as described above, a method for calculating the change in hydration free energy associated with the formation of a protein-ligand complex using a computing device is assuming the protein, the ligand, and the protein-ligand complex of the protein and the ligand as solutes, respectively, and calculating, by the computing device, morphological indices of each solute including the excluded volume, the exposed surface area, the integrated value of the mean curvature of the exposed surface, and the integrated value of the Gaussian curvature of the exposed surface based on the molecular structure parameters of each solute; Based on the above morphological indices of each solute, the following five parameters are calculated for each solute: (1) The first hydration energy (|ε1|), which is the first parameter; (2) The product of the absolute temperature and the first hydration entropy (|TS1|), which is the second parameter; (3) The van der Waals potential component in the second hydration energy (|ε 2,vdw |), which is the third parameter; (4) The electrostatic potential component in the second hydration energy (|ε 2,ES |), which is the fourth parameter; and (5) The product of the absolute temperature and the electrostatic potential component in the second hydration entropy (|TS 2,ES |), which is the fifth parameter; Calculating the above by the above computing device, Calculating the hydration free energy change by the above computing device by subtracting the sum of the above five parameters in the protein and the sum of the above five parameters in the protein from the sum of the above five parameters in the protein-ligand complex; Comprising.

[0072] [Processing flow of the method for calculating the hydration free energy change] The method for calculating the above hydration free energy change can be implemented according to the method for calculating the above binding free energy. Hereinafter, with reference to FIG. 3, the processing flow of the method for calculating the hydration free energy change in one embodiment of the present disclosure will be described in more detail.

[0073] First, for each solute of the protein, ligand, and protein-ligand complex, for example, molecular structure parameters (e.g., atomic coordinates, charges, diameters) obtained from a Protein Data Bank file and an Amber force field file are input (S201).

[0074] Based on the input molecular structure parameters, for each solute, morphological indices such as the excluded volume, the exposed surface area, the integral value of the mean curvature of the exposed surface, and the integral value of the Gaussian curvature of the exposed surface are calculated by the above computing device, for example, by AlphaMol (https: / / github.com / pkoehl / AlphaMol) (S202). Note that if the calculated values of the morphological indices already exist for all or part of each solute, S201 to S202 may be omitted as necessary, and these morphological indices may be directly input.

[0075] Based on these morphological indices, for each solute, the generalized Born energy is calculated by the above computing device, for example, by AlphaMol (https: / / github.com / pkoehl / AlphaMol) (S203). Note that if the calculated value of the generalized Born energy already exists for all or part of each solute, S203 may be omitted as necessary, and the calculated value of the generalized Born energy may be directly input.

[0076] Next, for each solute, based on the morphological indices and the generalized Born energy, five parameters (the first parameter (ε1), the second parameter (-TS1), the third parameter (ε 2,vdw ), the fourth parameter (ε 2,ES ), and the fifth parameter (-TS 2,ES )) are calculated by the above computing device, for example, by the following formula:

Equation

[0077] Based on the above five parameters obtained for each solute, for example, the following formula:

Number

[0078] Note that the above formula for calculating the hydration free energy change is the following formula:

Number

[0079] Also, in the above formula, since the part "Δε1 - TΔS1" is the sum of the difference of the first parameter (the change amount of the first parameter) and the difference of the second parameter (the change amount of the second parameter), it can also be said to mean the change amount of the first hydration free energy (Δμ1). Also, in the above formula, since the part "Δε 2,vdw +Δε 2,ES -TΔS 2,vdw " is the sum of the difference of the third parameter (the change amount of the third parameter), the difference of the fourth parameter (the change amount of the fourth parameter), and the difference of the fifth parameter (the change amount of the fifth parameter), it can also be said to mean the change amount of the second hydration free energy (Δμ2). Therefore, the above formula can also be expressed as "hydration energy change (Δμ)=Δμ1 + Δμ2".

[0080] According to the present disclosure, the calculation of the binding free energy between a protein and a ligand and the change in hydration free energy associated with the formation of a protein-ligand complex can be automated by a computing device such as a computer. Therefore, a program for causing a computing device to execute the method of the present disclosure is provided. According to the present disclosure, a computer-readable recording medium recording the program of the present disclosure is also provided. According to the present disclosure, further provided is a computing device having recorded the program of the present disclosure in its internal recording device or a calculation system for calculating the binding free energy between a protein and a ligand or the change in hydration free energy associated with the formation of a protein-ligand complex provided with the computing device of the present disclosure.

[0081] Hereinafter, with reference to FIGS. 4 to 7, an example of a functional block of a computing device and an example of a system of the present disclosure in one embodiment of the present disclosure will be shown.

[0082] [Functional Blocks of Binding Free Energy Calculation Device] FIG. 4 shows an example of a functional block diagram of a device for calculating the binding free energy (also referred to as a "binding free energy calculation device" in this specification) in one embodiment of the present disclosure, and FIG. 5 shows an example of a system for calculating the binding free energy (also referred to as a "binding free energy calculation system" in this specification) in one embodiment of the present disclosure.

[0083] In FIG. 4, the binding free energy calculation device 1 includes at least a processing unit 10 and a storage unit 11. In the example of FIG. 4, the binding free energy calculation device 1 further includes a communication unit 12 and an input / output interface unit 13.

[0084] The processing unit 10 is configured by a computer such as a CPU, for example, and can execute the necessary processing in the binding free energy calculation device. The processing unit 10 includes at least a calculation unit 10a.

[0085] The calculation unit 10a can perform various calculation processes. The calculation unit 10a includes at least a parameter calculation unit 103a and a binding free energy calculation unit 106a. In the example of FIG. 4, the calculation unit 10a further includes a morphological index calculation unit 101a, a generalized Born energy calculation unit 102a, a structural energy calculation unit 104a, and a free energy calculation unit 105a.

[0086] The morphological index calculation unit 101a can calculate at least four morphological indices, namely, the excluded volume, the exposed surface area, the integral value of the average curvature of the exposed surface, and the integral value of the Gaussian curvature of the exposed surface, for each solute of the protein, ligand, and / or protein-ligand complex. The morphological index calculation unit 101a can calculate the above four morphological indices for each of the above solutes based on the input or stored molecular structure parameters (such as atomic coordinates, diameters, charges, etc.) in the storage unit 11.

[0087] The generalized Born energy calculation unit 102a can calculate the generalized Born energy for each solute of the protein, ligand, and / or protein-ligand complex. The generalized Born energy calculation unit 102a can calculate the generalized Born energy for each of the above solutes based on the morphological indices input, stored in the storage unit 11, or calculated by the morphological index calculation unit 101a. The calculation method of the generalized Born energy may be, for example, the method described in Amber etc. (Annu. Rev. Biophys. 2019, 58, 275-296).

[0088] The parameter calculation unit 103a calculates, for each solute of the protein, ligand, and / or protein-ligand complex, a first parameter (|ε1|, preferably ε VH,1 ), a second parameter (|TS1|, preferably -TS1), a third parameter (|ε 2,vdw |, preferably ε 2,vdw ), a fourth parameter (|ε 2,ES |, preferably ε 2,ES ) and a fifth parameter (|TS 2,ES |, preferably -TS2,ES ) can be calculated. The parameter calculation unit 103a can calculate the first to fifth parameters based on the morphological indices input, stored in the storage unit 11, or calculated by the morphological index calculation unit 101a for each of the above solutes.

[0089] The structural energy calculation unit 104a can calculate the structural energy (e.g., E int , E int,vdw , E int,ES ) for each of the solutes of the protein, ligand, and / or protein-ligand complex. The structural energy calculation unit 104a can calculate the structural energy based on the molecular structure parameters (e.g., atomic coordinates, diameter, charge, etc.) input or stored in the storage unit 11 for each of the above solutes. The method for calculating the structural energy may be, for example, the method described in Amber or the like.

[0090] The free energy calculation unit 105a can calculate the free energy for each of the solutes of the protein, ligand, and / or protein-ligand complex. The free energy calculation unit 105a can calculate the free energy based on the five parameters calculated by the parameter calculation unit 103a and the structural energy input, stored in the storage unit 11, or calculated by the structural energy calculation unit 104a for each of the above solutes.

[0091] The binding free energy calculation unit 106a can calculate the binding free energy. For each of the above solutes, the binding free energy calculation unit 106a can calculate the binding free energy based on the five parameters calculated by the parameter calculation unit 103a and the structure energy calculated by the structure energy calculation unit 104a that is input, stored in the storage unit 11, or calculated by the structure energy calculation unit 104a. Alternatively, the binding free energy calculation unit 106a may calculate the binding free energy based on the free energy of each solute calculated by the free energy calculation unit 105a. In this case, the binding free energy may be calculated by subtracting the free energy of the protein and the free energy of the ligand from the free energy of the protein-ligand complex.

[0092] The storage unit 11 can store programs for processing executed by the processing unit 10, data necessary for the processing, etc. The storage unit 11 may also store a control program such as an OS (Operating System). Further, the storage unit 11 may store information such as chemical structural formulas, amino acid sequences, gene sequences, crystal structures, and molecular structure parameters (e.g., atomic coordinates, diameters, charges, etc.) for each solute of proteins, ligands, and protein-ligand complexes. The storage unit 11 may store a constant C or a combination thereof that can be used for the calculation of the first to fifth parameters.

[0093] The communication unit 12 can mutually transmit and receive data, etc. by wireless communication and / or wired communication via the communication network NW.

[0094] The input / output interface unit 13 is connected to the built-in and / or externally connected input unit 14 and / or output unit 15, and can control the input unit 14 and / or output unit 15 (note that the input unit 14 and output unit 15 are omitted in FIG. 4). Examples of the input unit 14 include a keyboard, a mouse, a microphone, etc. Examples of the output unit 15 include a display, a monitor, a speaker, etc.

[0095] Note that FIG. 4 shows an example in which each of the above-described functional blocks is provided in one computing device (computer). In this case, one computing device having each of the above-described functional blocks may constitute the binding free energy calculation system 2 in one embodiment of the present disclosure. Further, each of the above-described functional blocks may be divided and present in two or more computing devices, and these two or more computing devices may be communicably connected to each other via a communication network NW or the like. In this case, the two or more computing devices may constitute the binding free energy calculation system 2 in one embodiment of the present disclosure (FIG. 5).

[0096] The binding free energy calculation device and / or the binding free energy calculation system in one embodiment of the present disclosure may have other functional blocks as necessary in addition to the above-described functional blocks. For example, when calculating the conformational entropy change (ΔS C ), the binding free energy calculation device and / or the binding free energy calculation system in one embodiment of the present disclosure may include a conformational entropy change calculation unit. The conformational entropy change calculation unit may calculate the conformational entropy change (ΔS P-L ) by subtracting the conformational entropy (S P ) of the protein and the conformational entropy (S L ) of the ligand from the conformational entropy (S C ) of the protein-ligand complex. The conformational entropy change calculation unit or other components may calculate the conformational entropy (S P-L ) of the protein-ligand complex, the conformational entropy (S P ) of the protein, and / or the conformational entropy (S L ) of the ligand by the above-described Boltzmann-quasi-harmonic (BQH) method or the simple method developed by Kinoshita et al. Note that when calculating the conformational entropy change (ΔS C ) in the process flow diagram in FIG. 2, the conformational entropy change (ΔS C ) may be calculated at an arbitrary position (for example, between S105 and S106).

[0097] [Functional Blocks of Hydration Free Energy Change Calculation Device] Further, FIG. 6 shows an example of a functional block diagram of a device for calculating the hydration free energy change (also referred to as the “hydration free energy change calculation device” in this specification), and FIG. 7 shows an example of a system for calculating the hydration free energy change (also referred to as the “hydration free energy change calculation system” in this specification) according to an embodiment of the present disclosure.

[0098] In FIG. 6, the hydration free energy change calculation device 3 includes at least a processing unit 30 and a storage unit 31. In the example of FIG. 6, the hydration free energy calculation device 3 further includes a communication unit 32 and an input / output interface unit 33.

[0099] The processing unit 30 is configured by a computer such as a CPU, for example, and can execute necessary processing in the hydration free energy change calculation device. The processing unit 30 includes at least a calculation unit 30a.

[0100] The calculation unit 30a can perform various calculation processes. The calculation unit 30a includes at least a parameter calculation unit 303a and a hydration free energy change calculation unit 305a. In the example of FIG. 6, the calculation unit 30a further includes a morphological index calculation unit 301a, a generalized Born energy calculation unit 302a, and a hydration free energy calculation unit 304a.

[0101] The morphological index calculation unit 301a can calculate at least four morphological indices, namely, the excluded volume, the exposed surface area, the integral value of the mean curvature of the exposed surface, and the integral value of the Gaussian curvature of the exposed surface, for each solute of a protein, a ligand, and / or a protein-ligand complex. The morphological index calculation unit 301a can calculate the above four morphological indices based on the input or stored molecular structure parameters (such as atomic coordinates, diameters, charges, etc.) in the storage unit 31 for each of the above solutes.

[0102] The generalized Born energy calculation unit 302a can calculate the generalized Born energy for each solute of a protein, a ligand, and / or a protein-ligand complex. The generalized Born energy calculation unit 302a can calculate the generalized Born energy for each of the above solutes based on the morphological indices input, stored in the storage unit 31, or calculated by the morphological index calculation unit 301a. The calculation method of the generalized Born energy may be, for example, the method described in Amber or the like.

[0103] The parameter calculation unit 303a can calculate a first parameter (|ε1|, preferably ε1), a second parameter (|TS1|, preferably -TS1), a third parameter (|ε 2,vdw |, preferably ε 2,vdw ), a fourth parameter (|ε 2,ES |, preferably ε 2,ES ) and a fifth parameter (|TS 2,ES |, preferably -TS 2,ES ) for each solute of a protein, a ligand, and / or a protein-ligand complex. The parameter calculation unit 303a can calculate the first to fifth parameters for each of the above solutes based on the morphological indices input, stored in the storage unit 31, or calculated by the morphological index calculation unit 301a.

[0104] The hydration free energy calculation unit 304a can calculate the hydration free energy for each solute of a protein, a ligand, and / or a protein-ligand complex. The hydration free energy calculation unit 304a can calculate the hydration free energy for each of the above solutes based on the five parameters calculated by the parameter calculation unit 303a. The hydration free energy of each solute can typically be calculated by the sum of the above five parameters (the first to fifth parameters) of each solute.

[0105] The hydration free energy change calculation unit 305a can calculate the hydration free energy change. The hydration free energy change calculation unit 305a can calculate the hydration free energy change based on the above five parameters calculated by the parameter calculation unit 303a. Alternatively, the hydration free energy change calculation unit 305a may calculate the hydration free energy change based on the hydration free energy of each solute calculated by the hydration free energy calculation unit 304a. In this case, the hydration free energy change may be calculated by subtracting the hydration free energy of the protein and the hydration free energy of the ligand from the hydration free energy of the protein-ligand complex.

[0106] The storage unit 31 can store programs for the processes executed by the processing unit 30, data necessary for the processes, and the like. Further, the storage unit 31 may store a control program such as an OS (Operating System). Furthermore, the storage unit 31 may store information such as chemical structural formulas, amino acid sequences, gene sequences, crystal structures, and molecular structure parameters (for example, atomic coordinates, diameters, charges, etc.) for each solute of the protein, ligand, and protein-ligand complex. The storage unit 31 may store a constant C or a combination thereof that can be used for the calculation of the first to fifth parameters.

[0107] The communication unit 32 can transmit and receive data and the like to and from each other by wireless communication and / or wired communication via the communication network NW.

[0108] The input / output interface unit 33 is connected to the built-in and / or externally connected input unit 34 and / or output unit 35, and can control the input unit 34 and / or output unit 35 (note that the input unit 34 and output unit 35 are omitted in FIG. 6). Examples of the input unit 34 include a keyboard, a mouse, a microphone, and the like. Examples of the output unit 35 include a display, a monitor, a speaker, and the like.

[0109] Note that FIG. 6 shows an example in which each of the above-described functional blocks is provided in one computing device (computer). In this case, one computing device having each of the above-described functional blocks may constitute the hydration free energy change calculation system 4 in one embodiment of the present disclosure. Further, each of the above-described functional blocks may be divided and present in two or more computing devices, and these two or more computing devices may be communicably connected to each other via a communication network NW or the like (FIG. 7). In this case, the two or more computing devices may constitute the hydration binding free energy change calculation system 4 in one embodiment of the present disclosure.

[0110] The hydration free energy change calculation device and / or the hydration free energy change calculation system in one embodiment of the present disclosure may have other functional blocks as necessary in addition to the above-described functional blocks.

[0111] The present disclosure includes the following. [1] A method for calculating the binding free energy between a protein and a ligand using a computing device, assuming the protein, the ligand, and the protein-ligand complex of the protein and the ligand as solutes, respectively, and based on the morphological indices of each solute including the excluded volume, the exposed surface area, the integral value of the mean curvature of the exposed surface, and the integral value of the Gaussian curvature of the exposed surface, the following five parameters for each solute: (1) The first hydration energy (|ε1|), the first parameter; (2) The product of the absolute temperature and the first hydration entropy (|TS1|), the second parameter; (3) The van der Waals potential component (|ε 2,vdw |) in the second hydration energy, the third parameter; (4) The electrostatic potential component (|ε 2,ES |) in the second hydration energy, the fourth parameter; and (5) The product of the absolute temperature and the electrostatic potential component in the second hydration entropy (|TS 2,ES |), the fifth parameter; The step of calculating with the above computing device, The step of calculating, with the above computing device, the binding free energy between the protein and the ligand based on the structural energy of each solute and the above five parameters of each solute, A method comprising the above. 〔2〕The method according to 〔1〕, further comprising the step of calculating, with the above computing device, the above morphological index for each solute based on each molecular structure parameter in each of the above solutes. 〔3〕The method according to 〔1〕 or 〔2〕, further comprising the step of calculating, with the above computing device, the generalized Born energy based on the above morphological index. 〔4〕The above five parameters are represented by the following mathematical formula:

Number

Equation

[10] A program for causing a computing device to execute the method according to any one of [1] to [9].

[11] A computer-readable recording medium having recorded thereon the program according to

[10] .

[12] A computing device having recorded in its internal storage the program according to

[10] .

[13] A system comprising the computing device according to

[12] .

Examples

[0112] Hereinafter, the method and the like in one embodiment of the present disclosure will be described in more detail using examples, but it is not intended to limit the scope of the present disclosure in any way. In the following examples, the above process 1 was set as a constant volume condition.

[0113] The coordinates (x, y, z) of all atoms in each solute (ligand, protein, ligand-protein complex) used in the following examples were obtained from the PDB file (RCSB Protein Data Bank (RCSB PDB) (https: / / www.rcsb.org / )). Also, the diameter and charge of each solute used in the examples were obtained from the Amber force field (University of California, San Francisco Dept. of Pharmaceutical Chemistry, https: / / ambermd.org).

[0114] Also, in the following examples, each parameter (the first parameter (ε VH,1 ), the second parameter (-TS VH,1 ), the third parameter (ε 2,vdw ), the fourth parameter (ε 2,ES ), and the fifth parameter (-TS 2,ES )) was calculated using the following formula:

Equation

[0115]

Table 1

[0116] [Example 1: Evaluation of Binding Affinity between Membrane Protein and Ligand] Using the adenosine A2A receptor, which is a membrane protein with a known three-dimensional structure, and its known ligands caffeine, XAC, and ZM241385, the binding free energy between the protein and the ligand was calculated by the method of the present disclosure, and the binding affinity of these ligands for the adenosine A2A receptor was evaluated.

[0117]

Chemical formula

[0118] 1. Calculation of the First to Fifth Parameters For each solute of the adenosine A2A receptor, each ligand, and the adenosine A2A receptor - each ligand complex, based on the obtained molecular structure parameters (all-atom coordinates, diameter, and charge), the Connolly's algorithm was combined with the TINKER program package (https: / / dasher.wustl.edu / tinker / ) to calculate the excluded volume (V ex ), the exposed surface area (A), the integral value (X) of the mean curvature of the exposed surface, and the integral value (Y) of the Gaussian curvature of the exposed surface, and the generalized Born energy (E GB ) was calculated by the computer using the Amber program. Next, based on the above four morphological indices and the above generalized Born energy of each solute, for each solute, five parameters (the first to fifth parameters) were applied to the above formula and calculated using a computing device.

[0119] 2. Calculation of Structural Energy For each solute, the structural energy (E int 、E int,vdw 、E int,ES ) was calculated using a computing device with the Amber program.

[0120] 3. Calculation of Free Energy Based on the five parameters (the first to fifth parameters), the structural energy (E int 、E int,vdw 、E int,ES ), and the generalized Born energy (E GB ) obtained above for each solute, the following formula:

Equation

[0121] 4. Calculation of Binding Free Energy Based on the free energy (G) of each solute obtained, the following formula:

Equation

[0122] 5. Calculation of Other Thermodynamic Indicators For the first parameter (ε VH,1 ) obtained for each solute, the difference (Δε VH,1 ) in the first parameter before and after ligand binding was calculated using the following formula:

Equation

[0123]

Table 2

[0124] [Example 2: Prediction of the Binding Pose of the Crystal Structure 1] Using the adenosine A2A receptor, a membrane protein with a known three-dimensional structure, and its known ligand, ZM241385, the binding free energy between the protein and the ligand was calculated by the method of the present disclosure, and it was evaluated whether it was possible to predict the binding pose (the binding pose considered to be correct) of the adenosine A2A receptor-ZM241385 complex.

[0125] 1. Creation of Binding Pose Candidates Based on the obtained molecular structure parameters of the adenosine A2A receptor and ZM241385, a plurality of binding pose candidates (i.e., a plurality of candidates for the adenosine A2A receptor-ZM241385 complex) were created using a known program (Integrated Computational Science System (MOE)). Then, using a known program (Amber), the structure of each obtained binding pose candidate was optimized. The data of the structure thus optimized includes information on the molecular structure parameters (atomic arrangement, diameter, and charge) of the structure.

[0126] 2. Calculation of Binding Free Energy For the adenosine A2A receptor, ZM241385, and the candidate of the binding pose with optimized structure (the candidate of the adenosine A2A receptor-ZM241385 complex), the binding free energy (ΔG) was calculated by a computing device in the same manner as in Example 1. This calculation was performed for each of the plurality of binding pose candidates obtained above. The time required for the calculation of the binding free energy by the method of the present disclosure was about 1 to 2 seconds for each binding pose.

[0127] 3. Calculation and Plotting of Root Mean Square Deviation For each binding pose, the root mean square deviation (RMSD, calculated by Molecular Operating Environment (MOE)) was calculated. Then, for each binding pose, the binding free energy (ΔG) was plotted on the vertical axis and the root mean square deviation (RMSD) was plotted on the horizontal axis. The results are shown in FIG. 8.

[0128] [Example 3: Prediction of the Binding Pose of the Crystal Structure 2] Using the adenosine A2A receptor, a membrane protein with a known three-dimensional structure, and its known ligand, 5'-(N-ethylcarboxamido)adenosine (NECA), the binding free energy between the protein and the ligand was calculated by the method of the present disclosure, and it was evaluated whether it was possible to predict the binding pose (the binding pose considered to be the correct answer) of the adenosine A2A receptor-NECA complex.

[0129] The same method as in Example 2 was carried out. The results are shown in Fig. 9.

[0130] [Example 4: Prediction of the binding pose of the crystal structure 3] Using the serotonin 2A receptor, a membrane protein with a known three-dimensional structure, and its known ligand, risperidone, the binding free energy between the protein and the ligand was calculated by the method of the present disclosure, and it was evaluated whether it was possible to predict the binding pose (the binding pose considered to be the correct answer) of the serotonin 2A receptor-risperidone complex.

[0131] The same method as in Example 2 was carried out. The results are shown in Fig. 10.

[0132] [Example 5: Prediction of the binding pose of the crystal structure 4] Using the serotonin 2A receptor, a membrane protein with a known three-dimensional structure, and its known ligand, zotepine, the binding free energy between the protein and the ligand was calculated by the method of the present disclosure, and it was evaluated whether it was possible to predict the binding pose (the binding pose considered to be the correct answer) of the serotonin 2A receptor-zotepine complex.

[0133] The same method as in Example 2 was carried out. The results are shown in Fig. 11.

[0134] [Example 6: Prediction of the binding pose of the crystal structure 5] In Example 2, the binding free energy (ΔG) was calculated by the following formula

Equation

[0135] [Example 7: Prediction of the binding pose of the crystal structure 6] In Example 3, except for calculating the binding free energy (ΔG) by the method of Example 6, the same method was implemented to evaluate whether it was possible to predict the binding pose (the binding pose considered to be correct) of the adenosine A2A receptor-NECA complex. The results are shown in Figure 13.

[0136] [Example 8: Prediction of the binding pose of the crystal structure 7] In Example 4, except for calculating the binding free energy (ΔG) by the method of Example 6, the same method was implemented to evaluate whether it was possible to predict the binding pose (the binding pose considered to be correct) of the serotonin 2A receptor-risperidone complex. The results are shown in Figure 14.

[0137] [Example 9: Prediction of the binding pose of the crystal structure 8] In Example 5, except for calculating the binding free energy (ΔG) by the method of Example 6, the same method was implemented to evaluate whether it was possible to predict the binding pose (the binding pose considered to be correct) of the serotonin 2A receptor-zotepine complex. The results are shown in Figure 15.

[0138] From the results of Test Example 1, the binding free energy (ΔG) of each ligand to the adenosine A2A receptor was calculated as caffeine > XAC > ZM241385 (Table 2). It is considered that the smaller the value of the binding free energy, the higher the binding affinity. Therefore, the binding affinity of these ligands to the adenosine A2A receptor was evaluated as caffeine < XAC < ZM241385, which was consistent with the publicly available information (Andrew S. Dore et al., Structure 19, 1283 - 1293, 2011, 10.1016 / j.str.2011.06.014).

[0139] Also, from each parameter calculated in Test Example 1, XAC has a smaller value than ZM241385 from the perspective of entropy (-ΔTS VH,1 difference, -ΔTS 2,ES difference) (i.e., more stable), and a larger value than ZM241385 from the perspectives of electrostatic potential and hydration energy (Δε 2,vdw difference, Δε 2,ES difference) (i.e., less stable) (Table 2). Thus, according to one embodiment of the present disclosure, not only the binding free energy (ΔG), which is considered to directly govern the binding affinity between a protein and a ligand, but also each parameter behind it (e.g., components related to the entropy change (ΔS), components related to the electrostatic potential change (Δε 2,ES ), etc.) can be calculated, and it is considered possible to quantitatively evaluate which parameter affects the binding affinity and so on. And from each obtained parameter, it is considered possible to determine the modification strategy of the ligand (e.g., increasing or decreasing the binding affinity by introducing substituents, etc.).

[0140] From the results of Example 2, there was a positive correlation that the binding free energy decreased as the root mean square deviation (RMSD) from the crystal structure (the binding pose considered to be the correct one) of each binding pose between the adenosine A2A receptor and ZM241385 decreased (Figure 8). Similarly, from the results of Example 3, there was a positive correlation that the binding free energy decreased as the RMSD from the crystal structure (the binding pose considered to be the correct one) of each binding pose between the adenosine A2A receptor and NECA decreased (Figure 9). Also, from the results of Example 4, there was a positive correlation that the binding free energy decreased as the RMSD from the crystal structure (the binding pose considered to be the correct one) of each binding pose between the serotonin 2A receptor and risperidone decreased (Figure 10). Further, from the results of Example 5, there was a positive correlation that the binding free energy decreased as the RMSD from the crystal structure (the binding pose considered to be the correct one) of each binding pose between the serotonin 2A receptor and zotepine decreased (Figure 11). Therefore, according to one embodiment of the present disclosure, it is considered possible to predict the binding pose (or the binding pose considered to be close to the correct one) considered to be the correct one for the protein-ligand complex. That is, according to one embodiment of the present disclosure, even for a protein-ligand complex with an unknown crystal structure, by combining it with a known program (for example, the integrated computational science system (MOE) or Amber), among a plurality of binding pose candidates, it is considered possible to predict the one with a small binding free energy (for example, the smallest one) as the binding pose considered to be the correct one.

[0141] From the results of Examples 6 to 9, compared with Examples 2 to 5, the binding pose considered to be the correct one could be predicted with higher accuracy. This is presumably because by considering E int,L in the calculation of ΔG, it is possible to exclude binding poses that are bent in a strange way and are forced to bind to the receptor.

[0142] According to one embodiment of the present disclosure, it becomes possible to calculate the binding free energy, the change in hydration free energy, etc. for a protein (e.g., a membrane protein) containing a hydrophilic region and a hydrophobic region within a molecule or a molecular unit. According to one embodiment of the present disclosure, it becomes possible to calculate the binding free energy, the change in hydration free energy, etc. even for a protein, a ligand, or a protein-ligand complex having various sizes, shapes, charges, etc. According to one embodiment of the present disclosure, even for a protein (e.g., a membrane protein) containing a hydrophilic region and a hydrophobic region within a molecule or a molecular unit, it becomes possible to calculate the binding free energy, the change in hydration free energy, etc. with surprisingly high precision in a short time (e.g., within several minutes, preferably within about 60 seconds, more preferably within about 30 seconds, still more preferably within about 10 seconds, even more preferably within several seconds). According to one embodiment of the present disclosure, not only the binding free energy and the change in hydration free energy, but also other components (e.g., hydration energy component, entropy component, etc.) can be decomposed, and physical factors that can promote or inhibit ligand binding and the relative strengths thereof can be quantified. Therefore, according to one embodiment of the present disclosure, it is considered possible to establish a strategy for modifying to a ligand with higher or lower binding affinity. Further, according to one embodiment of the present disclosure, by taking the statistical thermodynamics theory as the main axis instead of molecular dynamics simulation, it becomes possible to provide a method capable of calculating in a short time the binding free energy, the change in hydration free energy, etc. applicable to a protein containing a hydrophilic region and a hydrophobic region within a molecule or a molecular unit. Conventional free energy perturbation methods (or thermodynamic integration methods) + molecular dynamics simulations and energy representation methods, when applied to proteins containing hydrophilic and hydrophobic regions within the molecule or molecular unit, such as membrane proteins, result in extremely high computational costs, making it extremely difficult to calculate the binding free energy for numerous binding poses or types of ligands. Additionally, in the energy representation method, accurate calculations can only be performed when the total charge of the protein and the total charge of the ligand are zero. On the other hand, according to the method of the present disclosure, even when applied to proteins containing hydrophilic and hydrophobic regions within the molecule or molecular unit, it is advantageous in that the computational cost is low and calculations for many types of ligands (including proteins with non-zero total charges) are possible.

[0143] Taking proteins as an example, the denatured state is an ensemble of numerous different microscopic structures and is considered to have a very low degree of order. On the other hand, in the native state, it only fluctuates slightly around one defined structure and is considered to have a very high degree of order. With the binding to a ligand, the structural fluctuations of the side chains of the protein within the binding interface become smaller, the number of microscopic structures that the protein can adopt decreases, and the degree of order increases. Thus, the thermodynamic index S of the number g of microscopic structures that the protein can adopt is called the structural entropy of the protein (S = k B ln(g); k B is the Boltzmann constant). With the folding of the protein, the structural entropy significantly decreases, and with the binding to a ligand, the structural entropy of the protein is considered to decrease somewhat. The structural entropy loss (ΔSc) can be estimated by the Boltzmann - quasi - harmonic (BQH) method or the simplified method of Kinoshita et al., but there may be cases where accuracy cannot be ensured. On the other hand, according to one embodiment of the present disclosure, it is advantageous in that the calculation of ΔSc can be avoided (that is, ΔSc can be considered to be 0) without calculating ΔSc.

Explanation of Symbols

[0144] 1 Binding free energy calculation device 10 Processing Unit 10a Calculation Unit 101a Morphological Index Calculation Unit 102a Generalized Born Energy Calculation Unit 103a Parameter Calculation Unit 104a Structural Energy Calculation Unit 105a Free Energy Calculation Unit 106a Binding Free Energy Calculation Unit 11 Memory Unit 12 Communication Unit 13 Input / Output Interface Unit 2 Binding Free Energy Calculation System 3 Hydration Free Energy Change Calculation Unit 30 Processing Unit 30a Calculation Unit 301a Morphological Index Calculation Unit 302a Generalized Born Energy Calculation Unit 303a Parameter Calculation Unit 304a Hydration Free Energy Calculation Unit 305a Hydration Free Energy Change Calculation Unit 31 Memory Unit 32 Communication Unit 33 Input / Output Interface Unit 4 Hydration Free Energy Change Calculation System NW Communication Network

Claims

1. A method for calculating the binding free energy between a protein and a ligand using a computing device, comprising: assuming the protein, the ligand, and the protein-ligand complex of the protein and the ligand as solutes respectively, and based on the morphological indices of each solute including the excluded volume, the exposed surface area, the integral value of the mean curvature of the exposed surface, and the integral value of the Gaussian curvature of the exposed surface, for each solute, the following five parameters: (1) A first parameter that is a first hydration energy (|ε １ |); (2) A second parameter, which is the product of the absolute temperature and the first hydration entropy (|TS １ |); (3) The third parameter, which is the van der Waals potential component (|ε ２，ｖｄｗ |) in the second hydration energy; (4) The fourth parameter, which is the electrostatic potential component (|ε ２，ＥＳ |) in the second hydration energy; and (5) A fifth parameter that is the product of the absolute temperature and the electrostatic potential component in the second hydration entropy (|TS ２，ＥＳ |); calculating, by the computing device; calculating, by the computing device, the binding free energy between the protein and the ligand based on the structural energy of each solute and the five parameters of each solute; A method comprising the above steps.

2. The method according to claim 1, further comprising calculating, by the computing device, the morphological index for each solute based on each molecular structure parameter in each solute.

3. The method according to claim 1, further comprising calculating, by the computing device, the generalized Born energy based on the morphological index.

4. The five parameters are represented by the following mathematical formula: 【Number 1】 (wherein ε １ represents the first hydration energy, S １ represents the first hydration entropy change, T represents the absolute temperature, ε ２，ｖｄｗ represents the van der Waals potential component in the second hydration energy, ε ２，ＥＳ represents the electrostatic potential component in the second hydration energy, S ２，ＥＳ represents the electrostatic potential component in the second hydration entropy, V ｅｘ indicates the excluded volume, A represents the exposed surface area, X represents the integral value of the mean curvature of the exposed surface, Y represents the integral value of the Gaussian curvature of the exposed surface, E ＧＢ represents the generalized Born energy and each C represents a constant that may be the same or different) The method according to claim 1, calculated by the above formula.

5. The method according to claim 1, wherein the protein is a membrane protein.

6. A method for calculating the change in hydration free energy accompanying the formation of a protein-ligand complex using a computing device, comprising: assuming the protein, the ligand, and the protein-ligand complex of the protein and the ligand as solutes respectively, and calculating, by the computing device, the morphological index of each solute including the excluded volume, the exposed surface area, the integral value of the mean curvature of the exposed surface, and the integral value of the Gaussian curvature of the exposed surface based on the molecular structure parameter of each solute; based on the morphological index of each solute, for each solute, the following five parameters: (1) The first hydration energy (|ε １ |), the first parameter; The second parameter, which is the product of the absolute temperature and the first hydration entropy (|TS １ |); (3) The third parameter, which is the van der Waals potential component (|ε ２，ｖｄｗ |) in the second hydration energy; (4) The fourth parameter, which is the electrostatic potential component (|ε ２，ＥＳ |) in the second hydration energy; and (5) A fifth parameter, which is the product of the absolute temperature and the electrostatic potential component in the second hydration entropy (|TS ２，ＥＳ |); calculating, by the computing device; calculating, by the computing device, the change in hydration free energy by subtracting the sum of the five parameters in the protein and the sum of the five parameters in the ligand from the sum of the five parameters in the protein-ligand complex; A method comprising the above steps.

7. The method according to claim 6, further comprising the step of calculating a generalized Born energy by the computing device based on the morphological index.

8. The five parameters are given by the following formula: 【Number 2】 (wherein ε １ represents the first hydration energy, S １ represents the first hydration entropy, T represents the absolute temperature, ε ２，ｖｄｗ represents the van der Waals potential component in the second hydration energy, ε ２，ＥＳ represents the electrostatic potential component in the second hydration energy, S ２，ＥＳ represents the electrostatic potential component in the second hydration entropy, V ｅｘ represents the displacement volume, A represents the exposed surface area, X represents the integrated value of the mean curvature of the exposed surface, Y represents the integrated value of the Gaussian curvature of the exposed surface, E ＧＢ represents the generalized Born energy, each C represents a constant that may be the same or different) The method according to claim 6, which is calculated by the formula.

9. The method according to claim 6, wherein the protein is a membrane protein.

10. A program for causing a computing device to execute the method according to claim 1 or 6.

11. A computer-readable recording medium having recorded thereon the program according to claim 10.

12. A computing device having recorded in an internal storage unit the program according to claim 10.

13. A system comprising the computing device according to claim 12.

Images