Binding residence time through scaled molecular dynamics: a

Subscriber access provided by REGIS UNIV

Computational Chemistry

Binding residence time through scaled molecular dynamics: a prospective application to hDAAO inhibitors Mattia Bernetti, Elena Rosini, Luca Mollica, Matteo Masetti, Loredano Pollegioni, Maurizio Recanatini, and Andrea Cavalli J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00518 • Publication Date (Web): 19 Oct 2018 Downloaded from http://pubs.acs.org on October 21, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Binding Residence Time through Scaled Molecular Dynamics: a Prospective Application to hDAAO Inhibitors Mattia Bernetti,1,2 Elena Rosini,*3 Luca Mollica,*4,5 Matteo Masetti,1 Loredano Pollegioni,3 Maurizio Recanatini1 and Andrea Cavalli1,2 1

Department of Pharmacy and Biotechnology, Alma Mater Studiorum – Università di Bologna,

Via Belmeloro 6, 40126, Bologna, Italy 2

Computational & Chemical Biology, Istituto Italiano di Tecnologia, Via Morego 30, 16163,

Genova, Italy 3

Department of Biotechnology and Life Sciences, Università degli Studi dell’Insubria, Via J.H.

Dunant 3, 21100, Varese, Italy

4

Istituto Nazionale Genetica Molecolare "Romeo ed Enrica Invernizzi", Via F. Sforza 35, Milan,

20122, Italy

5

Department of Medical Biotechnology and Translational Medicine, Università degli Studi di

Milano, Milano 20129, Italy 1

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 44

ABSTRACT Traditionally, a drug potency is expressed in terms of thermodynamic quantities, mostly Kd, and empirical IC50 values. While binding affinity as an estimate of drug activity remains relevant, it is increasingly clear that it is also important to include (un)binding kinetic parameters in the characterization of potential drug-like molecules. Herein, we used standard in silico screening to identify a series of structurally related inhibitors of hDAAO, a flavoprotein involved in schizophrenia and neuropathic pain. We applied a novel methodology, based on scaled molecular dynamics, to rank them according to their residence times. Notably, we challenged the application in a prospective fashion for the first time. The good agreement between experimental residence times and the predicted residence times highlighted the procedure’s reliability in both predictive and refinement scenarios. Additionally, through further inspection of the performed simulations, we substantiated a previous hypothesis on the involvement of a protein loop during ligand unbinding.

1. INTRODUCTION In recent times, binding and unbinding kinetics have emerged as key concepts in drug discovery.1–3 Indeed, while thermodynamic binding affinity has been widely used as the major predictor of drug potency in vivo, recent studies have pointed to kinetics as a key parameter to consider during lead optimization and preclinical development. Specifically, the focus has shifted towards the inclusion of kon and koff, the association and dissociation rate constants, respectively.1 Of particular interest is the unbinding kinetics, often referred to as residence time 2


Page 3 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(1/koff), which is a measure of the time a drug spends in contact with the target. Unbinding kinetics is often seen as a key parameter to consider in drug discovery programs.2,4 Indeed, while there is no doubt about the significance of the binding affinity as an estimate of drug potency, there are cases in which it is more convenient to integrate this information with residence time (tr). One such example is when the duration of the pharmacological effect plays a substantial role in in vivo efficacy.2 Many experimental approaches to determining tr have been reported,5–8 including surface plasmon resonance (SPR),9 nuclear magnetic resonance spectroscopy (NMR),10 and fluorescence methods.11 In contrast, the literature contains few computational tools that can predict tr with a good level of accuracy. Clearly, computational approaches to tr would considerably improve the efficiency of the optimization strategy, as they routinely do for binding affinity. Molecular dynamics (MD) is emerging as a key tool for drug discovery,12 and tr estimation is a field where MD-related methods are increasingly being applied.8,13,14 The atomic-level determinants underlying binding and unbinding kinetics can in principle be determined via molecular simulations. However, this task can only be accomplished by mapping the full binding/unbinding process, which is far beyond the common limits of ordinary MD simulations. Therefore, viable alternatives have been purposely designed during recent years. Enhanced sampling methods are currently the most effective tools to access rare events, including protein-ligand dissociation. As such, methodologies based on these approaches allow providing estimates of tr. Metadynamics (MetaD),15,16 originally intended to reconstruct free energy surfaces, has been quite extensively explored for this purpose. In this respect, either means to compute kinetic rates directly17 or protocols to prioritize ligands according to their unbinding times18 have been introduced. 3



