Binary and Ternary Complexes of Cucurbit[8]uril with Tryptophan, Phenylalanine, and Tyrosine: A Computational Study.

Selective binding of amino acids, peptides, and proteins by synthetic molecules and elucidation of the geometry and dynamics of the resulting complexes and their strengths are active areas of contemporary research. In recent work, we analyzed via molecular dynamics (MD) simulations the complexes formed between cucurbit[7]uril (CB7) and three aromatic amino acids: tryptophan (W), phenylalanine (F), and tyrosine (Y). Herein, we continue this line of research by performing MD simulations lasting 100 ns to investigate the formation, stabilities, binding modes, dynamics, and specific host-guest noncovalent interactions contributing to the formation of the binary (1:1) and ternary (2:1) complexes in aqueous solution between W, F, and Y amino acids and cucurbit[8]uril (CB8). All complexes were found to be stable, with the binding in each complex dominated by one mode (except for the F-CB8 complex, which had two) characterized by encapsulation of the aromatic side chains of the amino acids within the cavity of CB8 and the exclusion of their ammonium and carboxylate groups. Using the molecular mechanics/Poisson-Boltzmann surface area method to estimate the individual contributions to the overall free energies of binding, results revealed that the key role is played by the amino acid side chains in stabilizing the complexes through their favorable van der Waals interactions with the CB8 cavity and the importance of favorable electrostatic interactions between the carbonyl portal of CB8 and the ammonium group of the amino acid. Visual analysis of structures of the ternary complexes indicated the presence of π-π stacking between the aromatic side chains of the included amino acids. The insights provided by this work may be of value for further efforts aiming to employ the recognition properties of CB8 toward amino acids in applications requiring more elaborate recognition of short peptides and proteins.

Introduction

Molecular recognition is a major aspect and the driving force behind many of the interactions involving biological molecules, peptides, and proteins in particular. The selective binding of free amino acids by synthetic receptors in aqueous solution is an important line of research;[1−4] this, however, is complicated by the similarities existing between some of the side chains of the amino acids (AAs).[5]

AAs serve as important biomarkers for many diseases of the human body, and methods have been developed for detection and measurement of AAs in biological samples such as blood or serum.[6−10] Various synthetic molecules with high binding affinity and selectivity toward AAs and peptides have been synthesized, including cyclodextrins,[11] cavitands,[12] pillararene,[13] calixarenes,[14,15]p-sulfonatocalix[n]arene,[16−18] cucurbit[n]urils.[11,19]

Cucurbiturils (CBn, n = 5, 6, 7, 8, and 10) are a family of synthetic macrocycles whose members are distinguished by the number n of glycoluril units comprising their pumpkin-like cavities, with 2n methylene bridges linking the units. All CBn compounds have polar carbonyl portals that can bind cationic species and a hydrophobic cavity capable of binding hydrophobic guests in water.[20] Features of CBn compounds that make them attractive candidates to pursue as receptors of biomolecules include high stability, availability in different sizes, low cost and relative ease of synthesis and functionalization, and high binding affinity and selectivity toward some biomolecules in water, including AAs,[11,19,21−24] peptides,[25−28] and proteins.[29−31]

The unique molecular recognition features of CBn are responsible for their broad impact and wide range of applications, as in the regulation of biological catalysis (e.g., by CB7 and CB8),[32,33] detection and reaction monitoring of AAs, peptides, and proteins,[34] self-assembly and modulation of proteins,[3] enzyme regulation,[35] and so forth.

Investigations of the binding of CBn toward AAs started in 1998 with the calorimetric titrations employed by Buschmann et al. to determine the binding constants and thermodynamic parameters of complexes formed by CB6 with four AAs and four dipeptides.[11] However, CB6 analogues and the larger members—CB7 and CB8—proved to be more successful at forming inclusion complexes with AAs and short peptides. For example, Isaacs and his co-workers synthesized an expanded analogue of CB6[21] and studied its interaction with aromatic AAs in an aqueous buffer using fluorescence titration. Nau and co-workers developed a tandem-enzyme assay incorporating CB7 for the detection of the presence of AAs.[23] Using UV–visible spectroscopy, Tao et al. reported that CB7 formed inclusion complexes with a series of AAs, binding the aromatic AAs more strongly than the rest.[36]

