Computational Alanine Scanning with Interaction Entropy for Protein

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Computational Alanine Scanning with Interaction Entropy for Protein−Ligand Binding Free Energies Xiao Liu,†,# Long Peng,†,# Yifan Zhou,† Youzhi Zhang,† and John Z. H. Zhang*,†,‡,§,∥ †

State Key Laboratory for Precision Spectroscopy, Shanghai Engineering Research Center of Molecular Therapeutics & New Drug Development, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062, China ‡ NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China § Department of Chemistry, New York University, New York, New York 10003, United States ∥ Collaborative Innovation Center of Extreme Optics, Shanxi University, Taiyuan, Shanxi 030006, China ABSTRACT: In protein−ligand binding, only a few residues contribute significantly to the ligand binding. Quantitative characterization of binding free energies of specific residues in protein−ligand binding is extremely useful in our understanding of drug resistance and rational drug design. In this paper, we present an alanine scanning approach combined with an efficient interaction entropy method to compute residue-specific protein−ligand binding free energies in protein−drug binding. In the current approach, the entropic components in the free energies of all residues binding to the ligand are explicitly computed from just a single trajectory MD simulation by using the interaction entropy method. In this approach the entropic contribution to binding free energy is determined from fluctuations of individual residue−ligand interaction energies contained in the MD trajectory. The calculated residue-specific binding free energies give relative values between those for ligand binding to the wild type protein and those to the mutants when specific results mutated to alanine. Computational study for the binding of two classes of drugs (first and second generation drugs) to target protein ALK and its mutant was performed. Important or hot spot residues with large contributions to the total binding energy are quantitatively characterized and the mutation effect for the loss of binding affinity for the first generation drug is explained. Finally, it is very interesting to note that the sum of those individual residue-specific binding free energies are in quite good agreement with the experimentally measured total binding free energies for this protein−ligand system.

1. INTRODUCTION The fundamental understanding of protein−ligand interaction and the quantitative characterization of binding affinity are essential in rational drug design. Accurate computation of protein−ligand binding free energy is of critical importance in theoretical design of ligands and reliable prediction of protein− ligand binding affinities. Drug-like small molecules typically prefer to bind to a few hot spots1−7 in protein−ligand and protein−protein bindings. Thus, identification of key protein residues or hot spots that provide dominant interactions to the ligands in a quantitatively measurable manner is desirable to help reveal protein−ligand binding mechanism and design mutation-resistant drugs. Once these hot spots are identified, binding interaction with adjacent areas of the protein surface can be subsequently explored to increase selectivity and improve affinity. To evaluate the quantitative contribution of a particular residue to the protein−ligand binding free energy, computational alanine scanning approach, often used in protein−protein binding,8−12 can be slightly modified to apply to protein−ligand systems. Rigorous approaches such as free energy perturbation (FEP)13−18 and thermodynamic integration methods (TI)19,20 © XXXX American Chemical Society

are computationally expensive when evaluation of a large number of protein−ligand bindings is required. As a result, the MM/GBSA (PBSA)21−25 method, which uses an implicit solvation model to compute solvation energy, is the method of choice for most practical applications in computing protein− ligand binding free energy. The MM/PBSA (GBSA) method has been successfully used to predict hot spots in protein− protein interaction in recent years.26−30 Although a large amount of computational study has focused on the calculations of total protein−ligand interaction energies with limited success, relatively few studies have focused on quantitative characterization of individual residue’s contribution to binding free energy in a consistent manner. A popular approach for computing free energy due to residue mutation is the SAAMBE approach, which is based on the MM/PBSA method.31−33 In addition, in the standard MM/PBSA (GBSA) approach to protein−ligand binding study, one often neglects the computation of entropic change due to difficulties in the calculation of standard normal mode method. Received: December 30, 2017 Published: February 6, 2018 A

DOI: 10.1021/acs.jctc.7b01295 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation In this paper, we develop a computational alanine scanning (AS) approach with the interaction entropy (IE)33 method to calculate residue-specific binding free energies in protein− ligand binding. The current method is based on a previous approach for computational alanine scanning in protein− protein system26 but specifically modified for the protein− ligand system. In the current AS-IE approach, a standard single trajectory MD simulation is performed starting from the crystal structure of the protein−ligand complex. Then, mutation is performed along the MD trajectory for all the residues of interest to generate virtual MD trajectories for mutated protein−ligand complexes. The MM/GBSA method is used to compute solvation energies along these virtual trajectories to obtain averaged enthalpy component of the binding free energy. The fluctuations of individual residue-ligand interaction energies are utilized to obtain entropic contribution to the binding free energies.26,3 This AS-IE method is applied to study bindings of three existing drugs to anaplastic lymphoma kinase (ALK) and its mutant to characterize the binding mechanisms and to quantitatively measure contributions of important residues to the total protein−ligand binding affinity. ALK has been associated with nonsmall cell lung cancer (NSCLC) and several other cancers.34−37 An ALK inhibitor crizotinib38 was approved by the U.S Food and Drug Administration (FDA) for the treatment of advanced NSCLC containing ALK rearrangements. However, many patients eventually develop resistance to crizotinib. The tumors that progress on crizotinib treatment have exhibited a range of resistance mechanisms, with ALK gene point mutations and insertions being well characterized. In 2014, two secondgeneration inhibitors (TKI), ceritinib39 and alectinib40 were developed to help overcome drug resistance. More ALK inhibitors are recently proposed to overcome drug resistance. Thus, quantitatively characterize the hot spots of ALK should be extremely useful to help understand the mechanism of drug resistance and guide the design of new inhibitors.