Page 4 of 44

However, these approaches require careful and specific setup by the user, where the choice of suitable collective variables (CV) to guide ligand binding/unbinding plays the major effort. τRAMD represents an alternative strategy to rank drug candidates according to their computational unbinding times.19 The procedure is based and named after random accelerated MD (RAMD), a biased method that does not require the definition of CVs to enhanced the sampling of rare events. Multiscale schemes have also been explored to estimate rate constants. Taking advantage of a dual regime, Brownian Dynamics (BD) is employed at high protein-ligand separation, while plain MD is used when the two species are in closer vicinity.20 Thus, being the protein-ligand encounter modelled as a rigid diffusion of the two bodies, a remarkable amount of computational time can be saved. As a result, taking advantage of these frameworks and accounting on typically limited computational resources, researchers can manage just one to a few ligands at a time, at most. Moreover, such procedures do not produce results quickly enough for most drug optimization programs. Recently, a more practical and accessible methodology has been introduced to tackle the limitations described above.21 Instead of aiming to calculate Kd or IC50 values, the goal is to rank ligands according to their tr and to distinguish them as faster and slower in terms of unbinding simulation time.21,22 The methodology is based on scaled molecular dynamics (scaled MD),23,24 an enhanced sampling method that improves the exploration of the configurational space by smoothing the potential energy surface (PES) of the entire system. Through scaled MD, several unbinding events can be observed at considerably lower computational costs. By applying the procedure to a series of ligands, the compounds can be assigned an average unbinding time and prioritized accordingly. Thus, although these simulations are out of equilibrium, their 4


Page 5 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interpretation in relative terms is meaningful and reliable. Moreover, valuable indications about the molecular features involved in the (un)binding process can be retrieved. These features can then be examined further via more expensive methods to achieve a more detailed characterization. Herein, we used a virtual screening campaign to identify new potential inhibitors of a pharmaceutically relevant biomolecular target, human D-amino acid oxidase (hDAAO, EC 1.4.3.3). Then, we carried out scaled MD simulations on a small set of analogs from among the identified compounds in order to achieve a kinetic characterization of their binding properties and validate the virtual results. hDAAO is a flavoprotein that catalyzes the oxidative deamination of D-amino acids with excellent stereoselectivity.25 Its pharmacological relevance in humans stems from its role in regulating D-serine levels in several brain regions, from cerebellum to frontal cortex.26,27 D-serine acts as a synaptic co-agonist at the NR1 subunit of the N-methyl-Daspartate receptors (NMDARs), central regulators of synaptic plasticity, learning and memory, and pain sensation. A decreased D-serine concentration is correlated with hypofunction of the NMDAR mediated neurotransmission, a core pathway deficit in schizophrenia.25,28–30 The physiological concentration of D-serine can be restored by inhibiting the DAAO activity.31,32 In the last 10 years, several academic and industrial research groups focused on the therapeutic potential for small molecule inhibitors of hDAAO to treat disorders like schizophrenia or neuropathic pain.33 Hundreds of new hDAAO inhibitors have now been discovered.31,34 Here, by performing multiple scaled MD simulations on bound complexes of hDAAO, five inhibitors were ranked according to their computational unbinding times. The predictions were then compared with experimental thermodynamic (Kd) and kinetic constants (kon and koff). There 5



Page 6 of 44

was a good agreement between the kinetic trends obtained by computational methods and by experimental analyses. Moreover, since multiple unbinding events were sampled, we explored the obtained trajectories to retrieve potential insights into the molecular rearrangements taking place during ligand detachment from the protein. Interestingly, this allowed us to identify two major conformations that can be adopted by a loop located at the entrance of the active site when ligands unbind.