In 2001, Kim et al. reported that CB8 can bind two guests, one of them being methyl viologen (MV);[37] and simultaneously, the CB8–MV complex binds W or Y to form a ternary complex. Urbach et al. verified via ITC measurements that CB8–MV formed ternary complexes with the aromatic AAs[26] and in a later study found that no such complexes were formed for the rest of the AAs.[24]

In a previous study, we used molecular dynamics (MD) simulations to investigate the binding between CB7 and the zwitterionic forms of the three aromatic AAs, namely tryptophan (W), phenylalanine (F), and tyrosine (Y), and investigated the binding of F in different protonation states and in its deaminated and decarboxylated forms, and we found that all three AAs formed stable inclusion complexes with CB7.[38] We also found that the cationic form of F formed a more stable complex than its anionic form; deamination of F reduced complex stability, while decarboxylation did not result in a noticeable change.[38] Zheng et al. also used MD simulation to investigate the binding of tripeptide that contains N-terminal F; their results predict a sequence selectivity with a high binding affinity of CB7 toward tripeptide that contains N-terminal F.[39]

In this work, we carry out MD simulations and quantum mechanical calculations to investigate the binding of CB8 with the aromatic AAs in aqueous solution and analyze the binding modes and dynamics of the resulting complexes. Furthermore, we employ the molecular mechanics/Poisson–Boltzmann surface area (MM–PBSA) and density functional theory (DFT) methods to determine the noncovalent host–guest interactions and estimate their relative contributions to the binding free energies of the complexes and compare them with experiment.

Computational Methods

The starting geometry of CB8 was obtained from the experimental X-ray diffraction data.[40] MD simulations were performed using the CUDA version of the pmemd module of the AMBER 16.0 program.[41−44] A parameterization by Horn for zwitterionic amino acids[45] was used with the ff14SB force field.[46] GAFF (general Amber force field) parameters were used for CB8.[47] The restrained electrostatic potential (RESP) charges were calculated from the electrostatic potential using an ab initio (HF/6-31G*) method.[48] All systems were solvated in a truncated octahedral periodic box of TIP3P water molecules.[49] The nonbonded interaction cutoff distance was set to 10.0 Å. Prior to starting the MD simulation, each system was subjected to energy minimization using the conjugate gradient algorithm, followed by gradual heating up to 298.15 K for 60 ps and equilibration for 20 ns at 298.15 and 1.0 atm. During the minimization and production runs, the particle mesh Ewald method was employed to treat the long-range electrostatic interactions under periodic boundary conditions.[50] A time step of 2 fs was used in all simulations, updating the nonbonded pair list every 25 steps and restraining all bond lengths involving hydrogen atoms with the SHAKE Algorithm.[51] Production runs were carried out for 120 ns at 298.15 K and 1 atm. Trajectories were analyzed with the CPPTRAJ[52] module of AMBER 16.0.

Visualization of the molecular structures was done using VMD 1.9.3 and 1.2.5 UCSF ChimeraX.[53,54] Hydrogen bonds (HBs) were estimated using a cutoff distance of ≤3.2 Å and a donor–hydrogen–acceptor cutoff angle of ≤120°. The binding free energies, (ΔG), were estimated with the MM–PBSA[55] method via the MMPBSA.py script[56] as implemented in AMBER 16.0, applied to 5000 frames using the protocol described previously.[57]