x→a a x ΔΔG bind = ΔG bind − ΔG bind x→a x→a = ΔΔGgas + ΔΔGsol

ΔGxbind

(1)

ΔGabind

where and are, respectively, the binding free energies of the WT (wild type) and alanine-mutated protein− ligand binding. The relative binding free energy can be partitioned into the gas-phase and solvation components defined by x→a a x ΔΔGgas = ΔGgas − ΔGgas

(2)

and x→a a x ΔΔGsol = ΔGsol − ΔGsol

(3)

The gas-phase component in eq 2 can be simplified by the approximation x→a a ΔΔGgas ≈ ΔGgas (a−L) − ΔGgxas(x−L)

(4)

where ΔGxgas(x−L) represents the gas-phase component of the free energy resulting from the interaction between residue x and the ligand L, and similarly for ΔGagas(a−L). Equation 4 is based on the standard alanine scanning approach in which the configuration of all other residues (excluding x residue) remain the same before and after the x residue is mutated to alanine. Thus, the gas-phase interaction energy between the ligand and the rest of all other residues cancels out in eq 4. 2.2. Interaction Entropy Method for Gas-Phase Component ΔΔGx→a gas . To calculate the gas-phase component of the binding free energy in eq 4, we employ a recently developed interaction entropy (IE) method.26,33,41−43 In the IE approach, the second term of free energy on the right side of eq 4 is calculated by26 x x x ΔGgas = ⟨E int ⟩ − T ΔSint x

x = ⟨E int ⟩ + KT ln⟨e βΔEint⟩

2. THEORETICAL METHODS 2.1. Alanine Scanning for Residue-Specific Binding Free Energies. In computational alanine scanning approach to protein−ligand binding, one mutates a specific residue to alanine and computes the difference in binding free energies before and after the mutation as shown schematically in Figure 1. One generally assumes that the mutated alanine does not

(5)

and similarly for the mutant a

a a ΔGgas = ⟨E int ⟩ + KT ln⟨e βΔEint⟩

(6)

Exint

In the above equations, is the interaction energy between the ligand and residue x and the exponential argument is ΔExint = Exint − ⟨Exint⟩, and there is a similar definition for the alanine mutant. The exponential average can be evaluated by discrete time averaging, x

⟨e β ΔEint⟩ =

1 N

N

∑ e βΔE

x int(ti)

(7)

i=1

where N is the number of discrete MD time steps. For the alanine mutant, the same MD trajectory obtained for the WT ALK is used except that the residue-ligand interaction energy Exint is replaced by Eaint. Finally, eq 4 gives

Figure 1. Binding free energy difference in alanine scanning mutagenesis process.

x→a x→a x→a ΔΔGgas = ΔΔEgas − T ΔΔSgas a

x

a x = ⟨E int ⟩ − ⟨E int ⟩ + KT[ln⟨e β ΔEint⟩ − ln⟨e β ΔEint⟩]

make measurable contribution to the free energy of protein− ligand binding. Thus, the difference in binding free energy before and after mutation gives a quantitative measure of the free energy contribution of the specific residue to the protein− ligand binding free energy. This difference in binding free energy for mutating residue x to a (alanine) is defined by

(8)

2.3. MM/GBSA Method for Solvation Component ΔΔGx→a sol . In our approach, the solvation energy is calculated by using the MM/GBSA method, ΔGsol = ΔGgb + ΔGnp B

(9) DOI: 10.1021/acs.jctc.7b01295 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Amber 14. A total of 200 snapshots equally distributed in the selected 1000 configurations are used to calculate the enthalpy. Interaction entropy is calculated using all of the selected configurations by time averaging of eq 7. The last 10 ns trajectory is also selected to calculate entropy contribution by normal mode method. Thus, a total of 50 configurations or snapshots are extracted from the MD trajectories for the calculation of entropy contribution. It should be mentioned that the standard MM/GBSA calculation of the total binding free energies employed the default dielectric constant of 1 for protein interior. The single trajectory approach is used to calculate the binding free energy difference which has been used by several groups.28−32 By using a single trajectory, error cancellation could overcome the insufficient sampling of the conformational space.31,32 The mutant trajectories were obtained by a single truncation of the amino acid side chain, replacing Cγ with a hydrogen atom and setting the Cβ-H direction to that of the former Cβ − Cγ based on the WT trajectory.

for both the WT and the alanine mutant. Here the polar component ΔGgb is the electrostatic solvation free energy obtained by Born (GB) model and the nonpolar component ΔGnp is estimated by using the empirical solvent accessible surface (SASA) formula available in Amber14 suite:

ΔGnp = γ SASA + β

(10)

The constants γ and β used in our calculation are the standard values of 0.00542 kcal/(mol·Å2) and 0.92 kcal/mol, respectively. Thus, eq 3 becomes x→a x→a x→a ΔΔGsol = ΔΔGgb + ΔΔGnp a x a x = ΔGgb − ΔGgb + ΔGnp − ΔGnp

(11)

