Quantitative assessment of protein interaction with methyl-lysine analogues by hybrid computational and experimental approaches.

In cases where binding ligands of proteins are not easily available, structural analogues are often used. For example, in the analysis of proteins recognizing different methyl-lysine residues in histones, methyl-lysine analogues based on methyl-amino-alkylated cysteine residues have been introduced. Whether these are close enough to justify quantitative interpretation of binding experiments is however questionable. To systematically address this issue, we developed, applied, and assessed a hybrid computational/experimental approach that extracts the binding free energy difference between the native ligand (methyl-lysine) and the analogue (methyl-amino-alkylated cysteine) from a thermodynamic cycle. Our results indicate that measured and calculated binding differences are in very good agreement and therefore allow the correction of measured affinities of the analogues. We suggest that quantitative binding parameters for defined ligands in general can be derived by this method with remarkable accuracy.

Fine-tuned regulation of gene expression in eukaryotic cells relies on packaging of DNA into different chromatin contexts.[1] The repetitive unit of chromatin, the nucleosome, is formed by wrapping short stretches of DNA around a proteinacous core of histones (H2A, H2B, H3, and H4). Chemical modification of various histone amino acids determines distinct functional chromatin states of transcriptional activation or repression.[2−4] Different methylation states (mono- (me1), di- (me2), and tri- (me3) methylation) of many lysine residues are of high interest, as specific binding proteins that regulate chromatin structure recognize these.[5,6]

Analysis of the different histone methyl-lysine binding regulatory proteins as well as their functional roles in a chromatin context has been hampered by lack of simple experimental tools that allow introduction of defined histone methyl-lysines into recombinant templates. Native chemical ligation and genetic code expansion using stop codon suppression are in many cases cumbersome. These also allow access to only a limited pool of the many methyl-lysine sites.[7−10] Therefore, methyl-lysine analogues (KC) derived from alkylation of genetically introduced cysteine residues have gained particular interest.[11,12] Here, incorporation of sulfur in the γ-position changes the overall geometry by causing longer bond distances but smaller bonding angles compared to the tetrahedral chemical bond geometry of carbon (Figure 1a,b).

Figure 1

Comparing interaction of methyl-lysine binding domains (BD) with methyl-lysine (Kme) and methyl-lysine analogues (KCme) containing peptides. (a, b) Structural parameters of Kme and KCme. (c) Experimental scheme for measurement and calculation of ΔΔG values in a thermodynamic cycle. The difference in binding free energy between the BD*Kme-pep and BD*KCme-pep complexes can be assessed by alternate routes: experimentally via ΔG2 – ΔG1 or computationally via ΔG4 – ΔG3. In the nonequilibrium thermodynamic integration runs the λ = 0 state of the complex is represented by BD*Kme-pep and the λ = 1 state is represented by BD*KCme-pep.

However, it is unclear how accurately methyl-lysine analogues actually mimic the methyl-lysine binding affinities, a situation that is prevalent in many fields. Often, quantitative assays exist (or are much easier/faster/cheaper to set up) only for analogues rather than for the biomolecule of interest. Whether the analogue is close enough to justify quantitative interpretation of an experiment has so far been an issue of chemical intuition and frequent controversial discussions.

To systematically address this question, we developed, applied, and assessed a hybrid computational/experimental approach that extracts the binding free energy difference between the native ligand (methyl-lysine) and the analogue (methyl-amino-alkylated cysteine) from a thermodynamic cycle (Figure 1c). The result is used to correct the measured affinity of the analogue.

We started by calculating differences in free binding energy using atomistic force field free energy simulations for methyl-lysine and methyl-lysine analogue ligands for histone methyl-lysine binding domains where structural details from crystallization and X-ray analysis are available. Then, we used quantitative isothermal titration calorimetry (ITC) and fluorescence polarization (FP) methods to measure the dissociation constants of methyl-lysine binding domain (BD) histone peptide complexes at thermodynamic equilibrium (Figure 1c). Four interaction pairs of different methyl-lysine binding domains interacting with distinct histone methyl-lysine sites and at different methylation level were initially analyzed: the Chromo domain of Su(Var)205 binding to histone H3 trimethylated on lysine 9 (H3K9me3), the 3xMBT domains of L3MBTL1 binding to histone H4 monomethylated on lysine 20 (H4K20me1), the PHD finger of ING1 binding to histone H3 trimethylated on lysine 4 (H3K4me3), and the pseudo Tudor domain of ICBP90 binding to H3K9me3 (Figure 2).

Figure 2

