A Systematic Analysis of the Binding Affinity ... - ACS Publications

Subscriber access provided by UQ Library

Article

A Systematic Analysis of the Binding Affinity Between the Pim-1 Kinase and Its Inhibitors Based on the MM/3D-RISM/KH Method Takeshi Hasegawa, Masatake Sugita, Takeshi Kikuchi, and Fumio Hirata J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00158 • Publication Date (Web): 11 Oct 2017 Downloaded from http://pubs.acs.org on October 17, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

A Systematic Analysis of the Binding Affinity Between the Pim-1 Kinase and Its Inhibitors Based on the MM/3D-RISM/KH Method Takeshi Hasegawa1, Masatake Sugita1, Takeshi kikuchi1 and Fumio Hirata2* 1

Ritsumeikan University, College of Life Science, Department of Bioinformatics, 1-1-1,

Noji-higashi, Kusatsu, Shiga, 525-8577, Japan 2

Toyota Physical and Chemical Research Institute, 41-1, Yokomichi, Nagakute, Aichi

480-1192, Japan *Corresponding author. Email: [email protected] Abstract A systematic study of the binding affinity for sixteen lead compounds targeting the Pim-1 kinase is reported based on the 3D-RISM/KH theory and the MD simulation. The result shows the correlation coefficient R = 0.69 between the theoretical and experimental values of the binding free energy. It demonstrates that the method is applicable to the problem of compounds screening and lead optimization, for which relative values of the free energy among the compounds have significance. We elucidate the contribution of the ligand fragments to the binding free energy. Our results

indicate

that

the

interaction

between

the

residues

and

the

triazolo[4.3-b]pyridazine scaffold as well as the phenyl ring of the ligand molecule makes significant contribution for stabilizing the complex. We further analyze, based on the 3D-RISM/KH theory, the probability distribution of a ligand fragment around the protein-ligand complex in which the substituent around the phenyl ring is removed from the ligand. The result demonstrates that the 3D-RISM/KH theory is capable of predicting the position of substitution on a ligand, which has higher affinity to a target protein.

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 29

1. Introduction The objective in the computer aided drug design (CADD) is to find a compound that tightly binds to an active site of a protein to enhance or to inhibit intrinsic functions of protein. One of the central issues of the CADD is to screen a molecule from many candidate compounds, which has highest affinity to a target protein. The other issue is optimization of the candidate compound to increase the affinity and specificity. The problem is closely related to the molecular recognition, and both of these issues can be solved by accurate prediction of the binding affinity.

Figure 1. Standard thermodynamic cycle for the Protein-Ligand binding in aqueous solution.

The binding affinity is determined thermodynamically by the three factors as shown in figure 1, the direct interaction between protein and ligand, entropy change of the solute molecules, and so called “desolvation free energy” of the both ligand and target protein, i.e., a difference of the solvation free energy between the binding state and isolated ligand and receptor. Among the three factors, a method to evaluate the direct interaction has been more or less well established, that is the molecular mechanics (MM) method. On the other hand, evaluation of the entropy change of the solute molecule immersed in solution has remained as a challenging issue. The origin of the entropy is the structural fluctuation of protein. The earlier attempt to evaluate the fluctuation based on the normal mode analysis 1,2 does not resolve the problem, because the fluctuation in vacuum is essentially different from that in solvent. Attempts to take solvent effect on the structural fluctuation into consideration have been reported. One of those is based on the principal component analysis of MD trajectories.

3-5

The method

provides the variance-covariance matrix of structural fluctuation of protein in solvent, from which one may deduce the entropy change of a protein in solution. Numbers of


Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


other methods are also introduced to increase accuracy of estimating the entropy change. 6,7

It is the desolvation free energy that has become a touchstone of ability of various methods, proposed for screening drug compounds. So called “Docking Simulations,” which rely on the bioinformatic database and/or the empirical force field, take virtually no account or highly limited account for the desolvation free energy.

8,9

So are the

methods based on the continuum model of solvent, such as PBSA and GBSA, which have been implemented in the commercial software.

10-13

Those methods rely on the

phenomenological theory treating electrostatic field created by a charge distribution. Success of the theory depends heavily on how well the boundary conditions as well as empirical parameters such as the dielectric constant are tuned. The methods for estimating the desolvation free energy based on the molecular dynamics simulation seem to be more reliable compared to those two methods described above,

14

because