2.4. Total Protein−Ligand Binding Free Energy ΔGbind. Once the residue-specific protein−ligand binding free energies ΔΔGx→a bind have been obtained by the above-described methods for all contributing residues x, one can use the following approximation to derive the total protein−ligand binding free energy by the relation x→a ΔG bind = −∑ ΔΔG bind x

3. RESULTS AND DISCUSSION In the following, we applied the IE method to alanine scanning calculations involving a total of 4 protein−ligand complexes. These protein−ligand complexes (Protein Databank ID: 2xp2, 2yfx, 4mkc, 3aox) from the PDB protein−ligand complexes are taken as the starting structures. Figure 2 shows the ALK inhibitors studied in this paper.

(12)

where the summation is over all contributing (nonalanine) residues. We call the above-described method leading to eq 12 for calculating total protein−ligand binding free energy the ASIE method, i.e., the combined method using alanine scanning and interaction entropy approach. It will be shown later in this study that eq 12 provides a much more accurate and reliable prediction of protein−ligand binding free energies than those predicted by standard MM/PBSA approach. 2.5. Molecular Dynamics Simulations. MD simulations starting from the crystal structures of protein−ligand complexes were performed using pmemd.cuda in Amber14 with the ff14SB force field. Simulations are performed in TIP3P44 explicit water with a periodic rectangular box. The distance from the surface of the box to the closet atom of the solute was set to 10 Å. The particle mesh Ewald (PME)45 is used to treat the long-range electrostatic interactions. The nonbonded interactions are rectangular with 8 Å cutoff. Periodic boundary condition (PBC) is imposed on the system during the calculation of nonbonded interactions. The time step is set at 2 fs, and SHAKE46 is used to constrain the bonds involving a hydrogen atom. A Langevin thermostat with collision frequency 2.0 is applied to control the temperature. Two minimization steps were carried out to optimize the complex structure. In the first step, only the solvent molecules and hydrogen atoms were optimized using the steepest descent algorithm followed by the conjugate gradient method. In the second step, the entire system was energy-minimized until convergence was reached. The system is then slowly heated to 300 K, followed by a 40 ns equilibration of the whole system in an NPT ensemble. A total of 4000 snapshots are saved for further analysis for the four protein−ligand complexes. The trajectory of the last 10 ns is used to calculate enthalpy and entropy contribution to the binding energies. Thus, a total of 1000 configurations or snapshots are extracted from the MD trajectories for the calculation of residue-specific interaction. The perl script mm_pbsa.pl script integrated into AMBER 14 package was used to calculate the enthalpy contribution to binding free energy. The solvation energies are calculated using the GB model47 with “igb” set to 2 in MM/GBSA protocol in

Figure 2. 2D structures of the first generation (crizotinib) and second generation (alectinib and ceritinib) inhibitors of ALK.

3.1. Protein’s Internal Dielectric Constant. The proper use of protein’s internal dielectric constant is important and its optimal value is still in debate. In the application of MM/GBSA method, the low dielectric constant of 1 or 2 is used in most applications study although larger values are applied as well in some applications. It is easy to understand that no common value can be appropriate for all systems and in all methods. In C


Article


Table 1. Residue-Specific Binding Free Energies of Top Residues for Binding of the First Generation Inhibitor Crizotinib to the WT ALK (pdbid 2xp2) from Alanine Scanning Calculationsa 2xp2

ΔΔEelec

ΔΔEvdw

ΔΔE

IE

ΔΔGx→a bind

ΔGNM bind

ΔGexp

L1256A L1122A M1199A L1196A L1198A D1203A V1130A V1180A V1270A K1150A E1197A N1254A L1204A total(−)b

−0.10 −0.41 0.55 −0.27 0.04 0.33 −0.11 0.10 0.20 −0.40 0.31 −0.01 0.03 −0.26

3.71 3.48 0.87 1.43 1.04 0.68 1.26 0.66 0.34 1.20 0.09 0.19 0.02 −14.97

3.61 3.07 1.42 1.16 1.08 1.01 1.15 0.76 0.54 0.80 0.40 0.18 0.05 −15.08

−0.86 −0.61 −0.05 −0.14 −0.11 −0.29 −0.43 −0.05 0.00 −0.27 0.00 −0.04 0.00 2.85

2.75 2.46 1.37 1.02 0.97 0.72 0.72 0.71 0.54 0.53 0.40 0.14 0.05 −12.38b

−22.53

−12.50c

⟨ΔΔEgas elec⟩

All values are in kcal/mol. ΔΔEelec includes both the gas-phase electrostatic energy and the solvation energy ΔΔEelec = + ΔΔGgb and b NM ΔΔEvdw = ⟨ΔΔEgas vdw⟩ + ΔΔGnp. The normal mode calculation of ΔGbind is described in the paper. The total binding energy is obtained by summing over all columns and multiplied by a minus sign. cThe experimental value is obtained by using the relation ΔGexp = −RT ln Ki at T = 298 K and Ki is from ref 38. a

Figure 3. Co-crystal structures of the first generation inhibitor crizotinib with (a) WT ALK (pdbid 2xp2) and (b) L1196 M mutant (pdbid 2yfx).