2. RESULTS AND DISCUSSION In the initial effort to identify new potential drug-like inhibitors of hDAAO, a virtual screening campaign was conducted on a set of 916,000 compounds (see Methods for details). In particular, we sought competitive inhibitors targeting the protein active site, where the oxidation of the substrate takes place through reduction of a molecule of FAD. This led to the identification and purchase of 24 compounds (Suppl. Table 1, subgroup A), which were characterized in terms of their binding properties through an activity assay suitable for highthroughput screening (see Methods).31 The inhibition efficiency was investigated by a coupled enzyme assay in 96-well plate format. For all concentration-response curves, the regression equation allowed us to estimate an IC50 value. Of these, six compounds were active in the µM regime (Suppl. Table 1). The IC50 values of the ligands ranged from about 10 to about 500 µM. To improve their potency, we used the SciFinder engine to seek purchasable analogs of these compounds.35 As a result, a further 19 compounds were identified, purchased, and assayed (Suppl. Table 1, subgroup B). As shown in Suppl. Table 1, however, no significant improvements in inhibitory activity were observed. Nevertheless, a group of four compounds (4, 6


Page 7 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


7, 9, 15, see Figure 1) was particularly interesting because subtle differences in chemical substituents were responsible for significantly different IC50 values (see Table 1). Moreover, we added the scaffold without any substituent (100, see Figure 1) to the group of the considered compounds because it was previously reported to be active towards hDAAO.36 We applied the procedure described in previous studies to prioritize these compounds in terms of residence times.21,22 Here, with no a priori knowledge of their experimental kinetic profile, the method was applied in a prospective manner.21 Notably, the crystal structure of hDAAO in complex with the specific scaffold of the compounds that we selected was not present in the Protein Data Bank (PDB), which made the application even more challenging. However, a crystallographic structure with a closely related molecule was available (PDB code 3CUK), presenting the same pyrrole-2-carboxylic acid moiety responsible for the driving interactions in the binding site. Therefore, this crystal structure was used as a template to locate the ligands shown in Figure 1 inside the binding pocket (see Methods). As depicted in Figure 2A, the carboxylic group forms a salt bridge with the positively charged side chain of Arg283 and a hydrogen bond with the side chain hydroxyl group of Tyr228, while the pyrrole nitrogen acts as a hydrogen bond donor to bind the backbone carbonyl of Lys313. Notably, in different crystal complexes, there is conservation of these interactions between hDAAO and ligands made of double-ring aromatic scaffolds. Being all these features maintained in the compounds chosen for the present work, we assumed the generated binding modes to be reasonable and to most likely correspond to the real poses. Indeed, we highlight that the initial configuration is expected to be a critical aspect for the methodology, since it can affect the dynamics of unbinding. Another classic feature characterizing the binding mode of most of the known inhibitors,31,32,37 depicted in 7



Page 8 of 44

Figure 2B, is their sandwich-like placement between the isoalloxazine ring of the FAD cofactor and the phenol ring of Tyr224. The latter belongs to the so-called active site lid, a loop located at the entrance of the binding pocket, which hinders its accessibility (Figure 2C). For each compound shown in Figure 1, starting from the bound complexes (Figure S1), 20 scaled MD simulations were performed to reproduce ligand unbinding from the active site of hDAAO. Runs were stopped once release took place or if no detaching of the compounds was observed within a computational time of 100 ns. A scaling factor of 0.45 was used, which allowed the observation of unbinding events between 3 and 100 ns. Aggregating the entire set of trajectories, a total simulation time of about 3.4 µs was gathered, including 82 unbinding events recorded out of the 90 runs carried out on the 5 ligands. The distribution of the computational time of unbinding for each ligand is reported in Figure 3: the unbinding times obtained were recorded and, according to the procedure followed in previous works,21,22 the corresponding average, median, standard deviation, and standard error values were calculated over 18 runs (the fastest and the slowest unbinding events for each series of simulations for each ligand were excluded) performed for each ligand. Moving from ligand 9 to 100, there was a shift in the distributions coupled with their broadening (as indicated by the standard deviation reported in Table 2), with a gradual reduction of the population of fast dissociation events and an increase of that for the slow ones. Despite the small structural differences, the ligands were characterized by considerably different average unbinding times. When considered more comprehensively, two major regimes could be distinguished. Ligands 9, 4 and 7, with average unbinding times of about 18, 26 and 28 ns respectively, can be grouped in a fast regime, and ligands 15 and 100, approaching about 60 ns, can be grouped in a slow regime. Notably, while the first group is characterized by ligands 8