Structural details of methyl-lysine (derived from the PDB entries, in gray) and corresponding methyl-lysine analogues (in purple with the sulfur in yellow) containing BD*pep complexes as used for the λ = 0 and λ = 1 states of the calculation are depicted for (a) Su(Var)205 Chromo domain/H3K9me3 (PDB code 1kne), (b) L3MBTL1 3xMBT/H4K20me1 (PDB code 2rhy), (c) ING1 PHD finger/H3K4me3 (PDB code 2qic), and (d) ICBP90 pseudo Tudor/H3K9me3 (PDB code 3db3). Histone peptide backbones are represented by ribbon diagrams. Protein domains are given as van der Waals surfaces. Note that the peptide backbone of H4K20me1 (residues 15–25) in the L3MBTL1 3xMBT complex was added in the absence of structural details.

Simulation parameters for the different methyl-lysine and methyl-lysine analogues were obtained according to the generalized amber force field procedure with partial charges derived from quantum mechanical calculations.[13,14] Binding free energy differences were then derived using a non-equilibrium thermodynamic integration protocol in which the free energy change for an alchemical transition from methyl-lysine to methyl-lysine analogues is assessed.[13,14] The free energy calculations clearly indicated that the interaction strengths of the PHD finger of ING1 with H3K4me3 or H3KC4me3 and of the pseudo Tudor domain of ICBP90 binding to H3K9me3 or H3KC9me3 are not considerably affected by changing of the methyl-lysine into a methyl-lysine analogue (ΔΔG values of −0.6 and −1.2 kJ/mol, respectively). However, significant differences in free binding energies were calculated for the Chromo domain of Su(Var)205/H3K9me3 and the 3xMBT domains of L3MBTL1/H4K20me1 interaction pairs. In both cases the methyl-lysine analogues reduced binding by more than 6 kJ/mol compared to interaction with the methyl-lysine containing sequence (Table 1).

Table 1

Measured and Calculated Binding Free Energy Differences (ΔΔG) of Interaction of the Indicated Proteins/Domains with Different Methylated Histone Tail Peptides

protein	PDB	domain	site	methoda	Kd(Kme)[μM]	ΔG(Kme)[kJ/mol]	Kd(KCme)[μM]	ΔG(KCme)[kJ/mol]	ΔΔG(exp)[kJ/mol]	ΔΔG(calc)[kJ/mol]
Su(Var)205	1kne	Chromo	H3K9me3	ITC	0.3 ± 0.0	–37.6 ± 0.2	1.6 ± 0.1	–33.1 ± 0.2	4.5 ± 0.2	7.3 ± 0.4
ING1	2qic	PHD finger	H3K4me3	FP	16.2 ± 4.4	–27.3 ± 0.7	11.4 ± 2.3	–28.2 ± 0.5	–0.9 ± 0.8	–0.6 ± 0.4
L3MBTL1	2rhy	3xMBT	H4K20me1	FP	6.2 ± 0.5	–29.7 ± 0.2	48.2 ± 3.0	–24.6 ± 0.2	5.1 ± 0.3	6.1 ± 0.5
ICBP90	3db3	pseudo Tudor	H3K9me3	FP	0.4 ± 0.1	–36.7 ± 0.5	0.4 ± 0.1	–36.7 ± 0.7	0.0 ± 0.8	–1.2 ± 0.3
53BP1	2ig0	Tudor	H4K20me2	FP	18.1 ± 2.8	–27.0 ± 0.4	b	b	c	2.7 ± 0.4
hHP1α(CBX5)	3dm1	Chromo	H3K9me3							7.3 ± 0.4
hHP1β(CBX1)		Chromo	H3K9me3	FP	1.9 ± 0.5	–32.7 ± 0.7	11.6 ± 3.9	–28.2 ± 0.8	4.6 ± 1.1
hHP1γ(CBX3)	3fdt	Chromo	G9aK165me3							7.3 ± 0.4

Apparent molar dissociation constants (Kd) for the different methyl-lysine binding domains were derived from isothermal calorimetry (ITC) or fluorescence polarization (FP) measurements for methyl-lysine (Kme) and methyl-lysine analogues (KCme) containing peptides.

Not binding.

Not determined.

When we measured the dissociation constants of the same complexes we could verify the differences derived from the theoretical calculations (Table 1 and Supplementary Figures S1 and S2). The differences in free binding energy of the Su(Var)205/H3K9me3 and L3MBTL1/H4K20me1 complexes when comparing methyl-lysine with corresponding methyl-lysine analogues target peptides were 4.5 and 5.1 kJ/mol, respectively. In contrast, the differences for the ING1/H3K4me3 and ICBP90/H3K9me3 interaction pairs were below 1 kJ/mol. Importantly, as Figure 3 shows, the calculated data were in very good agreement with the differences derived from the experimental measurements of the Kd of the same interaction pairs.

Figure 3