the study of protein−protein binding by computational alanine scanning, it was found that more accurate result is obtained by using different protein dielectric constants for different types of residues.26,28,29,31,32 In the present method, we adopt dielectric constants found in alanine scanning study for protein−protein binding with values of 1, 3, 6, 7, respectively, for nonpolar, polar, basic, and acidic residues. 3.2. First Generation Inhibitor Crizotinib. Binding of Crizotinib to WT ALK. In 2011, crizotinib became the first ALK inhibitor approved by the U.S. Food and Drug Administration (FDA) as a first-line treatment for ALK-positive lung cancer patients. Starting from the crizotinb−ALK complex structure (pdbid 2xp2), a 40 ns MD simulation was run and the last 10 ns trajectory was used to perform alanine scanning calculations. Table 1 lists top protein residues with quantitative and specific contributions to the protein−ligand binding free energy. The residue-specific binding free energies are decomposed into electrostatic (account for solvation) and van der Waals energies as well as entropic change. Table 1 tells a lot about the nature of crizotinb−ALK binding in a quantitative manner. First, the dominant contribution to the binding free energy comes from L1256 (∼3.0 kcal/mol) and L1122 (∼2.5 kcal/mol) in ALK, with the next 2−3 residues contributing more than 1 kcal/mol. each. Second, crizotinb−ALK binding is exclusively dominated by van der Waals interaction as shown in Table 1. In particular,

these van der Waals interactions all come from ligand interactions with Leu residues (4 Leu and 1 Met). The above quantitative analysis of residue-specific binding free energies can be easily explained from the crystal structure. As shown in Figure 3a, Leu1256 forms a strong π−alkyl interaction with the R3 pyridine ring of crizotinib with a distance of 4.4 Å between them. In addition, it also forms a π− alkyl interaction with R2 ring of crizotinib with a distance of 4.4 Å between them. This explains why L1256 gives the largest contribution to the crizotinib−ALK binding free energy and is therefore the most important residue. Figure 3a also shows that Leu1122 forms a strong π−alkyl interaction with the R1 ring of crizotinib with a distance of 4.3 Å between them and is the second highest ranking residue in crizotinib−ALK binding. Besides providing residue-specific binding free energies in this approach, it is important to note in Table 1 that the total sum of these residue-specific binding free energies using eq 12 is in very good agreement with the measured experimental binding free energy. For comparison, the total binding free energy from the standard MM/GBSA calculation with the default dielectric constant equal to one is about10 kcal/mol larger than the experimental value. Thus, the current AS-IE approach also provides a reliable method for obtaining the total protein−ligand binding free energy in addition to residuespecific binding free energies. D


Article

Journal of Chemical Theory and Computation Table 2. Same as Table 1 except for Binding of Crizotinib to the L1196M Mutant of ALK (pdbid 2fyx) 2yfx

ΔΔEelec

ΔΔEvdw

ΔΔE

IE

ΔΔGx→a bind

ΔGNM bind

ΔGexp

L1256A L1122A M1199A L1198A V1130A V1180A M1196A D1203A E1197A D1276A K1150A E1200A N1284A N1287A L1204A C1255A total(−)a

−0.16 −0.36 0.49 0.02 −0.08 0.12 0.17 0.34 0.32 0.20 −0.35 0.37 0.00 0.03 0.03 0.00 −1.14

3.72 3.17 0.96 0.95 1.31 0.75 1.84 0.55 0.16 0.31 0.95 0.02 0.18 0.05 0.01 0.02 −14.95

3.56 2.81 1.45 0.97 1.23 0.87 2.01 0.89 0.48 0.51 0.60 0.39 0.18 0.08 0.04 0.02 −16.09

−0.87 −0.97 −0.21 −0.11 −0.38 −0.06 −1.23 −0.23 0.03 −0.02 −0.12 −0.16 −0.03 −0.01 0.00 0.00 4.37

2.69 1.84 1.24 0.86 0.85 0.81 0.78 0.66 0.51 0.49 0.48 0.23 0.15 0.07 0.04 0.02 −11.72a

−21.95

−11.04b

a

b

The total binding energy is obtained by summing over all columns and multiplied by a minus sign. The experimental value is obtained by using the relation ΔGexp = −RT ln Ki at T = 298 K and Ki is from ref 38.

Figure 4. Main contributing residues from alanine scanning calculation for the first generation inhibitor crizotinib with WT ALK (pdbid 2xp2) and its L1196 M mutant (pdbid 2yfx).

Binding of Crizotinib to L1196 M Mutant ALK. Although crizotinib initially demonstrates robust efficacy in treating ALK positive tumors, patients eventually develop resistance. Table 2 lists the main residues from AS-IE calculation for binding of crizotinib to the L1196 M mutant of ALK. It is interesting to note that the dominant effect of L1196 M mutation on binding is the reduced binding interaction between L1122 and the R1 ring of the ligand, not due to the reduction of interaction between M1196 and the ligand! This is due to the fact that L1122 is a hot residue with a large contribution (2.46 kcal/ mol) to the binding free energy as shown in Table 1. This L1122 happens to be spatially very close to the mutating 1196 residue (cf. Figure 3) and thus its interaction with the ligand is more affected than other residues. Figure 3 shows that the distance between L1122 and the R1 ring is increased from 4.3 Å in the WT ALK to 4.7 Å in L1196 M mutant. This explains the reduction in binding free energy from AS-IE calculation. From Table 1 for L1122A in wild type 2xp2 and mutant 2yfx, the relative enthalpy change is 0.26 kcal/mol, whereas the relative