Page 9 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bearing methyl or ethyl substituents on the pyrrol ring, this is not the case for ligands 15 and 100. This suggests that the ligand dissociation is favored by hindrance caused by chemical groups inserted in that position of the molecules, pointing towards the bottom of the binding pocket. Notably, previous structure-activity relationship (SAR) studies suggested how the presence of bulkier substituents, such as methyl and ethyl groups, reduced the inhibition efficiency of hDAAO.38 The consistency between our results and the behavior of known drug-like molecules inhibiting hDAAO supported the reliability of the initial theoretical complexes constructed for the 5 ligands. Our findings can be considered in the context of the recently introduced concept22 of structure-kinetics relationships (SKR), i.e. an accurate description of the structural determinants responsible for variations of the residence time that should be of paramount importance when designing new compounds with improved in vivo profiles. The presented method can separate the available ligands into two distinct kinetic classes (in terms of their residence times) and, at the same time, provide a structural/chemical perspective about kineticsbased drug design even based on subtle differences between ligand substituents. The steric hindrance of the aliphatic groups in position 3 of the indole moiety seems to be the most avoidable choice in the synthetic route in order to obtain ligands with increased tr. We performed an in silico prediction of residence times of congeneric ligands of hDAAO, based on the assumption that the ligand pose in hDAAO binding site was retained with respect to a reference crystal structure in complex with pyrrole-2-carboxylic acid. In this scenario, the sensitivity of the method was tested for a series of closely chemically related compounds (i.e. presenting only minor chemical substitutions on the same scaffold), as is routinely the case during the optimization phase of drug discovery projects. Notably, these substituent groups did 9



Page 10 of 44

not introduce the possibility of achieving any additional hydrogen bonds or salt bridges. To the best of our knowledge, this is the first report of prospective kinetic predictions carried out with scaled MD simulations as they might feasibly be applied in “real-life” industrial scenarios. To validate our predictions, the compounds were assayed for binding, inhibition, and kinetics of association and dissociation using recombinant human enzyme (Table 1). All the tested compounds inhibited hDAAO activity. A full inactivation was apparent for ligands 4, 7, and 15, with the latter also showing the higher potency (IC50 = 5.4 ± 0.3 µM, Table 1). The binding affinity of each ligand was estimated by titrating the oxidized form of the enzyme and by following the changes in the visible absorption spectrum of the FAD cofactor, in the absence of substrate. In all cases, there was a perturbation at above 500 nm indicative of a charge-transfer complex formation.39 Moreover, the kinetics of binding/unbinding was monitored by a stoppedflow apparatus. The time course of the complex formation was followed by rapid mixing of hDAAO and inhibitor; the observed rate constant (kobs) plotted vs. inhibitor concentration yields a linear behavior (Figure 4). The rates of inhibitor binding ranged from 1.2 x 105 to 8 x 105 M-1 s1

for all compounds (Table 1) with the exception of ligand 100 which binds relatively slowly (kon

≈ 3 x 104 M-1 s-1), yielding a dynamic Kd value (≈ 37 µM, as calculated by the koff/kon ratio) similar in magnitude to that estimated by static titration (44 ± 2.5 µM, Table 1). If we assume that, in most cases, more potent compounds will have longer residence times, then the proportionality between the measured IC50 values and the predicted unbinding times (Figure 5A) ensures a fair robustness of the proposed method. However, to provide a proper figure of merit, the predicted dissociation times need to be compared to the experimental koff values, which in turn correlate inversely to the residence time. The computational unbinding 10


Page 11 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


times for ligands 9, 4, 7, and 15 reported in Table 2 are in good agreement with the koff values shown in Table 1. In particular, increasing tr predicted by the simulations corresponded to decreasing experimental koff values. Thus, there was excellent agreement for these compounds in terms of ranking (see Figure 5B). Surprisingly, notwithstanding a very subtle structural difference between 9 and 7, no binding was observed for the former in our experiments. This behavior could be due to an increase in steric hindrance when replacing a fluorine with a methyl group, coupled to a change in local polarity of the molecule. Indeed, typical inhibitors of hDAAO are double aromatic ring scaffolds that almost completely fill the active site. The only exception was 100, for which the highest tr was estimated from simulations. However, the experimental value of koff for 100 is very close to 4. Indeed, Figure 5B shows how the correlation between predicted and experimental data would worsen remarkably when including 100 (blue point in the figure). It is interesting to notice how the distribution of the dissociation times for this ligand, as reported in Figure 3, shows a separate tail at very high values, close to 100 ns. These events have a significant weight on the average, shifting the outcome towards higher values. Indeed, this behavior is also apparent in Table 2, where, conversely to all the other compounds, the median differs from the average by more than 10 ns, indicating that most of the samples are found at lower values. On the one hand, this could require more runs. Moreover, previous studies have demonstrated that applying softer scaling factors (that is, closer to 1) can improve the discrimination between different ligands.22 This would mean an increased computational burden. On the other hand, since faster solutions are required during the optimization of novel drug-like molecules, the qualitative interpretation of the unbinding simulations can be used instead. Figure 3 shows the distribution of the computational unbinding 11