they try to determine the local chemical potential of water molecules right at an active-site of target protein. Provided availability of a high performance computer, it may in principle become a powerful tool to evaluate the desolvation free energy, thereby the binding free energy. However, there remains still a big question for the method to be a practical tool for drug discovery, which usually requires screening of a drug molecule from a huge number of candidate compounds. The problem becomes more serious when realistic solution conditions, such as ionic solutions, are concerned with. The difficulty is inherent to the method that samples the configuration space of solvent molecules around and inside the protein molecule numerically. There are some methods that do not separately estimate the three factors determining the binding free energy but directly estimate the binding free energy based on the MD simulation.

15,16

Although these

methods are rigorous when the all configuration space including the solute and solvent molecules are completely explored, the numerical sampling takes a huge computational cost sometimes: for instance, when a big structural change is induced by ligand binding, the problem becomes serious. In this regard, the method based on the 3D-RISM/KH theory17,18 deserves for a special attention, because it enables to sample, or integrate, the configuration space of solvent molecules analytically, by means of the statistical mechanics. That provides a great advantage over the methods based on the molecular simulation to calculate the binding affinity of drug candidates to a target protein. Few applications of the theory to



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 29

evaluate the binding affinity of ligands to target molecules in combination with the MD simulation, including Neuraminidase, have demonstrated soundness of the method. 19-21 The most important finding from the study is that the desolvation penalty contributes to the binding free energy almost to the same amount with the direct interaction energy between protein and ligand, and they largely cancel out each other to make the net affinity marginal, or a few kcal/mol.

20

The finding suggests that any attempt

disregarding atomistic description concerning both solute and solvent for the interactions may easily produce an erroneous result, since subtle balance between the direct interaction energy and the desolvation penalty essentially determine the binding free energy. In the present paper, we report a systematic study of the binding affinity of drug candidates based on the 3D-RISM/KH theory in combination with the MD simulation, targeting the Pim-1 kinase. In order to clarify if the theory can discriminate such subtle differences among the chemical compounds, we depict a correlation of the binding affinity between experimental and theoretical results for sixteen different compounds. We further analyze the effect of each substituent of the ligand, and the probability distribution of a ligand fragment around the protein-ligand complex in which the substituent is removed from the ligand, based on the 3D-RISM/KH theory and the molecular-mechanics (MM) method. The Pim-1 kinase belongs to Pim (Proviral Integration-site MulV) family along with Pim-1, 2, and it is a Serine/Threonine Kinase. The 3D structure of a Pim-1 kinase (PDBID: 3BGQ) with an inhibitor is shown in Figure 2(a). The enzyme is expressed widely in our body due to its phosphorylation activity of substrates concerning a variety of bio-functions such as cell cycle, apoptosis, and differentiation.

22-26

In particular, the

enzyme is expressed excessively in malignant tumor such as leukemia, lymphoma, and prostatic carcinoma.

26-28

On the other hand, a mouse, the Pim-1 kinase of which is

knocked out, does not show any indication as a phenotype.

29,30

Therefore, it is

considered that such tumor may be treated by an inhibitor of the enzyme, with minor side-effect. Although there have been some reports about the compounds that tightly bound to the Pim-1 kinase,

31-34

they have not been commercialized yet as an actual

cancer drug.


Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure2. (a) A 3D structure of the Pim-1 kinase and triazolo pyridazine inhibitor termed as VX2 in the PDB. The PDB code is 3BGQ. We refer the VX2 as a ligand 1, in this study, for the sake of simplicity of

the terminology. (b) Chemical structure of the VX2. Circled part of the ligand is termed as triazolo[4.3-b]pyridazine scaffold. This is a common part of the all ligands which is applied in this study.

2. Methods 2.1. 3D-RISM/KH theory The 3D-RISM/KH theory describes the solvation structure in terms of the spatial density distribution of water molecules around a solute. Here, we provide a brief outline of the 3D-RISM/KH theory, concerning just the physics. Detailed derivation of the equation can be found elsewhere. 17,18 Let us consider the average density of solvent molecules at a position around a solute molecule. When the position is far from the solute molecule so as to be regarded as the bulk, the density will be constant which is the same as in the pure solvent. On the other hand, when it is nearby solute, the density will be “perturbed” significantly by the field due to the solute molecule, and be different from that in the bulk, depending on the strength of the perturbation. Physical meaning of the 3D-RISM equation can be readily understood in terms of a non-linear response theory as follows. Let us denote the average density of solvent atom γ at the bulk, the density nearby solute, and the density response to the perturbation as ργ, ργ(r), and ∆ργ(r), respectively. Then, statement above can be expressed as, ∆

