Most residues in the active sites of CDK2, CDK4 and CDK6 are remarkably conserved (Figure 2). A key difference is the presence of a histidine residue in CDK4/6 (His95CDK4 and His100CDK6) while CDK2 comprises a phenylalanine (Phe82) in the equivalent position. The His95CDK4/His100CDK6 side-chain is in a position where it potentially can donate or accept a H-bond from an inhibitor. Other differences are in Val96CDK4 and Val101CDK6 corresponding to Leu83CDK2. This residue is capable of forming H-bonds to inhibitors with both backbone NH and carbonyl group, but as its side chain is pointing away from the binding site and is not in direct contact with inhibitors the Val/Leu variation appears to be less relevant for selectivity. Other differences in the binding site are residues Thr120CDK4 and Thr107CDK6, these threonines correspond to Lys89CDK2. The negatively charged residues Asp97CDK4 and Asp102CDK6 have His84CDK2 in the equivalent position, and finally glutamate Glu144CDK4 is corresponding to Gln131CDK2 and Gln149CDK6 ?the latter being the only position where CDK4 and CDK6 have different residues. Interestingly, in all three of these positions CDK4 gains a negative charge relative to CDK2. The potential role of charge as a determinant of CDK4 inhibitor specificity has been pointed out originally by McInnes et al. [40] and more recently by Mascarenhas et al. [41]. In this work, we have studied this example of charge-determined protein-ligand interactions using a variety of methods from the molecular modelling and drug design fields. The binding of inhibitors to protein receptors with high affinity and specificity is central to structure-based drug design applications. The quest for the calculation of binding affinities remains one of the main goals of modern computational biophysical methods [42?8]. The most accurate methods for calculating binding free energies are based on molecular dynamics simulations which predict the physical properties of the protein-ligand complexes based on atomistic structural models. The energetic consequences of small structural changes in inhibitor complexes have been successfully studied using thermodynamic integration
[49?4]. An added benefit of TI calculations, as compared to empirical ligand docking algorithms is that the former include accurate estimates of binding entropy as well as enthalpy, based on rigorous statistical thermodynamics. In this work, we specifically address the contribution of the positive charge of fascaplysin (FAS, see Figure 1A) to selectivity by applying thermodynamic integration calculations. In silico, fascaplysin can be modified easily by the iso-electronic substitution of the positively charged nitrogen to a charge neutral carbon atom, resulting in a compound, which for clarity and simplicity we refer to as carbofascaplysin (CRB, indolo[1,2-c]-fluorene-9-one, see Figure 1B). By calculating the energetic effect of this substitution for the protein-inhibitor complexes of both CDK2 and CDK4, we can quantify the impact of the positive charge of fascaplysin on its specificities towards CDK2 and CDK4.
Methods Structure preparation, homology modelling and ligand docking
The X-ray crystal structure of CDK2 in its active form (PDBID: 1FIN [55]) was used as starting structure for molecular modelling. For CDK4 ligand docking, molecular dynamics simulations and free energy calculations a different strategy was employed. As the five published X-ray structures for CDK4 [38,39] have sections missing and are in an inactive state, a `hybrid’ model was constructed in modeller v9.1 [56]. In this `hybrid’ model the core of the structure is based on the experimentally solved CDK4 structure 2W96, but information from CDK2 (PDB-ID: 1FIN) is used for missing regions and the positioning of the activation loop. ProSa-web [57] and WhatCheck [58] were used for model validation. The coordinates of the `hybrid model’ are provided as supporting information (File S1). The inhibitor fascaplysin (FAS) and the hypothetical compound `carbofascaplysin’ (CRB) were built in Hyperchem 8.0 [59], and energy minimised prior to docking. Ligand docking calculations were performed in GOLD v4.1 [60] using the ChemScore scoring ?function with the binding site defined with a radius of 12 A around the backbone-N of Leu83 for CDK2 and Val96 for CDK4. PyMOL and VMD were used for molecular visualisation. [61,62]
Molecular dynamics simulations
All molecular dynamics simulations were performed using the Amber 10 package [63] with the ff99SB [64] force field for proteins and GAFF [65] for ligands. RESP partial charges for the two ligands FAS and CRB were derived using GAUSSIAN03 [66] at the HF/6-31G* theory level and the antechamber program. CDK2 and CDK4 were solvated in a cubic solvent box so that the distance between every solute atom and the box boundary was at ?least 12 A and neutralised by adding counter ions. Water molecules were treated using the TIP4P-Ew water model, a reparameterization of TIP4P [67] with Ewald summation [68]; five buried crystal waters for CDK2 (PDB-ID 1FIN) and four crystal waters present in equivalent positions in CDK4 were kept in the simulations. Before the simulations the systems were energy minimized, initially by steepest descent followed by conjugate gradient minimization. Then the energy-minimised complexes were equilibrated by 50 ps heating from 0 K to 300 K followed by 50 ps density equilibration, both with weak restraints for all ?protein residues (2.0 kcal/mol A22) and 500 ps constant pressure equilibration. Production runs were run for 5 ns with 2 fs time steps. For all simulations the SHAKE algorithm was used to constrain bonds between hydrogens and heavy atoms [69], Langevin dynamics were used for temperature control. Amber