Comparison of experimental and calculated ΔΔG values for different BD*Kme-pep and BD*KCme-pep complexes.

The prediction of the binding differences of methyl-lysine and methyl-lysine analogues targets depends, however, on the structural details provided in the structural X-ray analysis. We therefore also investigated the interaction of the Tudor domain of 53BP1 with H4K20me2. While this domain bound the methyl-lysine peptide with a Kd of 18 μM, no binding to the corresponding methyl-lysine analogue containing target was observed. The force field free energy calculation based on the single amino acid visible in the available structure provided here a ΔΔG of only 2.7 kJ/mol. Nevertheless, we found this value to somewhat vary with the sequence context of the methl-lysine or methyl-lysine analogues. Addition of residues flanking the single methyl-lysine for the force field calculations resulted in large variations depending on how and where these sequences were placed in the structure (data not shown). Interestingly, a similar procedure for the 3xMBT domain of L3MBTL1 resulted in converging results, independent of the addition of amino acids flanking the structurally resolved single methyl-lysine residue.

Additional calculations on more detailed structures such as the Chromo domains of hHP1α/H3K9me3 and the non-histone methyl-lysine interaction pair of hHP1γ/G9aK165me3 also indicated large differences in ΔΔG. Although the hHP1β/H3K9me3 complex could not be theoretically analyzed due to lack of a high resolution structure, our measurements nevertheless indicated preferential binding of methyl-lysine over methyl-lysine analogues (Table 1). Together with the data obtained on the Drosophila HP1 ortholog Su(Var)205, these findings suggest that Chromo domains in general are particularly sensitive to methyl-lysine analogues versus methyl-lysine interaction.

Our results caution against simple use of methyl-lysine analogues in quantifying histone methyl-lysine binding domain protein interaction and likely in analyzing the activity of methyl-lysine modifying enzymes. Where detailed structural data are available, our hybrid computational/experimental approach allows the remarkably accurate determination of interaction differences from analogues. In other instances direct experimental comparison of methyl-lysine and methyl-lysine analogues paradigm targets seems crucial. Considering that several histone methyl-lysine binding domain proteins such as the HP1 factors multimerize and are thought to have multivalent nucleosome binding,[15] such differences might amplify and need to be factored into any calculation and modeling of higher order chromatin interactions.[16]

We envision the generic methodology described and tested here will be applicable to a wide range of complex protein–ligand pairs from basic biochemical research to pharmaceutical lead compound development. Often such interactions are not accessible to experimental analysis, e.g., due to difficulties, time constraints, and/or high costs associated with purification of natural sources or chemical synthesis, or accuracy and/or throughput is limited. Here, our hybrid computational and experimental approach will not only provide more accurate results from measurements using analogues. It will also make it possible to shift the point of view and to relax the requirement of stereochemically highly similar analogues, e.g., in favor of easier chemistry. Since the difference in binding affinity is calculated and fully accounted for, significant dissimilarities can be tolerated. Accordingly, we expect that our approach will not only provide more accurate results in many fields but also open up access to systems where so far no consensus analogues are available or accepted by the scientific community.

Methods

Binding Free-Energy Calculations

All simulations were carried out using the Gromacs molecular dynamics package[17−19] (version 4.0.7) and the AMBER99SB force field.[20] Simulations were carried out in explicit solvent with the tip3p water model at 150 mM NaCl.[21] Simulation parameters for the methylated lysine derivatives and their methyl-amino alkylated cysteine analogues were obtained according to the generalized amber force field (GAFF[14]) procedure with partial charges derived from quantum mechanical calculations with Gaussian03 (Hartree–Fock/6-31G* basis set).[22] Hybrid residues representing either the lysine derivative or its sulfur analogue as a function of λ were constructed as described.[13] Binding free energy differences were calculated using a non-equilibrium thermodynamic integration scheme. Each system was sampled for 20 ns at λ = 0 and λ = 1. Sampling was carried out at 298 K using a leapfrog stochastic dynamics integrator, with pressure kept at 1 atm using a Parrinello-Rahman barostat.[23] Electrostatic interactions were calculated at every step with the particle-mesh Ewald method,[24] short-range repulsive and attractive dispersion interactions were described by a Lennard-Jones potential with a cutoff of 1.1 nm, and a switching function was used between 1.0 and 1.1 nm. Dispersion correction for energy and pressure was applied. The SETTLE algorithm was used to constrain bonds and angles of water molecules, and LINCS was used for all other bonds, allowing a time step of 2 fs.[25,26] From the equilibrium trajectories snapshots were taken every 40 ps, and thermodynamic integration runs from λ = 0 to λ = 1 and λ = 1 to λ = 0, respectively, with a switching time of 50 ps were performed from each snapshot. The derivative of the Hamiltonian with respect to λ was recorded at every step. Free energies were subsequently calculated using the Crooks-Gaussian-Intersection method.[27]