(1)

The density response to the perturbation can be expressed in terms of a non-linear response theory as



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 29

Δ ∑ , ,

(2)

where the cα(r’) is a non-linear perturbation due to the solute molecule, and χαγ(r,r’) is a response function. (The reason why we call it “non-linear” will be clarified in a moment.) The equation can be viewed as the 3D-RISM equation

17,18,35

by identifying

cα(r) as the direct correlation function and χαγ(r,r’) as the correlation function of density fluctuation in bulk solvent, that is, , ′ 〈 ′〉,

(3)

where is the density fluctuation of atom α in the pure liquid defined by

− .

χαγ(r,r’) is the density pair correlation function in pure

solvent, that can be calculated based on the 1D-RISM theory. 35 Several approximate equations have been devised for the relation between direct

correlation function and total correlation function. The KH closure-relation, which we have employed in the present paper, reads

ℎ

− !" ℎ − # − 1 ℎ ≤ 0 − ! ℎ − ℎ > 0

(4)

In the expression, uα(r) is the direct interaction potential exerted on the solvent molecules from the solute molecule, and hα(r) is the density fluctuation in solvent at position r, normalized by the bulk density, namely,

ℎ ( Δρ (/ρ .

(5)

The three dimensional distribution function used in this study is defined from hα(r) by

+ ( h ( 1.

(6)

It is not only the direct interaction uα(r) with solute that perturbs the density of solvent at a certain position, but also that from solvent molecules at the other positions, whose density is also perturbed by the existence of the same solute. Such “indirect” perturbations are renormalized into the terms including (hα(r) - cα(r)). Such renormalization makes the perturbation highly “non-linear.” The solvation free energy ∆µ can be obtained from hα(r) and cα(r) using the following equation, 8

8

23 Δ-./01 45 6 ∑ 7 ℎ Θ−ℎ − − ℎ :,

where Θ(x) denotes the Heaviside step function. 35


(7)

Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.2. Preparation of the initial coordinate for the MD simulation We prepare the initial condition based on a PDB structure of a Pim-1 kinase (PDBID: 3BGQ) with an triazolo pyridazine inhibitor termed as VX2 in the PDB. The chemical structure of the VX2 is shown in the Figure 2(b). Other 15 ligands 34 are modeled based on the structure of VX2 and initial binding mode of the 15 ligands is determined as superimposing the ligands to VX2. We refer the VX2 to Ligand 1, in this study, and other ligands are referred as Ligand 2 to 16 for the sake of simplicity of the terminology. The ligands are listed in Table 1.

All the ligands commonly have a

triazolo[4.3-b]pyridazine scaffold. Cyclo-alkane and phenyl ring show some variety. The cyclo-alkane is substituted by cyclo-propane, cyclo-butane, and cyclo-hexane. The meta or para postion of phenyl ring are subsitituted by H, F, CF3, and OMe. The ff12SB and GAFF force field is applied to the protein and ligand structures and partial charges of the ligands are parameterized with AM1-BCC method.

36,37

For determining the

initial coordinates of the solvent molecules around and inside the protein-ligand complex, we employ the 3D-RISM/KH theory and a software called “Solutionmap.” Solutionmap is software to place the water molecules to appropriate position and orientations, in consistent with the predicted distribution function of the water molecules.

38,39

We first calculate the 3D-distribution of water molecules around the

complex based on the 3D-RISM/KH theory. We then apply “Solutionmap” to place 150 water molecules at positions where the distribution has higher peaks. The complex of Pim-1 kinase and its inhibitor, solvated by water molecules with “Solutionmap”, is immersed in a TIP3P truncated octahedral water box with a margin of 10Å from the complex with the pre-distributed water molecules. Thirteen sodium ions are also distributed as counter ions to neutralize the negative charges of residues. 2.3. Protocol for the MD simulation In this study, we minimize the system along three steps. The harmonic constraint is applied to coordinates of the Pim-1 kinase, changing the force constant in each step; 500, 50, and 20 kcal/molÅ2. Each minimizing step includes 1000 steps of a steepest descent method and 2000 steps of the conjugate gradient method. The minimized structure is heated in the NVT ensemble from 0K to 300K taking 20 ps with applying the harmonic restraint for a coordinate of a Pim-1 kinase with a force constant of 10 kcal/molÅ2. After that the system is relaxed in the NPT ensemble at 1 bar for a 1 ns. After the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

relaxation process, 200ns of a MD simulation in the NPT ensemble at 1 bar is applied. All hydrogen atoms are restrained by a SHAKE algorithm. This protocol is applied to the sixteen protein-ligand complexes prepared in the preceding procedure. All simulations are carried out using the AMBER 14 software package. 40 2.4. Free energy calculation The binding free energy consists of the potential energy change of the solute molecules,

∆Esolute, the contribution from the entropy change of the solute molecules to the free energy change, -T∆Ssolute, and the desolvation free energy, ∆Gdesolve. The desolvation free energy is estimated by the 3D-RISM/KH theory as the difference of the solvation free energies between the bound and isolated states. The 3D-RISM/KH calculation is preceded by the 1D-RISM calculation. The DRISM theory and the KH closure are employed for the calculation. We use the TIP3P model as a solvent with a slight modification for the van der Waals parameters of a hydrogen atom, which are set to ε = 0.0152 kcal/mol, σ = 1.236Å. Temperature, the number density and the dielectric constant are respectively set as, T = 300 K, ρ = 0.033 329 Å−3 and ε = 78.4. All the 3D-RISM/KH calculations are carried out with a grid spacing of 0.5Å. The box size is defined with a margin of 16Å from solute molecules. The binding energy and desolvation free energy are calculated using 2000 snap shots extracted from a trajectory of the MD simulation for 200 ns. The standard error is calculated as a standard deviation of the average values of the ∆Esolute + ∆Gdesolve in each 10 ns of the trajectory. The entropy change of the solute molecules consists of two contributions, the external entropy and the vibrational entropy, both of which are calculated using the tools implemented in AMBER. The external entropy, or the translational and rotational entropies are estimated by the classical statistical mechanical formulas, based on a rigid rotor assumption.

42

The vibrational part is estimated by means of the quasi-harmonic

method using 100000 snapshots extracted from the trajectory of the MD simulation for 200 ns.

41

We define a hypothetical state in which the solvent and the solutes (protein,

ligand and their complex) are separated in space as the standard (or reference) state of the thermodynamic cycle. The solvation free energy of each solute species is defined by the free energy change associated with the process of transferring the species from the gas phase to the infinite dilution of a solution phase, at a fixed position and orientation. So, the standard state used in the thermodynamic cycle and that employed in the


Page 8 of 29

Page 9 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


experimental data are different concerning the external entropy. The correction to the difference in the standard state is included in the external entropy calculated from the AMBER

tool,

that

takes

the

experimental

standard

state,

1

mol/l,

into

consideration. ∆Esolute and -T∆Ssolute are calculated using MMPBSA.py in the AMBER 14 software package. ∆Gdesolve is calculated using the in-house software. 2.5. Prediction of the distribution of a ligand fragment In this study, we calculate the distribution of the mixed solution of water and ligand fragments around the Pim-1 kinase and its inhibitor, whose hydrogen atom and substituent on the phenyl ring are removed, for elucidating the relationship between the position of the substituent and binding affinity. The consistency between the predicted distribution and the binding free energy of the ligand 9, 10, 13 and 14 is elucidated since those ligands have the same structure, except for the type and position of the substituent of the phenyl ring. For this objective, the modified ligand molecule is modeled by removing the hydrogen atom and the substituent from the phenyl group from ligand 9. As the ligand fragment, CF3 and OMe that are substituted around the phenyl group of the inhibitors listed in the Table 1 are selected. For the 1D-RISM calculation, the concentration of the ligand fragment is defined as 0.000005Å−3 (8.30 × 10-3 mol/L). The concentration is just for mimicking the infinite dilution condition in the 1D-RISM calculation. Detail of the number does not have much sense, except for that the concentration is very low. It has been already shown that the distribution function of the co-solvent is not sensitive to its concentration. 43



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 29

Table 1. The list of ligands and its inhibition constant. Ki is the inhibition constant.34 ∆Gbind_exp is the experimental value of the binding free energy estimated from -RTlnKi where the temperature is regarded as 300K. ∆∆Gbind_exp is the binding free energy relative to that of Ligand 1. No.

Structure

Ki

∆Gbind_exp ∆∆Gbind_exp

(nM) (kcal/mol) (kcal/mol)

No.

Structure

Ki

∆Gbind_exp ∆∆Gbind_exp

(nM) (kcal/mol) (kcal/mol)

1

11

-10.93

0.00

9

18

-10.63

0.30

2

44

-10.10

0.83

10

320

-8.92

2.01

3

94

-9.65

1.28

11

430

-8.74

2.19

4

21

-10.54

0.39

12

210

-9.17

1.76

5

160

-9.33

1.60

13

410

-8.77

2.16

6

49

-10.03

0.89

14

980

-8.25

2.68

7

100

-9.61

1.32

15

50

-10.02

0.91

8

1800

-7.89

3.04

16

54

-9.98

0.95


Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3. Results and Discussions

3.1. Comparison between the theoretical and experimental values Table 2.

Binding free energy and its components obtained from MM/3D-RISM/KH method.

∆Gbind_cal, ∆Esolute, T∆Ssolute, ∆Gdesolve denote the binding free energy, the potential energy change of the solute molecule, the entropy change of the solute molecule, and the desolvation free energy, respectively. The predicted binding free energy is shown with its average value and standard error. The binding free energy estimated from Ki, ∆Gbind_exp, is also shown.

34

∆Gbind_cal No.

∆Esolute

T∆Ssolute

∆Gdesolve

(kcal/mol)

(kcal/mol)

(kcal/mol)

(kcal/mol)

∆Gbind_exp (kcal/mol)

Ave.

Std. err.

1

-56.18

-8.06

48.02

0.26

1.33

-10.93

2

-51.81

-8.05

45.39

1.62

1.94

-10.10

3

-42.22

-12.85

34.76

5.38

1.27

-9.65

4

-50.34

-7.95

42.56

0.17

2.96

-10.54

5

-57.74

-12.83

50.09

5.17

1.69

-9.33

6

-60.98

-12.57

52.86

4.45

1.41

-10.03

7

-60.04

-15.29

52.37

7.62

1.37

-9.61

8

-44.7

-9.96

39.27

4.53

2.13

-7.89

9

-42.49

-6.54

36.64

0.69

1.99

-10.63

10

-54.48

-11.11

48.08

4.71

2.46

-8.92

11

-43.69

-9.79

39.51

5.61

1.27

-8.74

12

-44.39

-8.87

39.81

4.29

1.07

-9.17

13

-55.57

-9.6

50.61

4.63

1.67

-8.77

14

-52.04

-10.11

46.95

5.02

1.16

-8.25

15

-59.06

-13.3

49.45

3.7

1.35

-10.02

16

-53.51

-9.88

46.16

2.53

1.86

-9.98



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. Correlation between calculated and experimental values of binding free energy for all the ligands. Correlation coefficient R = 0.69. The figures are colored with different manner based on the structural feature of the ligands. (a) Ligands which include CF3 on the meta position of the phenyl ring are colored with red. Other ligands are colored with blue. (b) Ligands which have cyclo-hexane, cyclo-butane, and cyclo-propane are colored with purple, orange, and green, respectively.

The binding free energies of the sixteen compounds to the target protein, calculated based on the MM/3D-RISM/KH method, are shown in table 2 and plotted against the corresponding experimental values

34

in Figure 3 (a). A fair correlation is observed

between the theoretical and experimental results, although the theoretical results show some systematic deviation from the experimental values toward the positive side. There are several causes conceivable for the systematic deviation: insufficient sampling of the conformational space of the protein, insufficient accuracy of the 3D-RISM/KH theory for estimating the solvation free energy, inadequate solution conditions, and so on. In contrast to the method that directly estimates the binding free energy, the method based on a thermodynamic cycle requires several theoretical methodologies for estimating each component of the binding free energy; inter-atomic interactions of solute, solvation free energy, configurational entropy, and external entropy. Since each method for estimating the component of the binding free energy have own approximation, each method may systematically under- or over-estimate the thermodynamic quantity. These systematic errors seem to give rise to the unphysical positive value of the binding free energy. In the following, we discuss possible origins of the error for each component of the binding free energy. As shown in Figure S1, the sampling of the entropy change of the solute molecules


Page 12 of 29

Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


is hardly converged after 200 ns MD simulation, and it contributes systematically in positive side to the binding free energy. Extended sampling of the conformational space will definitely improve the agreement with experiments. It is known that the 3D-RISM/KH theory overestimates the solvation free energy.

44

It is also known that the difference of the calculated and experimental solvation free energies is proportional to the partial molar volume of a solute molecule. 44,45 We expect from the following argument that the overestimated part of the solvation free energy is largely canceled out in the binding process of ligand due to protein. Let’s express the ?@ V

A Systematic Analysis of the Binding Affinity ... - ACS Publications

Recommend Documents