entropic change is 0.36 kcal/mol. Explanations of the mutation of this effect on crizotinib binding can also be found in ref 39. Figure 4 shows the difference in individual residues’ free energies for binding of the first generation drug crizotinib to the wild type ALK and L1196 M mutant ALK. As shown in Figure 4, the main effect of L1196 M mutation is to weaken the binding interaction between L1122 and crizotinib. 3.3. Second Generation Inhibitors FDA Approved (Ceritinib and Alectinib). Binding of Ceritinib (4mkc) to ALK. Table 3 shows the alanine scanning result of the second generation inhibitor ceritinib, the top five residues from the alanine scanning calculation are all higher than 1.0 kcal/mol. As shown in Table 3, the top two residues are the same as in binding of the first generation inhibitor crizotinib. The two top ranking residues Leu1122 (over 4.0 kcal/mol) and Leu1256 (∼3.1 kcal/mol) are hydrophobic ones with more contributions coming from the van der Waals interaction than that of the electrostatic energy. The above quantitative analysis of residuespecific binding free energies can be easily explained from the crystal structure. As shown in Figure 5a, Leu1122 forms a E


Article


Table 3. Same as Table 1 except for Binding of the Second Generation Inhibitor Ceritinib to the WT ALK (pdbid 4mkc) 4mkc

ΔΔEelec

ΔΔEvdw

ΔΔE

IE

ΔΔGx→a bind

ΔGNM bind

ΔGexp

L1122A L1256A K1150A V1130A L1198A D1203A M1199A E1132A D1270A L1196A E1200A V1180A E1197A S1206A N1254A H1124A total(−)a

−0.29 −0.22 −0.2 −0.03 0.09 0.45 0.18 0.31 0.1 0 0.43 0.04 0.25 0.02 −0.03 −0.01 −1.09

4.71 3.69 2.42 2.12 1.61 1.06 0.97 0.75 0.75 0.72 0.08 0.45 0.1 0.3 0.18 0.14 −20.05

4.42 3.47 2.22 2.09 1.7 1.51 1.15 1.06 0.85 0.72 0.51 0.49 0.35 0.32 0.15 0.13 −21.14

0.04 −0.41 −0.47 −0.2 −0.5 −0.75 −0.31 −0.54 −0.21 −1.03 −0.25 −0.04 0.02 −0.27 −0.04 −0.05 5.01

4.46 3.06 1.75 1.89 1.2 0.76 0.84 0.52 0.64 −0.31 0.26 0.45 0.37 0.05 0.11 0.08 −16.13a

−21.79

−14.04b

a

b

The total binding energy is obtained by summing over all columns and multiplied by a minus sign. The experimental value is obtained by by using 4mkc 2xp2 the relation RT ln(IC502xp2)/RT ln(IC504mkc) = ΔG2xp2 exp /ΔGexp at T = 298 K. The value of IC50 for 4mkc is from ref 39 and ΔGexp = −12.50 kcal/ mol is obtained from ref 38.

Figure 5. Co-crystal structures of the second generation inhibitors (a) ceritinib (pdbid 4mkc) and (b) alectinib (pdbid 3aox) with WT ALK.

Table 4. Same as Table 1 except for Binding of the Second Generation Inhibitor Alectinib (3aox) to the WT ALK 3aox

ΔΔEelec

ΔΔEvdw

ΔΔE

IE

ΔΔGx→a bind

ΔGNM bind

ΔGexp

L1122A L1256A V1130A L1198A R1120A M1199A D1203A E1210A E1132A L1196A K1150A V1180A E1197A I1171A D1270A E1167A total(−)a

0.03 −0.47 −0.17 0.13 −0.69 −0.05 0.38 0.44 0.37 −0.32 −0.28 −0.08 0.27 −0.1 0.19 0.17 0.18

3.27 3.13 2.06 1.21 1.93 1.1 0.58 0.43 0.49 0.98 0.93 0.66 0.14 0.42 0.09 0.09 −17.51

3.30 2.66 1.89 1.34 1.24 1.05 0.96 0.87 0.86 0.66 0.65 0.58 0.41 0.32 0.28 0.26 −17.33

−0.30 −0.36 −0.44 −0.25 −0.47 −0.2 −0.28 −0.74 −0.26 −0.55 −0.29 −0.37 −0.02 −0.04 −0.05 −0.01 4.63

3.00 2.30 1.45 1.09 0.77 0.85 0.68 0.13 0.6 0.11 0.36 0.21 0.39 0.28 0.23 0.25 −12.70a

−14.70

−12.4b

a

b

The total binding energy is obtained by summing over all columns and multiplied by a minus sign. The experimental value is obtained by using the relation ΔGexp = −RT ln Ki at T = 298 K and Ki is from ref 40.

F


Article


Figure 6. Main contributing residues from alanine scanning calculation for the second generation inhibitors ceritinib and alectinib with WT ALK.