Quantum calculations were performed using the GAUSSIAN 16 package;[58] the structure of the CB8 and the AAs and their corresponding complexes were first geometry-optimized by the DFT method using the B3LYP function[59,60] and 6-31G* basis set for all atoms and applying the implicit universal solvation model based on density.[61] Vibrational frequency calculations were conducted and resulted in no negative eigenvalues to verify the correspondence to a real minimum. The Minnesota M06-2X functional[62] with the 6-31G** basis set was applied for the free energy calculations. The free energy of complexation in water is obtained by eq , in which ΔE is the gas-phase binding energy calculated from the difference in electrostatic energy between the complex and the free host and guest(s) molecules. ΔGcorrec is the summation of correction from energy to free energy in the rigid rotor harmonic oscillator approximation. ΔGsolv is the free energy term.[63]

Results and Discussion

The average structures of the 1:1 complexes as computed from the corresponding 100 ns trajectories are shown in Figure , which confirm the formation of binary inclusion complexes, with the side chain of each AA encapsulated within the CB8 cavity. Inspecting the structures of Y and W binary complexes reveals the presence of favorable ion–dipole interactions and hydrogen bonding between their ammonium groups and the carbonyl portal of CB8, while their carboxylate groups are seen to be exposed to the bulk water, maintaining a distance from the CB8 portal and minimizing their mutual repulsions with it. The abovementioned information is also clear from inspecting the dynamics of each AA within its binary complex, as depicted by the superposition of 20 snapshots extracted at equal time intervals from its trajectory and displayed in Figure .

Figure 1

MD-average structures and dynamics of the AA-CB8 complexes.

Cluster analysis of the F–CB8 complex trajectory revealed the existence of two binding modes (Figure ): the first (occurring with probability 57%) is similar to the binding modes of W and Y, while in the second, F is fully encapsulated, with the carboxylate group interacting with water molecules trapped in the CB8 cavity (Figure ) and the ammonium group still interacting with the CB8 portal. Presumably, the absence of a second binding mode for W and Y is due to the conformational restriction imposed by HBs formed by the NH and OH groups of their side chains with the CB8 portal, hindering the reorientation of W and Y within the cavity. Further evidence for the full encapsulation of the carboxylate group of F inside the cavity of CB8 comes from the probability distributions of the distances between the carboxylate group of the AAs and the center of mass of CB8 in binary complexes (Figure ).

Figure 2

F–CB8 complex (second binding mode) with the water molecules inside the cavity.

Figure 3

Probability distributions of the distance between the carboxylate group of AAs and the center of mass of CB8 in binary complexes.

Further insights into the conformational dynamics of the binary complexes can be gained from the dynamics of the ammonium group of the AA in each complex, as represented by the distance between the nitrogen atom in the ammonium group and the center of mass of CB8, as depicted in Figure , which reveals a distribution with approximate rotational symmetry around the axis of CB8 and a most probable value of ∼3.3 Å. The probability curves for W and Y exhibit a second peak (around 6 Å), indicating larger conformational freedom for these AAs.

Figure 4

(a) Top and side views of the positions of the nitrogen atom in the ammonium group of the AAs in binary complexes. (b) Probability distributions of the distances between the ammonium nitrogen of AAs and the center of mass of CB8 in binary complexes.

MD simulations of the homogenous ternary complexes were performed for two starting geometries involving two relative initial orientations (parallel and antiparallel) of the two AAs within the CB8 cavity (Scheme ) and revealed that only the antiparallel orientation resulted in stable inclusion complexes, in agreement with the X-ray structure of the Y–CB8 ternary complex.[64] The discussion given below will be limited to the stable complexes. The average structures (Figure ) depict the ammonium groups lying exterior to the CB8 cavity and close to the carbonyl portal and reveal the presence of π–π stacking between the aromatic side chains in the studied AAs.

Scheme 1

Parallel (Left) and Antiparallel (Right) Orientations of AAs within the Ternary Complexes

Figure 5

Average structures (top) and dynamics (bottom) of the ternary complexes.

Similar to binary complexes, the dynamics of the ammonium group in each AA within its ternary complex (Figure ) reveals rotational symmetry around the CB8 axis and shows that the probability curves of W1 and W2 within their ternary complex display two peaks at approximately the same distance but with different heights, as opposed to the identical curves in the cases of Y and F.

Figure 6