Page 12 of 44

times for each ligand. In terms of prioritizing ligands, a valuable indication is gained by focusing on the separation between faster and slower unbinding events. Accordingly, during a drug optimization campaign, 9 would have been discarded and 15 preserved. However, 100 would have arisen as a false positive. In this sense, as indicated by the experimental data, 100 displays a kon one order of magnitude smaller and a Kd much larger than that of the other ligands, suggesting that 100 falls into a different binding kinetic regime. This can be attributed to the alteration of the solvation/desolvation enthalpy of the un/binding processes due to a change in size and/or solvability (i.e. the aromatic rings are more hydrophobic in 100 than in the other compounds) within a series of congeneric ligands of the same protein. As noted above, a peculiar feature of hDAAO is the presence of a loop at the entrance of the binding site, including residues 216 to 228,25,37 referred to as the active site lid.32 Hypotheses have been previously advanced that this loop possibly plays a role in ligand binding.32 Indeed, conformational rearrangements of flexible loops located in the vicinity of ligand binding sites in response to local conditions are common to other proteins. One example is provided by βsecretase (BACE1), comprising a β-hairpin loop, known as the flap, which controls accessibility to the binding pocket. Due to protonation state variations of protein titratable residues facing the flap in response to pH changes, multiple conformational states, associated with different binding site accessibility, can be visited by the loop.40–42 Most of the hDAAO ligands possess an aromatic, double-ring moiety that tends to dock in the active site, stacked between the isoalloxazine ring of FAD and the side chain phenol of Tyr 224 (Figure 2B).36–38,43 This stacked π-π interaction was observed in the bound state of several crystal complexes.43,44 In this configuration, access to the binding site is hindered (Figure S2). The hDAAO binding pocket has 12


Page 13 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


a relatively small volume, and this feature is confirmed by the size of various known inhibitors. However, we also discovered larger ligands that were able to inhibit the enzyme. The corresponding structures showed an open state of the active site lid.32,43 We monitored the active site lid behavior during ligand release in the scaled MD simulations, since no position restraints were applied to the loop. Transitions along all degrees of freedom of the system are significantly enhanced under scaled MD conditions. However, a reliable trend can be expected, particularly for the four compounds displaying predicted kinetics in agreement with the experimental data. As such, due to the smoothed potential, we were able to observe the loop transitioning between different conformations upon ligand unbinding during the relatively limited duration of our scaled MD simulations. Contrarily, reproducing spontaneous unbinding/binding process and the movement of the active site loop would likely require extremely long plain MD simulations. We monitored the overall stability of the protein and the behavior of the loop by carrying out three plain MD simulations with no ligands in the binding site. In particular, as the active form of the enzyme is a homodimer,25 two 100-ns-long runs for the monomer, and one 100-ns-long run for the dimer were performed. In terms of the Root Mean Square Fluctuation (RMSF) on the α-carbons (Figure S3), higher RMSF were observed for residues 216-228, in addition to high values for residues 295-305 in a fully solvated loop located a significant distance from the binding site. Although no active site lid opening/closure event was observed in the aggregated 400 ns of plain MD simulations for the monomer, this evidence supports the intrinsic flexibility of the active site loop. Focusing on the active site lid in the scaled MD runs during the dissociation process, we noticed two major behaviors, which we refer to as pathway A and pathway B. In pathway A, 13



Page 14 of 44

Tyr224 is fully solvated before the ligand is able to leave the active site. With the loop opened in this configuration, the entrance to the cavity is no longer shielded and it can be accessed by water molecules. In turn, this facilitates the ligand detachment, and unbinding takes place. In contrast, when pathway B is followed, Tyr224 points towards the binding pocket and is bent towards the bottom of the site. While the change with respect to the initial configuration is not as drastic as in pathway A, the cavity is also more accessible as a result. Thus, the ligand is able to squeeze through the available volume that connects with the bulk, once the stable interactions with protein residues are broken. Notably, our simulations started with the loop in the closed conformation, as in the PDB 3CUK that was used a template to reconstruct the starting complexes for our scaled MD runs. To monitor the two possible cases, we calculated the distance between Tyr224 and the binding site centers of mass (COM) as a function of the simulation time. Sample runs with unbinding in pathway A or pathway B are shown in Figure 6. We considered the four ligands whose predicted unbinding times agreed with the experimental koff. For these ligands, we recorded the frequency with which the two possible pathways were followed over the performed scaled MD runs (Figure 7). At first glance, pathway A was clearly preferred in most cases i.e. the opening of the lid was the limiting step necessary to allow the inhibitors to jump into the bulk. Despite the lower frequency with which pathway B was followed, it is interesting that all the ligands were able to access it. One relevant aspect is the different timescales related to the two paths. Notably, when pathway B was followed, dissociation was generally faster. Conversely, pathway A was more likely to be followed, but was slower. Although these observations are in agreement with the involvement of the loop in