strong π−alkyl interaction with the R3 pyridine ring of ceritinib with a distance of 4.5 Å between them. This explains why L1122 gives the largest contribution to the ceritinib−ALK binding free energy and is therefore the most important residue. Figure 5a also shows that Leu1256 forms a strong π− alkyl interaction with the R2 ring of ceritinib with a distance of 4.3 Å between them and is the second highest ranking residue in ceritinb−ALK binding. And from Figure 5a we can also learn that there is a strong π−alkyl hydrophobic interaction between the V1130 and the R1 ring of the ceritinib. The distance between the two hydrophobic centers is only 3.8 Å, which can explain the difference between these two generation inhibitors. Binding of Alectinib (3aox) to ALK. Table 4 shows that top four ranking residues binding to ALK are Leu1122 (∼3 kcal/ mol), Leu1256 (∼2.3 kcal/mol), Val1130 (∼1.45 kcal/mol), and Leu1198 (∼1.09 kcal/mol) whose binding energies are all larger than 1.0 kcal/mol. As shown in Figure 5b, these residues form hydrophobic interactions with the ligand. We can see from Figure 5b that the side chains of Leu1256 and Val1130 all have hydrophobic π−alkyl interaction with the R1 ring of the ligand and the distancse between them are 4.4 and 4.5 Å, respectively. These two residues are located on the opposite side of the R1 ring and are thus favorable to the binding affinity. There is also a π−alkyl interaction between the side chain of Leu1122 and the R3 ring of the ligand with a distance of 4 Å between them. This explains why Leu1122 (∼3 kcal/mol) gives the largest contribution to the alectinib−ALK binding free energy and is therefore the most important residue. Figure 6 shows the difference in individual residues’ free energies for binding of the second generation inhibitors ceritinib and alectinib to the wild type ALK. As shown, the binding of ceritinib to residues L1122, L1256, V1130, and K1150 of ALK are stronger than that of alectinib to these residues. 3.4. Total Binding Free Energies. The last lines of Table 1−4 give the sum of the hot spot residues’ total contribution to binding free energy. It is encouraging to see that such total binding energies obtained from eq 12 are generally in excellent agreement with the experimentally measure binding affinity. For example, −12.38 vs −12.50 (Table 1), −11.72 vs −11.04 (Table 2), −16.13 vs −14.06 (Table 3), and −12.70 vs −12.40 (Table 4) All units are kcal/mol. Except for Table 3, the computed total binding free energies for the remaining three

systems are generally in excellent agreement with the experimental data. For comparison, standard MM/GBSA calculations with normal mode method for entropy give absolute total binding free energies that are about 10 kcal/ mol larger than the corresponding experimental values. One should note, however, that the total binding free energy in eq 12 should in theory be corrected by a correction term which corresponds to the conformational energy of the reactants, i.e. corr ΔG bind = ΔG bind + ΔGcorr

(13)

where ΔGcorr is the free energy difference between the complex structures of protein and ligand and the “free” structure of the Lig protein (apo) and the ligand or ΔGcorr = ΔGpro corr + ΔGcorr. It should be noted that this correction energy in eq 13 is always positive because conformational energy upon ligand binding always raises the internal free energy of individual binding partners. However, calculation of this ΔGcorr may not be easier and thus is neglected in our study and our result is equivalent to the rigid body approach. If the protein structure is largely unchanged before and after the ligand binding and the ligand is relatively rigid, this conformational energy difference should be small and could be neglected. A recent work discussed contributions from conformational and rotational entropy in protein−ligand binding.48

4. CONCLUSIONS The present alanine scanning-interaction entropy method allows us to obtain quantitative contributions of individual residues to the total protein−ligand binding free energy. In our computational study for the three existing drugs of ALK complexes, residues Leu1256 and Leu1122 are identified to be the two most important residues that contribute to the binding free energy. The effect of L1196 M mutation on binding of the first generation drug crizotinib to ALK was explained. In addition, the total binding free energies of these three drugs to ALK and its mutant obtained by this AS-IE method are in very good agreement with the measured experimental values. The present approach should be useful to help quantitatively understand the mutational effect on drug binding and design new antiresistant drugs. G


Article


■