Protein Expression and Purification

Plasmids, expression, and purification of His6-53BP1 (residues 1484–1603), GST-ING1 (residues 200–279), His6-L3MBTL1 (residues 197–526), and His6-Su(Var)205 (residues 17–76) were as described previously.[28−31] The coding sequence of the human ICBP90 pseudo Tudor domain (residues 126–285) was cloned into pET16b with an N-terminal His10-tag followed by factor Xa and TEV protease cleavage sites. The coding sequence of human HP1β encompassing an N-terminal His6-tag was cloned into pET11a. ICBP90 and HP1β were expressed in E. coli BL21(DE3)RIL and purified by standard Nickel-NTA chromatography.

Peptides

Peptides were synthesized using Fmoc chemistry on an Intavis Respep XL synthesizer. TentaGel R RAM resin (cap.: 0.19 mmol g–1) served as solid support and the amino acid side chains were protected as follows: Arg(Pbf), Asn(Trt), Asp(OtBu), Gln(Trt), His(Trt), Lys(Boc), Ser(tBu), Thr(tBu), and Lys(Mtt) for orthogonal deprotection of the ε-amino group. Coupling reactions were performed with 2-(1H-benzotriazole-1-yl)-1,1,3,3-tetramethyluroniumhexafluorophosphate (HBTU) as coupling reagent and N-methylmorpholine (NMM) in DMF/NMP as base. 5,6-Carboxyfluorescein was coupled either at the N-terminus or at the ε-amino group of a C-terminal Lys residue using HBTU/hydroxybenzotriazole (HOBt) and NMM in DMF. Pseudoproline dipeptides were used for efficient synthesis of H3-derived peptides. At the indicated positions methylated lysine or cysteine (Kc) residues were introduced. Cysteine was coupled as pentafluorophenyl ester (Fmoc-Cys(Trt)-OPfp) without base to ensure enantiomerization-free introduction of this building block.[32] Peptides were cleaved off the resin with TFA/phenol/triisopropylsilane/H2O (85:5:5:5) for 4 hours. Cysteines were alkylated to produce methyl-lysine analogues as described previously.[12] All peptides were purified by reversed phase C18 HPLC and verified by ESI-MS (see Table 2).

Table 2

H3(1–15)K4me3-Fl	ARTK(me3)QTARKSTGGKA-fluoresceine
H3(1–15)KC4me3-Fl	ARTKC(me3)QTARKSTGGKA-fluoresceine
Fl-H3(1–15Y)K9me3	fluoresceine-ARTKQTARK(me3)STGGKAY
Fl-H3(1–15)KC9me3	fluoresceine-ARTKQTARKC(me3)STGGKA
Fl-H4(12–27Y)K20me1	fluoresceine-KGGAKRHRK(me1)VLRDNIQY
Fl-H4(12–27)KC20me1	fluoresceine-KGGAKRHRKC(me1)VLRDNIQ
Fl-H4(12–27Y)K20me2	fluoresceine-KGGAKRHRK(me1)VLRDNIQY
Fl-H4(12–27)KC20me2	fluoresceine-KGGAKRHRKC(me1)VLRDNIQ
H3(1–15Y)K9me3	ARTKQTARK(me3)STGGKAY
H3(1–15Y)KC9me3	ARTKQTARKC(me3)STGGKAY

Binding Measurements

Fluorescence polarization measurements were performed at 25 °C in 10 mM triethanolamine (pH 7.5), 0.1 mM EDTA, 20 mM NaCl in a 96-well format as described previously using a Plate Chameleon II plate reader (HIDEX Oy).[33] Curves were fitted using least-squares fitting (Kaleidagraph) to the equation F = Fmin + (Fmax – Fmin)·x/(Kd + x). F is the fluorescence polarization signal at concentration x; Fmin is the fluorescence polarization in the unbound state; Fmax is the fluorescence polarization in the bound state; and Kd is the apparent molar dissociation constant. Data were normalized as fraction bound (Fb) using the equation Fb = (F – Fmin)/(Fmax – Fmin). Multiple titration series were averaged after data normalization.

ITC measurements were performed at 25 °C in 10 mM triethanolamine (pH 7.5), 0.1 mM EDTA, 20 mM NaCl on an iTC200 calorimeter (Microcal). Heats of binding reactions were recorded by sequential injection of the binding protein into unlabeled H3 peptides. Raw data were integrated and normalized, and the apparent heat change was plotted using Origin software (OriginLab). For determination of the molar association constant, nonlinear least-squares fitting was performed using a one set of identical binding sites model.