14


Page 15 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ligand (un)binding, little can be concluded about the active site lid behavior because insufficient statistics were gathered and because of the significantly smoothed PES due to the scaling factor.

3. CONCLUSION Real-life scenarios in drug discovery and development programs demand practical, effective, and relatively fast procedures that can be applied to a group of small molecules. The established protocols focus on thermodynamic properties to estimate protein-ligand binding affinity. However, it is becoming clear that the inclusion of binding kinetics can be crucial in identifying promising drug-like molecules. Thus, researchers in academia and industry are gradually including protocols to evaluate kinetic parameters during ligand optimization. On the one hand, being able to discriminate ligands based on their kinetic profile in a routine manner would be highly useful in identifying ligands with desired kinetic behaviors. On the other hand, understanding the molecular features underlying the (un)binding mechanism would be helpful in achieving this goal. Here, we sought to address both matters by performing dissociation simulations for a series of ligands complexed to the protein hDAAO. To this end, we used a computational methodology based on scaled MD that was recently developed to prioritize ligands according to the average unbinding times obtained from the simulations. Based on scaled MD, this method relies on simulating several unbinding events for series of small molecules presenting a common scaffold. We applied the procedure to ligands with the same scaffold and, most importantly, subtle modifications in term of substituent groups. Without any a priori knowledge about their threedimensional structure and residence time, these ligands were ranked according to their average 15



Page 16 of 44

theoretical unbinding times. Experimental kinetic parameters, including off rates and inhibition of enzymatic activity, were determined and showed a good correlation with the computational results. Besides being able to distinguish between fast and slow binders, the ranking obtained through the computational procedure was in fair agreement with the experimental values, with the exception of 100. Assessing this consensus is a necessary step towards the goal of using the simulations as a reliable independent tool. Moreover, we investigated the scaled MD trajectories to extract potential information about the molecular mechanism of dissociation. Interestingly, we recognized two major conformations in which the active site lid, located at the entrance of the binding pocket, could be found during the ligand release process. Although no real quantification could be pursued with the current setup, the analysis undoubtedly provided interesting indications about the role of the loop. Thus, these could be further examined in future studies to achieve a more thorough characterization and understanding of the binding release mechanism. While the applied strategy is emerging as a promising tool, at the current status the methodology is not devoid of limitations. First of all, the possibility of handling larger sets of compounds will demonstrate its real potential for drug discovery. Indeed, since multiple runs, typically in the order of tens of nanoseconds, are required for each ligand in order to achieve an estimate of the times of unbinding, the computational costs increases significantly as more compounds are included. This mostly limited the applicability to a small set of ligands. Future developments of current hardware and software architectures will be critical in achieving this goal. Another crucial aspect is the reliability of the initial structure for the unbinding simulations, which may have an influence on the results of the unbinding simulations. While slightly altered configurations might not affect significantly the outcomes, it is reasonable to think that starting 16