(15) Kita, Y.; Arakawa, T.; Lin, T. Y.; Timasheff, S. N. Contribution of the Surface Free Energy Perturbation to Protein-Solvent Interactions. Biochemistry 1994, 33, 15178−89. (16) Rao, B. G.; Singh, U. C. A free energy perturbation study of solvation in methanol and dimethyl sulfoxide. J. Am. Chem. Soc. 1990, 112, 3803−3811. (17) Kollman, P. Free energy calculations: Applications to chemical and biochemical phenomena. Chem. Rev. 1993, 93, 2395−2417. (18) Jorgensen, W. L.; Thomas, L. L. Perspective on Free-Energy Perturbation Calculations for Chemical Equilibria. J. Chem. Theory Comput. 2008, 4, 869−876. (19) Beveridge, D. L. Free Energy Via Molecular Simulation: Applications To Chemical And Biomolecular Systems. Annu. Rev. Biophys. Biomol. Struct. 1989, 18, 431−492. (20) Zacharias, M.; Straatsma, T. P.; McCammon, J. A. Separationshifted scaling, a new scaling method for Lennard-Jones interactions in thermodynamic integration. J. Chem. Phys. 1994, 100, 9025−9031. (21) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 2000, 33, 889−897. (22) Massova, I.; Kollman, P. A. Combined molecular mechanical and continuum solvent approach (MM-PBSA/GBSA) to predict ligand binding. Perspect. Drug Discovery Des. 2000, 18, 113−135. (23) Kuhn, B.; Gerber, P.; Tanja Schulzgasch, A.; Stahl, M. Validation and Use of the MM-PBSA Approach for Drug Discovery. J. Med. Chem. 2005, 48, 4040−4048. (24) Nicholls, A.; Honig, B. A rapid finite difference algorithm, utilizing successive over-relaxation to solve the Poisson−Boltzmann equation. J. Comput. Chem. 1991, 12, 435. (25) Luo, R.; David, L.; Gilson, M. K. Accelerated PoissonBoltzmann calculations for static and dynamic systems. J. Comput. Chem. 2002, 23, 1244−1253. (26) Yan, Y.; Yang, M.; Ji, C. G.; Zhang, J. Z. H. Interaction Entropy for Computational Alanine Scanning. J. Chem. Inf. Model. 2017, 57, 1112−1122. (27) Huo, S.; Massova, I.; Kollman, P. A.; Huo, S.; Massova, I.; Kollman, P. A. Computational alanine scanning of the 1:1 human growth hormone-receptor complex. J. Comput. Chem. 23, 15−27. J. Comput. Chem. 2002, 23, 15−27. (28) Martins, S. A.; Perez, M. A.; Moreira, I. S.; Sousa, S. F.; Ramos, M. J.; Fernandes, P. A. Computational Alanine Scanning Mutagenesis: MM-PBSA vs TI. J. Chem. Theory Comput. 2013, 9, 1311−1319. (29) Moreira, I. S.; Fernandes, P. A.; Ramos, M. J. Computational alanine scanning mutagenesisAn improved methodological approach. J. Comput. Chem. 2007, 28, 644−654. (30) Massova, I. M.; Kollman, P. A. Computational Alanine Scanning To Probe Protein−Protein Interactions: A Novel Approach To Evaluate Binding Free Energies. J. Am. Chem. Soc. 1999, 121, 8133− 8143. (31) Petukh, M.; Dai, L.; Alexov, E. SAAMBE: Webserver to Predict the Charge of Binding Free Energy Caused by Amino Acids Mutations. Int. J. Mol. Sci. 2016, 17, 547. (32) Petukh, M.; Li, M.; Alexov, E. Predicting Binding Free Energy Change Caused by Point Mutations with Knowledge-Modified MM/ PBSA Method. PLoS Comput. Biol. 2015, 11, e1004276. (33) Duan, L.; Liu, X.; Zhang, J. Z. H. Interaction Entropy: A New Paradigm for Highly Efficient and Reliable Computation of Protein− Ligand Binding Free Energy. J. Am. Chem. Soc. 2016, 138, 5722−5728. (34) Chiarle, R.; Voena, C.; Ambrogio, C.; Piva, R.; Inghirami, G. The anaplastic lymphoma kinase in the pathogenesis of cancer. Nat. Rev. Cancer 2008, 8, 11−23. (35) Webb, T. R.; Slavish, J.; George, R. E.; Look, A. T.; Xue, L.; Jiang, Q.; Cui, X.; Rentrop, W. B.; Morris, S. W. Anaplastic lymphoma kinase: role in cancer pathogenesis and small-molecule inhibitor development for therapy. Expert Rev. Anticancer Ther. 2009, 9, 331− 356.

AUTHOR INFORMATION

Corresponding Author

*J. Z. H. Zhang. E-mail: [email protected]. ORCID

John Z. H. Zhang: 0000-0003-4612-1863 Author Contributions #

Xiao Liu and Long Peng contributed equally to this work.

Funding

This work was supported by National Key R&D Program of China (Grant no. 2016YFA0501700), National Natural Science Foundation of China (Grant nos. 21433004, 91753103), Shanghai Putuo District (Grant 2014-A-02), Innovation Program of Shanghai Municipal Education Commission (201701070005E00020), and NYU Global Seed Grant. We thank the Supercomputer Center of East China Normal University for providing us computer time. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Cheung, Luthur S.-L.; Kanwar, M.; Ostermeier, M.; Konstantopoulos, K. A Hot-Spot Motif Characterizes the Interface between a Designed Ankyrin-Repeat Protein and Its Target Ligand. Biophys. J. 2012, 102, 407−416. (2) Barillari, C.; Marcou, G.; Rognan, D. Hot-Spots-Guided Receptor-Based Pharmacophores (HS-Pharm): A Knowledge-Based Approach to Identify Ligand-Anchoring Atoms in Protein Cavities and Prioritize Structure-Based Pharmacophores. J. Chem. Inf. Model. 2008, 48, 1396−1410. (3) Gohlke, H.; Hendlich, M.; Klebe, G. Predicting binding modes, binding affinities and ’hot spots’ for protein-ligand complexes using a knowledge-based scoring function. Perspect. Drug Discovery Des. 2000, 20, 115−144. (4) Burgoyne, N. J.; Jackson, R. M. Predicting protein interaction sites: binding hot-spots in protein-protein and protein-ligand interfaces. Bioinformatics 2006, 22, 1335−1342. (5) Chen, J.; Ma, X.; Yuan, Y.; Pei, J.; Lai, L. Protein-protein interface analysis and hot spots identification for chemical ligand design. Curr. Pharm. Des. 2014, 20, 1192−1200. (6) Bauman, J. D.; Harrison, J. J.; Arnold, E. Rapid experimental SAD phasing and hot-spot identification with halogenated fragments. IUCrJ 2016, 3, 51−60. (7) Kozakov, D.; Grove, L. E.; Hall, D. R.; Bohnuud, T.; Mottarella, S. E.; Luo, L.; Xia, B.; Beglov, D.; Vajda, S. The FTMap family of web servers for determining and characterizing ligand-binding hot spots of proteins. Nat. Protoc. 2015, 10, 733−755. (8) Keskin, O.; Gursoy, A.; Ma, B.; Nussinov, R. Principles of protein-protein interactions: what are the preferred ways for proteins to interact? Chem. Rev. 2008, 108, 1225−1244. (9) Muegge, I.; Schweins, T.; Warshel, A. Electrostatic contributions to protein−protein binding affinities: Application to Rap/Raf interaction. Proteins: Struct., Funct., Genet. 1998, 30, 407−423. (10) Schreiber, G.; Fersht, A. R. Energetics of protein-protein interactions: analysis of the barnase-barstar interface by single mutations and double mutant cycles. J. Mol. Biol. 1995, 248, 478−86. (11) Jones, S.; Thornton, J. M. Principles of Protein-Protein Interactions. Proc. Natl. Acad. Sci. U. S. A. 1996, 93, 13−20. (12) Schreiber, G.; Haran, G.; Zhou, H. X. Fundamental Aspects of Protein-Protein Association Kinetics. Chem. Rev. 2009, 109, 839−860. (13) Bash, P. A.; Field, M. J.; Karplus, M. Free energy perturbation method for chemical reactions in the condensed phase: a dynamic approach based on a combined quantum and molecular mechanics potential. J. Am. Chem. Soc. 1987, 109, 8092−8094. (14) Rao, S. N.; Singh, U. C.; Bash, P. A.; Kollman, P. A. Free energy perturbation calculations on binding and catalysis after mutating Asn 155 in subtilisin. Nature 1987, 328, 551−554. H


