Identification of productive inhibitor binding orientation in fatty acid amide hydrolase (FAAH) by QM/MM mechanistic modelling {

Modelling of the mechanism of covalent adduct formation by the inhibitor O -arylcarbamate URB524 in FAAH shows that only one of the two possible inhibitor binding orientations is consistent with the experimentally observed irreversible carbamoylation of the nucleophile serine: this is a potentially crucial insight for designing new covalent inhibitors of this promising drug target.

catalytic site in two possible orientations, both of which place the carbamic group close to Ser241. In the first one (I), the m-biphenyl moiety of URB524 occupies the acyl chain binding (ACB) channel, while in the second orientation (II) the cyclohexyl ring occupies the ACB channel and the aryl group is placed in the cytoplasmic access (CA) channel. [10][11][12] Although the inhibitor potency of biphenyl-3-yl derivatives with long and lipophilic N-alkyl groups suggest that this moiety should be placed in the ACB channel (orientation II), 8 traditional computational tools employed in drug discovery such as docking and scoring (e.g. by classical interaction energies) failed to discriminate clearly between these two binding orientations. 10,11 Here we apply a combined quantum mechanics/molecular mechanics (QM/MM) 13 approach to model the entire pathway for carbamoylation of FAAH by URB524. This QM/MM method has previously been validated for FAAH. [4][5][6] The calculations clearly show that the carbamoylation in orientation II is energetically preferred, thus identifying II as the productive binding mode. This provides a theoretical basis for SAR interpretation of URB524 analogues. Furthermore, the proposed mechanism, modelled here for the first time, provides significant insights for inhibitor design.
Non-covalent complexes between FAAH and URB524 in orientations I and II were built following the docking procedure described in ref. 10. The resulting structures were solvated by a 25 Å radius sphere of TIP3P water molecules and equilibrated by 130 ps of molecular dynamics (MD), similar to the procedure reported in ref. 4, and described in the electronic supplementary information (ESI{). The complexes were then minimized to an energy gradient of 0.01 kcal mol 21 Å 21 . CHARMM (version 27b2) was used for these calculations. 14 A schematic representation of such complexes is reported in Scheme S1 (ESI{).
In the QM/MM modelling, the methylamine group of Lys142, the side chains of Ser217 and Ser241 and the whole inhibitor were treated at the PM3 QM level, 15 while the other atoms were treated with the CHARMM22 MM force field 16 (7511 atoms). The covalent bonds crossing the boundary between the QM and MM regions were treated by introducing three 'HQ' link atoms, 4,6,13 which are included in the QM system (62 atoms in total). The QM/ MM approach here includes bonded and non-bonded interactions between the QM and MM systems and accounts for the essential effect of the protein on the modelled reactions. 4 Van der Waals and bonded interactions were described by MM terms, with standard CHARMM22 parameters for the QM atoms. Electrostatic interactions were treated by including the MM atomic charges (as atomic 'cores') in the Hamiltonian for the QM a Dipartimento Farmaceutico, Università degli Studi di Parma, 43100, Parma, Italy. E-mail: alessio.lodola@unipr.it b system. A nonbonded cut-off of 12 Å was applied and atoms further than 14 Å from the Ser241 hydroxyl oxygen were fixed. With the exception of these boundary restraints, all the other atoms (1235 for I, 1260 for II) were free to move during the calculations. To correct for possible shortcomings in the energetics due to the known limitations of the semiempirical method, highlevel energy corrections 4,6,17 were also applied (see below).
The adiabatic mapping approach, which has been shown to perform well with FAAH 4,6 and to provide results in agreement with a free-energy based method, 5 was used to calculate potential energy surfaces (PESs), generating models of the transition states (TSs) and intermediates along the carbamoylation pathway.
The carbamoylating reaction of the nucleophile Ser241 ( R x was increased in steps of 0.2 Å ; R y , R z , R r and R s were increased in steps of 0.1 Å using the RESD 14 command of CHARMM applying harmonic restraints of 5000 kcal mol 21 Å 22 . Geometry optimization of the structures was performed at each reaction coordinate point, to an energy gradient of 0.01 kcal mol 21 Å 21 . The energy was then computed by a single point calculation, removing the energy contributions due to reaction coordinate restraints. High-level corrections (B3LYP-6-31+G(d)) 18 were applied to the crucial stationary points A-G from the PM3-CHARMM22 PESs for a reliable description of the reaction energetics. 19 The corrected values were obtained by subtracting from the total QM/MM energy the PM3 energy of the isolated QM region and adding the B3LYP energy. The B3LYP-corrected potential energy surfaces consist of the B3LYP (in vacuo) QM energy, the CHARMM22 MM energy and the PM3-CHARMM22 QM/MM interaction energy. Interactions within the enzyme should therefore be described reasonably. 4,6,17,19 Fig . 3 shows the resulting B3LYP/6-31+G(d)//PM3-CHARMM22 energy profiles for the carbamoylation reaction of Ser241 by URB524 in orientation I (pink) and II (blue).
In orientation I, the first step of carbamoylation (activation of Ser241 followed by nucleophilic attack on the carbonyl carbon forming the TI (C)), has an energy barrier of 35 kcal mol 21 . Although stabilized by the oxyanion hole (composed of Ile238, Gly239, Gly240 and Ser241 arranged in a hairpin loop), the TI is much less stable than the reactant complex A (by 29 kcal mol 21 ), as expected for a transient configuration. The TI has hydrogen bonds only from the backbone NH groups of Ile238 and Gly239 to the oxygen of URB524 (Table 1).
Breaking of the C-OAr bond produces E with a very low barrier, indicating that expulsion of the m-biphenate anion and TI formation are effectively concerted. During this process, the carbonyl carbon assumed a planar geometry, while the carbonyl oxygen kept its interaction with the oxyanion hole. The high energy of E is due to steric clashes between the m-biphenate and the ACB channel, as shown by comparing QM/MM to QM (in vacuo) calculations and van der Waals contributions to the QM/MM interaction energy.
A double proton transfer (E-G) terminates the catalytic cycle. Protonation of OAr by Ser217 is concerted with proton transfer between Lys142 and Ser217 and represents the rate-limiting step of the whole process, with a barrier of 44 kcal mol 21 (relative to the inhibitor complex). Steric hindrance also affects the transition state (F) and the final product (G), which is less stable than the starting complex A by 23 kcal mol 21 .
Carbamoylation occurs much more easily in orientation II. The barrier for the formation of the TI (C) (30 kcal mol 21 ) is 5 kcal mol 21 lower than in I. The TI (C) is a transient configuration, and is greatly stabilized by the oxyanion hole (the energy of C is only 15 kcal mol 21 above the inhibitor complex). Electrostatic interactions with backbone NH groups play a major role in stabilizing the system: there are two hydrogen bonds (Ile238, Gly239) and two strong interactions (Gly240, Ser241) between the negatively charged oxygen of URB524 and the oxyanion hole ( Table 1).
Breakdown of the tetrahedral intermediate takes place with a very low barrier, and so is effectively concerted with the first reaction step. The product of the reaction, E, is more stable than the starting structure A by 4 kcal mol 21 , in contrast to the findings for orientation I. This key difference arises from crucial interactions at the active site (Fig. 4). Indeed, when the cyclohexyl ring is placed in the ACB channel (orientation II), it assumes a position allowing the carbonyl oxygen of the carbamoylated  Ser241 to interact better with the oxyanion hole. Moreover, the charged oxygen OAr accepts a short hydrogen bond from Ser217 H2 (OAr-H2 = 1.68 Å ) and is also well positioned to 'feel' the field effect of the positively charged Lys142 (OAr-N = 4.82 Å ) which at this stage of the reaction is in its protonated form. This stabilization is weaker in orientation I as the m-biphenate oxygen, residing in the ACB channel, remains further away from the catalytic triad than in orientation II (OAr-H2 = 3.05 Å and OAr-N = 6.37 Å ). The absence of steric clashes between the m-biphenate and the active site contributes to lowering the energy of E in orientation II.
The third step of carbamoylation (E-G) takes place without a significant energy barrier in orientation II, as protonation of the m-biphenate is favoured by the proximity of Ser217 (see above), which is also well orientated to deprotonate Lys142. The resulting product G is very stable: it is the most stable configuration along the modelled pathway in II (218 kcal mol 21 ), consistent with the experimentally observed irreversible inhibition of FAAH. 20 Carbamoylation of Ser241 is favourable in orientation II, while in orientation I the reaction has a significantly higher barrier, and leads to a highly unstable product: the reaction is highly unlikely to proceed in orientation I. Energy profiles at the PM3-CHARMM22 level give similar conclusions (Fig. S1, ESI{), indicating that these findings are independent of the calculation method. Similar pathways were obtained using other snapshots taken from MD simulation (Fig. S2, ESI{), suggesting that the preference observed for binding orientation II is not affected by the starting conformation.
The rate-limiting step for FAAH carbamoylation in orientation II is TI formation. The energy profile in Fig. 3 shows that TI breakdown is not the rate-limiting step, at least for the O-arylcarbamate URB524. This is consistent with the unexpected observation that inhibitor potency of URB524 analogues is not correlated with the electron-withdrawing effect of substituents modifying the stability of the m-biphenate leaving group. 11 Our results are also consistent with the indication that the cyclohexyl ring of URB524 can effectively mimic the anandamide acyl chain, as its replacement with a long lipophilic chain, resembling that of the substrate oleamide, significantly improves inhibitory potency. 8 Finally, our findings highlight the useful contribution that mechanistic modelling of enzymes can make in identifying productive binding modes for covalent inhibitors and will help in designing new FAAH inhibitors. QM/MM modelling of reaction mechanisms of covalent inhibitors in enzyme targets has the potential to provide detailed information for drug design in cases where traditional docking approaches alone may fail.
A. J. M. thanks BBSRC (with CZC) and EPSRC for support.