Page 17 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from completely incorrect poses would lead to unreliable results. In the best scenario and remaining in a perspective view, the crystal complex of a closely related ligand or possessing the same scaffold would be desirable to reconstruct reasonable binding modes. Differently, when there is no clue about the bound state, dedicated strategies, mainly referring to ligand docking, can be envisaged to select most likely binding poses. A further non-trivial step in the application of the methodology is the definition of the unrestrained region of the protein. While recognizing the binding pocket is a more general problem and represents a pre-requisite, neglecting residues in the vicinity of the cavity that may be play a role during ligand release would most likely affect the results of the simulations. While not necessarily compromising the quality of the outcomes, the selection of λ factor also represents a critical tunable parameter. As already mentioned in the text, and clearly indicated in the previous works introducing the approach,21,45 there is currently no efficient and standard procedure for this step. Since the effect of the scaled potential displays a system-dependent response and is related to the specific driving interactions between the protein and the ligand, the identification of an appropriate λ is still empirical and should be tested before running the entire set of simulations. Finally, despites the scaled MD simulations of unbinding provide an appealing material for adducing mechanistic hypothesis, care should be taken when interpreting these out-of-equilibrium trajectory data. Being the overall PES altered, the ligand might undertake unlikely routes when aggressive scaling factors are applied. In this case, gentler λ factors could also be explored on representative ligands. Once an appropriate space is envisaged on which to attempt a reweighting procedure, researchers can reconstruct the underlying free-energy surface and gain valuable insights into the energetics and kinetics associated with the binding event. Notably, the latter scheme has broader application possibilities 17



Page 18 of 44

than just protein-ligand binding. For instance, it could be used to consider conformational changes and protein-protein association.

4. METHODS 4.1. Enzymatic Assay. The effect of different compounds on the enzymatic activity was evaluated by an Amplex Red-based assay.31 The DAAO activity was estimated by the production of hydrogen peroxide monitored via a coupled enzyme assay and reagent Amplex Red (Thermo Fisher Scientific). hDAAO (prepared as described previously),46 horseradish peroxidase (HRP), FAD, and the considered compound were incubated for 30 min. After this period, D-serine and Amplex Red were added, and the reaction proceeded for an additional 30 min. A fluorescent product caused by hydrogen-peroxide-dependent Amplex Red oxidation during hDAAOcatalyzed substrate turnover was measured in endpoint mode (excitation and emission wavelengths of 530 and 590 nm, respectively). The final concentration for the reactive components were as follows: 50 mM sodium phosphate, pH 7.4, 0.06 mg/mL human serum albumin, 7 nM hDAAO, 0.1 units/mL HRP, 4 µM FAD, 35 µM Amplex Red, 5 mM D-serine, 0.8% (w/v) dimethyl sulfoxide (DMSO), 0-1.25 mM compound inhibitor. All enzymatic assays were conducted at room temperature in 96-well plate format. Data were fit to a standard fourparameter equation to determine curve top, bottom, concentration producing 50% inhibition (the IC50), and Hill slope. 4.2. Determination of equilibrium binding constants. The binding of the different small molecules was investigated by adding increasing concentrations of the ligand to a fixed amount of hDAAO (10 µM) in the presence of 4 µM free FAD in 50 mM potassium phosphate buffer, 18


Page 19 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


pH 7.5, and at 15 ºC. Absorption spectra were collected in the 300-800 nm range using a Jasco V-560 spectrophotometer. The absorbance changes at specific wavelengths were obtained by the difference spectra calculated by using the Jasco software.47 The dissociation constant Kd, herein referred to as static Kd, was determined according to a saturation binding model. 4.3. Determination of kinetic constants. Association and dissociation kinetic constants kon and koff were determined by stopped flow spectrophotometry.48,49 Under the same experimental conditions as indicated in the previous paragraph, hDAAO (1.1 mg/mL) was rapidly mixed in a SFM-300 Biologic apparatus with a similar volume of increasing concentrations of inhibitor. The time course of spectral changes was recorded in the 300-700 nm wavelength range, and the absorbance versus time data sets were extrapolated at the wavelength identified by static titrations (see Table 1). Observed rate constants (kobs) at increasing inhibitor concentrations were determined from the time courses by nonlinear regression using single exponential equations. The rate constants of inhibitor association and dissociation were determined by linear regression of the equation kobs versus inhibitor concentration: kobs = koff + kon X where X is the inhibitor concentration. The ensuing dissociation constant, here referred to as dynamic Kd, was determined from the ratio of the experimental kinetic rates (i. e., koff/kon). 4.4. Virtual screening. The Asinex, Analyticon, and Life Chemicals libraries were taken from the ZINC database and aggregated.50 The resulting set, comprising 1,130,000 compounds, was then refined by seeking and removing multiple occurrences, rhodanines, and species possessing reactive groups, such as Michael acceptors, or more than two chiral centers. All tautomers and 19



Page 20 of 44