Article

Journal of Chemical Theory and Computation (36) Iwahara, T.; et al. Molecular characterization of ALK, a receptor tyrosine kinase expressed specifically in the nervous system. Oncogene 1997, 14, 439−449. (37) Sun, H.; Li, Y.; Li, D.; Hou, T. Insight into crizotinib resistance mechanisms caused by three mutations in ALK tyrosine kinase using free energy calculation approaches. J. Chem. Inf. Model. 2013, 53, 2376−89. (38) Johnson, T. W.; Richardson, P. F.; Bailey, S.; Brooun, A.; Burke, B. J.; Collins, M. R.; Cui, J. J.; Deal, J. G.; Deng, Y.-L.; Dinh, D.; Engstrom, L. D.; He, M.; Hoffman, J.; Hoffman, R. L.; Huang, Q.; Kania, R. S.; Kath, J. C.; Lam, H.; Lam, J. L.; Le, P. T.; Lingardo, L.; Liu, W.; McTigue, M.; Palmer, C. L.; Sach, N. W.; Smeal, T.; Smith, G. L.; Stewart, A. E.; Timofeevski, S.; Zhu, H.; Zhu, J.; Zou, H. Y.; Edwards, M. P. Discovery of (10R)-7-Amino-12-fluoro-2,10,16trimethyl-15-oxo-10,15,16,17-tetrahydro-2H-8,4-(metheno)pyrazolo[4,3-h][2,5,11]-benzoxadiazacyclotetradecine-3-carbonitrile (PF06463922), a Macrocyclic Inhibitor of Anaplastic Lymphoma Kinase (ALK) and c-ros Oncogene 1 (ROS1) with Preclinical Brain Exposure and Broad-Spectrum Potency against ALK-Resistant Mutations. J. Med. Chem. 2014, 57, 4720−4744. (39) Friboulet, L.; Li, N.; Katayama, R.; Lee, C. C.; Gainor, J. F.; Crystal, A. S.; Michellys, P. Y.; Awad, M. M.; Yanagitani, N.; Kim, S. The ALK inhibitor ceritinib overcomes crizotinib resistance in nonsmall cell lung cancer. Cancer Discovery 2014, 4, 662−673. (40) Sakamoto, H.; Tsukaguchi, T.; Hiroshima, S.; Kodama, T.; Kobayashi, T.; Fukami; Takaaki, A.; Oikawa, N.; Tsukuda, T.; Ishii, N.; Aoki, Y. CH5424802, a Selective ALK Inhibitor Capable of Blocking the Resistant Gatekeeper Mutant. Cancer Cell 2011, 19, 679−690. (41) Sun, Z. S.; Yan, Y. N.; Yang, M. Y.; Zhang, J. Z. H. Interaction entropy for protein-protein binding. J. Chem. Phys. 2017, 146, 124124. (42) Song, J.; Qiu, L.; Zhang, J. Z. H. An efficient method for computing excess free energy of liquid. Sci. China: Chem. 2018, 61, 135. (43) Qiu, L.; Yan, Y.; Sun, Z.; Song, J.; Zhang, J. Z. H. Interaction entropy for computational alanine scanning in protein−protein binding. WIREs Comput. Mol. Sci. 2017, e1342. (44) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (45) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N· log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089−10092. (46) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n -alkanes. J. Comput. Phys. 1977, 23, 327−341. (47) Onufriev, A.; Bashford, D. A. Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins: Struct., Funct., Genet. 2004, 55, 383−394. (48) Ben-Shalom, I. Y.; Pfeiffer-Marek, S.; Baringhaus, K. H.; Gohlke, H. Efficient Approximation of Ligand Rotational and Translational Entropy Changes upon Binding for Use in MM-PBSA Calculations. J. Chem. Inf. Model. 2017, 57, 170−189.

I


Computational Alanine Scanning with Interaction Entropy for Protein

Recommend Documents