(a) Top and side views of the dynamics of the ammonium groups in the ternary complexes. (b) Probability distributions of the distance between the ammonium group and the center of CB8.

Figure displays the probability distributions of the angle between the two planes of the aromatic groups within the ternary complexes and the distance between their aromatic centroids and the corresponding density plots. Evidence for π–π stacking comes from the averages for all complexes of the angle and distance (around 13° and 3.8 Å, respectively), which is further confirmed by the density plots.

Figure 7

Probability distributions of the angle between the two planes of the aromatic part in amino acids (a) and the distances between centers of mass of the aromatic parts of amino acids in ternary complexes (b) and their density dot plots (c).

For each complex, the average numbers of HBs formed by each AA with the carbonyl portals of CB8 and with water are listed in Tables and 2, respectively. The ammonium group in each AA in the 1:1 and 2:1 complexes forms ∼1 HB with the carbonyl portal. The contribution of the side chains in W and Y to H-bonding is slightly higher for each AA in 2:1 than in 1:1 complexes. Table shows that complexation results in loss of AA–water H-bonding. The reduction in H-bonding by each ammonium group in 2:1 complexes is roughly the same as in 1:1 complexes for F and Y but is slightly higher for the 2:1 complex of W. The carboxylate group suffers no loss in H-bonding upon complexation, except for the 1:1 complex of F, because one of the two binding modes involves encapsulation of the carboxylate group within the cavity; however, the loss is partially compensated by hydrogen bonding to the water molecules trapped in the cavity (Figure ).

Table 1

Average Numbers of Intermolecular Hydrogen Bonds Between the Studied AAs and CB8 Carbonyl Groups in 1:1 and 2:1 Complexes

	W		F	Y
	NH3+	side chain	NH3+	NH3+	side chain
1:1	1.00	0.48	0.90	0.94	0.35
2:1	1.25	0.76	0.97	0.95	0.54
	1.09	0.68	0.98	0.96	0.54

Table 2

Average Numbers of AA-Water HBs

	W			F		Y
	NH3+	COO–	NH	NH3+	COO–	NH3+	COO–	OH
free	2.36	4.63	0.91	2.49	4.79	2.44	4.77	1.49
1:1	1.60	4.67	0.24	1.64	3.79	1.67	4.84	0.86
2:1	1.20	4.74	0.04	1.53	4.85	1.53	4.85	0.39
	1.18	4.70	0.03	1.54	4.84	1.54	4.86	0.38

The MM–PBSA estimates of the (favorable) free energies of binding of the AAs in their binary complexes are listed in Table , with ΔEvdW providing the largest contribution presumably coming largely from the interaction between the side chains of the AAs and the cavity of CB8. Factors contributing to the differences in ΔEvdW among the different AA binary complexes include differences in the sizes of the side chains and the degree of AA encapsulation and orientation within the cavity. The electrostatic contribution, ΔEELE, is favorable for all complexes and includes the ion–dipole interactions between the carbonyl portal of CB8 and the ammonium (favorable) and carboxylate groups (unfavorable) of the AAs. The free energy of solvation is unfavorable for all three AAs in binary complexes, highest for F. One would expect a higher ΔGSOLV for W and Y, based on the desolvation of their polar side chains. The fact that F had the highest ΔGSOLV instead could be understood by recalling that F had a significant probability of binding in a mode where its carboxylate group is fully encapsulated (Figures and 3).

Table 3

MM–PBSA Estimates of the Binding Free Energies of AA–CB8 Complexes and Their Decompositions, All in kcal mol–1a

	ΔEvdW	ΔEELE	ΔGPB	ΔGNP	ΔGsolv	ΔG
W	–23.0	–16.2	26.6	–3.2	23.4	–15.8 (–15.6)
F	–22.3	–16.0	27.6	–3.0	24.6	–13.7 (−14.5)
Y	–19.5	–14.6	23.8	–2.9	20.9	–13.2 (−14.2)

Values in parenthesis correspond to the MM-PBSA estimates of the binding free energies of AA–CB7 complexes.[38]