stereoisomers were generated, and protonation states at pH 7.0 ± 1 assigned through the Epik software.51 The procedure resulted in a final set including 916,000 chemicals. As for the protein, 15 crystal structures were present in the PDB. In particular, 14 were complexes of hDAAO with known inhibitors (PDB IDs: 2DU8, 2E4A, 2E49, 2E82, 3CUK, 3G3E, 3W4I, 3W4J, 3W4K, 3ZNP, 3ZNN, 3ZNO, 3ZNQ, 4QFD, 4QFC) and one was the holo form of the enzyme (PDB ID: 2E48). Notably, different conformations (or absence) of the active site lid were observed in the crystals. Thus, using the co-crystallized ligand and all the available structures for the protein, a cross-docking procedure was carried out. Taking advantage of ROC (receiver operating characteristic) curves, the protein structure from PDBs 2DU8, 3W4I, 4QFD, and 2E82 were selected. All proteins were prepared through the Protein Preparation Wizard pipeline of Maestro, Schrödinger Suite.52 The Glide software53,54 was used to conduct the virtual screening, applying the SP scoring function and using the OPLS_2005 force field.55 The resulting top 1000 scoring molecules were retained and subjected to visual inspection. In particular, during this procedure, compounds sharing structural similarities with potent inhibitors present in the literature were selected as well as compounds displaying remarkably different structures in order to expand the exploration of the chemical space. Additionally, ligands displaying unrealistic poses, including for instance intra-molecular steric clashes, imperfect networks of hydrogen bonds and mere shape complementarity, were discarded. This allowed the selection of the final 30 compounds. Of these, only 24 were available from vendors. They were purchased for characterization with biological assays. 4.5. MD Simulation Setup. First, the complexes between ligands and hDAAO were constructed. As already mentioned, no crystal structure was available in the PDB for the protein in complex 20


Page 21 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with ligands possessing a common scaffold with the chosen compounds. Therefore, we used a closely related small molecule, which was present in the PDB in complex with hDAAO (PDB code 3CUK, resolution 2.49 Å). The compound’s structure is reported in Figure S3 and compared with the common scaffold found for the 5 compounds. All the chosen ligands maintained the identical moiety, thus we expected the same interactions to take place. Therefore, using 3CUK as a template, we created the complexes for our chemicals by placing them in hDAAO binding site with the Cartesian coordinates of the common moiety perfectly overlapping. The resulting superposition is shown in Figure S4. Protonation states, tautomers and flip in residue sidechains were optimized through the Protein Preparation Wizard workflow of Maestro, Schrödinger Suite.52 Once the protein-ligand complexes were constructed, each system was solvated in a cubic box with TIP3P water molecules.56 To neutralize the overall charge, an adequate number of Na+ ions was added. This resulted in systems comprising a total amount of 20129 waters and 8 Na+ ions. The Amber ff99SB-ILDN57,58 was used for the protein, while the ligands were modeled according to the GAFF,59 following the RESP procedure60 to determine the charges. As already mentioned, substrate oxidation by hDAAO is concomitant with the reduction of a molecule of the FAD cofactor.25 Notably, the cofactor flavin group is located at the very bottom of the substrate binding pocket. In all the crystal complexes reported in the literature, the group is in close vicinity of ligands and it is oriented so as to give rise to a pi-pi stacking interaction with them. Therefore, we needed to include it in our system. Cofactors are usually composite molecules, typically comprising diverse chemical groups with different properties, as is the case for FAD. Moreover, due to their relatively large size, they are characterized by many degrees of freedom. Thus, deriving proper reliable parameters for such 21



Page 22 of 44

chemical species is far from trivial. In the present case, we took parameters for the FAD cofactor from the R.E.DD.B. database, where multiple conformations of each building block in the molecules are considered in a systematic procedure to derive charges.61 According to the classic pipeline for system preparation, we first minimized and subsequently equilibrated each protein-ligand complex. In particular, position restraints of about 2.4 kcal/(mol Å2) were first applied to all system heavy atoms and then to α-carbons and FAD and ligand heavy atoms in two subsequent steepest descent runs, 5000-step-long each. After this, during the equilibration phase, temperature was first increased in the NVT ensemble in three subsequent steps lasting 200 ps each. This was followed by a 400-ps-long run in NPT to relax the volume. Finally, a 500-ps-long run in the NVT ensemble was carried out under the same conditions used for the subsequent production phase. We took advantage of a GROMACS 4.6.1 version62 appropriately modified in house to perform scaled MD. Under scaled MD conditions, the PES of the system is scaled by a factor λ, where 0

Binding residence time through scaled molecular dynamics: a

Recommend Documents