Table lists the MM–PBSA estimates for the overall binding free energies of the ternary complexes. ΔEVDW and ΔEELE were favorable for all ternary complexes. The complexes of W and Y had a higher ΔEELE than F due to the HBs formed between their side chains and the carbonyl portal. W had the highest favorable ΔEVDW due to the larger size of its side chain. The results show that CB8 can form stable ternary complexes with the aromatic AAs. Furthermore, the free energy of adding the second AA to the binary complex (ΔG1:1→1:2 = ΔG1:2 – ΔG1:1) is comparable to or slightly greater than the free energy of forming the binary complex, presumably because the second guest, in the antiparallel relative orientation is positioned to simultaneously form favorable host–guest (with CB8) and guest–guest (with the first AA) interactions.

Table 4

MM–PBSA Estimates of the Binding Free Energies of Ternary Complex AA2–CB8 Complexes and Their Decompositions, All in kcal mol–1a

	ΔEvdW	ΔEELE	ΔGPB	ΔGNP	ΔGsolv	ΔG	ΔG1:1→1:2
W	–47.1	–36.7	58.8	–6.1	52.7	–31.1	–15.3
F	–40.0	–30.3	47.1	–5.5	41.5	–28.8	–15.1
Y	–41.9	–35.0	53.4	–5.7	47.7	–29.2	–16.0

Free energy of adding the second AA to the 1:1 complex.

Employing UV–visible spectroscopy, Cong et al. found 1:1 stoichiometry with CB8 for W and Y and 2:1 stoichiometry for F,[36] while Rajgariah and Urbach using mass spectrometry reported 2:1 complexes for W, F, and Y. Using ITC, they arrived at similar conclusions for W and F; however, for Y, although they observed a binding transition, they were not able to obtain an accurate determination of the stoichiometry.[24] Another study reported an X-ray crystal structure for the 2:1 Y complex.[64] Our simulations confirm the formation of stable ternary complexes for all three AAs

In general, the obtained binding free energy values for the binary complexes with CB8 are lower than the ones reported for CB7,[38] except for W; see Table . For example, the binding free energy values for F and W are higher for CB7 by ∼1 kcal mol–1, which may be attributed to size complementarity[65] and higher energy gain due to the release of high-energy water from the cavity.[66] In contrast, the free energy of adding the second AA to the CB8 binary complex is higher than the 1:1 free binding energy reported with CB7, which can be explained as the more favorable interaction energy with the cavity (additional π–π interaction with the precomplexed AA).

Optimized geometries of the binary and ternary complexes as obtained from DFT calculations (Figure ) show qualitative agreement with their corresponding counterparts obtained from MD simulations. This agreement applies to both geometries of the binary complex formed by F. The results of the DFT calculations of the free energy of binding are shown in Table , which showed a good correlation with the energy obtained from the MM–PBSA, as shown in Figure .

Figure 8

DFT-optimized structures of binary and ternary complexes.

Table 5

Free Energies of Complex Formation (kcal mol–1) Computed Through DFT

	W	F	Y
1:1	–9.6	–9.8, −12.1a	–6.2
2:1	–17.5	–16.4	–18.3

The two values are for the two orientations of F in Figure .

Figure 9

Plot of ΔG values obtained using the MM-PBSA method vs DFT calculations.

Conclusions

We have employed MD simulations to obtain a detailed picture of the geometry and dynamics of the binary and ternary complexes formed between CB8 and the aromatic amino acids (W, F, and Y). We have verified that all binary and ternary complexes with the three AAs were stable, and we used the MM–PBSA method to estimate the contributions of the noncovalent interactions to the binding free energies. Our results show that the complexes are stabilized by the van der Waals interaction between the side chain of the amino acid and the CB8 cavity and the electrostatic interaction with the carbonyl portal of CB8. Furthermore, π–π stacking was observed in ternary complexes. This work should prove useful to efforts aiming to better understand the use of CB8 in binding peptides and proteins.