Rehabilitating drug-induced long-QT promoters: In-silico design of hERG-neutral cisapride analogues with retained pharmacological activity

Abstract

Background

The human ether-a-go-go related gene 1 (hERG1), which codes for a potassium ion channel, is a key element in the cardiac delayed rectified potassium current, I_Kr, and plays an important role in the normal repolarization of the heart’s action potential. Many approved drugs have been withdrawn from the market due to their prolongation of the QT interval. Most of these drugs have high potencies for their principal targets and are often irreplaceable, thus “rehabilitation” studies for decreasing their high hERG1 blocking affinities, while keeping them active at the binding sites of their targets, have been proposed to enable these drugs to re-enter the market.

Methods

In this proof-of-principle study, we focus on cisapride, a gastroprokinetic agent withdrawn from the market due to its high hERG1 blocking affinity. Here we tested an a priori strategy to predict a compound’s cardiotoxicity using de novo drug design with molecular docking and Molecular Dynamics (MD) simulations to generate a strategy for the rehabilitation of cisapride.

Results

We focused on two key receptors, a target interaction with the (adenosine) receptor and an off-target interaction with hERG1 channels. An analysis of the fragment interactions of cisapride at human A_2A adenosine receptors and hERG1 central cavities helped us to identify the key chemical groups responsible for the drug activity and hERG1 blockade. A set of cisapride derivatives with reduced cardiotoxicity was then proposed using an in-silico two-tier approach. This set was compared against a large dataset of commercially available cisapride analogs and derivatives.

Conclusions

An interaction decomposition of cisapride and cisapride derivatives allowed for the identification of key active scaffolds and functional groups that may be responsible for the unwanted blockade of hERG1.

Background

Several classes of potassium channels are involved in regulating the heart rate by setting the amplitude and duration of the action potential and the resting membrane potential. Abnormalities in the function of these ion channels due to inherited mutations or pharmacological blockage can prolong the duration of the action potential, leading to the development of severe arrhythmias (i.e., long QT syndromes - LQTS). Genetic analysis has revealed that mutations in potassium channels, such as the human ether-a-go-go related gene (hERG) and KvLQT1, establish a molecular basis for LQTS[1–4]. A growing number of diseases have also been linked to genetic mutations in potassium channels. The channelopathies related to potassium channels include various cancer types, type 2 Bartter’s syndrome, type 1 episodic ataxia, and hyper-insulinemic hypoglycemia[4]. The best-known feature of the hERG1 channel is its unique promiscuity in binding to a wide range of organic molecules. A broad panel of organic compounds used in common cardiac and non-cardiac medications (e.g., antibiotics, antihistamines and antibacterial agents) are thought to cause a reduction in the repolarizing current I_Kr by blocking the central cavity of hERG and similar channels, leading to ventricular arrhythmia[5]. Several Food and Drug Administration (FDA) approved drugs (i.e., terfenadine, cisapride, astemizole and grepafloxin) have been withdrawn from the market, while others like thioridazine, haloperidol, sertindole, and pimozide have been restricted due to their effect on the function of the hERG channel. The discovery of drug-related arrhythmias has led to mandatory drug screening for hERG1 blockage by both the FDA and the European Medicines Agency (EMEA). Because most of these drugs have high binding affinity profiles for their principal targets, “rehabilitation” studies aimed at reducing their side effects (i.e., decreasing their high hERG1 blocking affinities) while keeping them efficient at binding to their original targets have become increasingly common. These studies may allow these drugs to reenter the market.

We chose cisapride to study as a model hERG blocker that can potentially be rehabilitated and returned back to the market. Researchers at Janssen Pharmaceutica discovered cisapride (a gastroprokinetic agent; trade names: prepulsid, propulsid) in 1980. Cisapride, a prokinetic agent, increases gastrointestinal motility and acts as a selective serotonin agonist for 5HT-4 receptors. It also relieves gastrointestinal symptoms (i.e., constipation and bloating) by indirectly stimulating the release of acetylcholine in muscarinic receptors. In the gastrointestinal tract, the activation of these receptors stimulates smooth muscle contraction. While cisapride’s main target is thought to be the 5-hydroxytryptamine receptor 4 (5HT-4), it can selectively target a number of other G-protein-coupled receptors (GPCRs). For example, cisapride has binding affinities for 5HT-4 (~14 nM) and adrenergic receptors (16 nM at the α-1 adrenergic receptor). Cisapride was marketed in the USA from 1993 to 2000; the use of cisapride has been associated with 341 reports of cardiovascular problems, including LQTS and torsades-de-pointes, and 80 reports of death[6]. After seven years on the market, cisapride was withdrawn in the USA and was limited for use in many other countries due to its high hERG1 blocking affinity (IC₅₀ 6.5 nM, Mohammad et al.[7]; IC₅₀ 44.5 nM, Rampe et al.[8, 9]; IC₅₀ 15.0 nM Drolet et al.[10]). A goal of our proof-of-principle study is to establish the feasibility of an in silico cardiotoxicity assessment with a multi-target computational approach for rehabilitation that:

• (i) Assesses drugs for their hERG-blocking ability;

• (ii) Identify the active components responsible for the original target activity;

• (iii) Caps or modifies the moieties responsible for the hERG blockade.

The workflow chart is shown in Figure 1. We combined step-by-step ligand modification focusing on key functional groups (Figure 2), two-target receptor docking, and molecular dynamics (MD) simulations to achieve this goal. First, we determined the key-molecular fragments of cisapride responsible for its high-affinity binding to the A_2A receptor using all-atom MD simulations and Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) binding energy decompositions, which are similar to the approaches used in previous studies[11–13]. Next, small modifications of the sites responsible for the hERG blockade that played a less significant role in the stabilization of cisapride in the binding pocket of the A_2A receptor were made. Our previously developed atomistic models of the open-state hERG1 channel were used to guide the modifications of cisapride. The most potent drug variants were tested for their binding to the central cavity of the hERG1 channel. The compounds that showed low affinities were then selected for further analysis. Because no crystal structure of the 5HT-4 receptor is currently available, we used a recently crystallized agonist-bound active GPCR from the same family as the target (the human adenosine A_2A receptor, PDB ID: 3QAK)[14]. Since cisapride can selectively bind to a number of other GPCRs (including A_2A receptor) with a similar range of binding affinity at a similar binding pocket, we chose to work with the available agonist-bound active target site of the A_2A receptor to screen cisapride and its derivatives. Prior binding site[15] and comparative homology modeling studies for these receptors indicate a high similarity between the 5HT-4 and A_2A binding pockets. Finally, a list of possible modifications to the original cisapride molecule was generated. MD simulations with MM/GB-SA computations, database drug screening and de novo design studies clearly showed that the shorter alkyl chain in cisapride analogues are key to retaining their binding to the A_2A receptor while remediating the blockade of the hERG1 channel. To compare in silico results to de novo developments we screened a large panel of already synthesized cisapride analogs. Small molecule databanks (i.e., ZINC[16]) were screened for synthesized cisapride derivatives and the relevant literature was reviewed to identify their activity in both the GPCR and hERG targets. Dual target docking combined with all-atom MD simulations were used to establish the key interactions between cisapride and its derivatives responsible for their high-affinity binding to the A_2A and how they receptor triggering hERG1 blockade. Using this combination of techniques it is possible to assess novel compounds a priori for their cardiotoxicity risks associated with hERG1 blockade as well as to identify sets of functional groups responsible for on-target binding.

Figure 1

The block-scheme for the work-flow in the computational rehabilitation of cisapride derivatives.

Figure 2

The group definition of cisapride (i.e., R1 = -C ₆ H ₄ -F; R2 = -O-CH ₂ -; R3 = -CH ₂ -CH ₂ -; R4 = -C ₅ H ₃ N-O-CH ₃ ; R5 = -NH-CO-; R6 = C ₆ H ₂ -Cl-NH ₂ ; R7 = -O-CH ₃ ) for per-group decomposition and inter-monomer contacts from the MM/GBSA analysis.

Methods

Molecular docking

The structures were initially optimized using Schrodinger’s Macromodel module to perform an all-atom MM geometry optimization. The optimization was carried out with the OPLS 2005 force field and the Polak-Rebiere Conjugate Gradient (PRCG) energy minimization method with a 0.001 kcal mol⁻¹ Å⁻¹ energy gradient convergence criterion[17]. The resultant ligand structures were docked to the targets using the following docking algorithms: Glide/Induced Fit Docking (IFD)[18], FlexX[19], Autodock[20], and Generalized Optimized Ligand Docking (GOLD)[21]. The details of the docking algorithms are described below:

Glide/IFD: The Glide-XP (extra precision) (v.5.0) and Induced Fit Docking (IFD) modules of the Maestro suite were used for the docking calculations. The docking studies were performed with the following steps: (i) constrained minimization of the receptor with an root mean square deviation (RMSD) with a cutoff of 0.18 Å; (ii) initial Glide docking of each ligand using soft potentials; (iii) refinement of the derived docking poses (i.e., minimization of the docking poses within 20 Å of the ligand poses) with Schrodinger’s Prime module; and (iv) Glide re-docking of the protein-ligand complexes. GOLD: The GOLD program (v.5.0.1) was used with two default docking scores (GOLD Fitness and ChemScore). Part of the receptor's flexibility was accounted for by assigning flexibility to 10 specific amino acid residues at the active site. The side chains of these amino acid residues were selected as flexible rotamers. The rotamers progressed by 10° increments to cover a full 360° rotation. The default genetic algorithm parameters (100 for the population size, 5 for the number of islands, 100000 for the number of genetic operations and 2 for the niche size) were used. However, the maximum number of runs was set to 100 for each docking simulation. FlexX: The FlexX program (v.4.0) from BioSolveIT was also used. The default algorithm parameters were used for the docking and construction of the active sites of the receptor. The solutions per fragment and per iteration were both set to 2000. AutoDock (v.4.0): The number of grid points in each direction was 126 with a grid spacing of 0.4 Å. The number of hybrid Genetic Algorithm-Local Search (GA-LS) runs was 200.

Homology modeling

The SWISS-MODEL homology modeling program[22] was used for the development of the 5HT-4 receptor model. A multiple sequence alignment was performed using the CLUSTALW algorithm[23]. A β1 adrenergic receptor (A_2A with a carvedilol agonist, PDB ID: 4AMJ) was used as the template because it had the highest sequence identity percentage in the sequence alignment (41%). Protein models were generated from the alignment in a stepwise manner. The backbone coordinates for the aligned positions were extracted from the template and the regions of insertions/deletions in the alignment were found by searching either a loop library or a conformational space search using constraint space programming. The templates were weighted by their sequence similarity to the target sequence and outlier atomic positions were excluded. The scoring function used for assessing favorable interactions (hydrogen bonds, disulfide bridges) and unfavorable close contacts for determining side chain conformations was derived from a backbone-dependent rotamer library[24].

MD simulations

All-atom MD simulations were carried out using CHARMM (v. c36a2)[25]. All simulations were carried out at 323 K and 1 atm using periodic boundary conditions (PBC) with the NPT ensemble. The particle mesh Ewald (PME) algorithm was used for the long-range electrostatic interactions. Both the A_2A receptor from PDB coordinates[14] and a model of the hERG1 channel[1, 2, 26] were embedded into the DPPC membrane bilayer using the CHARMM-GUI membrane builder protocol[27]. Structures were minimized and equilibrated with gradually decreasing harmonic constraints (i.e., they were initialized with 10.0 and 5.0 kcal mol⁻¹ Å⁻² for the backbone and side chains, respectively, and gradually decreased to 0.5 and 0.1 kcal mol⁻¹ Å⁻², respectively) over 2 ns (for equilibration) and then subjected to a 50 ns production run.

MM/GBSA computations

The enthalpies of cisapride binding to both targets were computed as averages from an ensemble of structures (5000 for each binding pose) sampled from evenly distributed points over five independent 50 ns runs. The electrostatic contributions to the desolvation components of the binding enthalpies were obtained by using the Generalized Born Solvation Energy Module and an Implicit Membrane (GBIM) module as implemented in CHARMM[28, 29]. A dielectric constant of 2 was assigned to the protein, and the protein-solvent surface was defined using a set of optimized atomic radii from Nina et al.[30]. Following the numerical recipe of Chandra et al.[31], the membrane was represented as a 24 Å slab with a dielectric constant of 2. To better describe the lipid dynamics, we extracted the positions of the lipid’s heavy atoms in the bilayer and used them to represent the neutral slab around the protein and cisapride. The protein-ligand interaction components of the binding enthalpies were computed following a standard protocol with infinite cut-offs[12, 32].

Results and discussions

MD simulations

As stated in the Introduction, one of the major goals of this study was to identify the determinants of the cisapride blockade of hERG and its binding to the A_2A receptor which serves as a suitable model for the 5HT-4 receptor. Although our docking studies unambiguously identified binding pockets in both proteins, the orientation of the bound ligand was less well defined. To account for the limited resolution of our docking studies, we ran five independent 50 ns full-membrane all-atom molecular dynamics (MD) simulations using the top 5 poses from docking analysis as a starting point. All simulations were found to produce stable trajectories with membrane proteins displaying heavy atom RMSD values of ~3 to 4 Å, which is comparable to root mean squared fluctuation (RMSF) estimates from previous studies of K-channels with available crystal structures[33]. The averaged locations of the ligand after MD simulations are shown in Figures 3 and4. Cisapride bound in either membrane protein (A_2A receptor or hERG1 open-state channel model) displays significant conformational dynamics in the binding site with an average RMSD value of ~1.1 Å relative to the average structure. The backbone RMSD of the hERG tetramer channel after 50 ns was approximately 5.3 Å including deviations from the symmetric tetramer; the corresponding values for the monomers were from 3.5 to 4.2 Å. The positional fluctuations in backbone atoms of residues forming the pore domain (PD) of the hERG1 monomer plateaued at a RMSD of ~2.4-2.9 Å. These values are similar to those reported previously for MD simulations of Kv channels[32]. The backbone RMSDs of the A_2A receptor plateaued at ~3.1 Å. The considerable dynamics of the receptor and the bound ligand suggest that explicitly accounting for site and ligand flexibility in evaluation of the binding energies is necessary. To circumvent obvious limitations of the chosen docking strategy and to identify key interacting partners e.g. relevant amino-acid residues and key functional groups in the drug, we performed MM/GBSA computations to estimate the binding affinity of cisapride in the A_2A receptor and the open-state hERG1 model. The distribution of binding enthalpies from 5 independent simulations is shown in the Figure 5. The average binding enthalpies for cisapride to the hERG1 and A_2A are –21.3 ± 2.8 kcal/mol and -24.3 ± 1.9 kcal/mol, respectively.

Figure 3

(left) Dynamics of cisapride in the hERG1 cavity (Backbone RMSD after 50 ns ~5.0 Å for a tetramer, 4.2 Å for a monomer; the pore domain has an RMSD of ~2.9 Å). (right) Summary of the MM/GB-SA decomposition analysis (per-residue and per-group contributions) from the last 40 ns of the MD simulation of the hERG1-cisapride complex. (see Figure 2 for a definition of the R groups at cisapride).

Figure 4

(left) The dynamics of cisapride in the A _2A binding pocket (Backbone RMSD after 50 ns ~3.1 Å). (right) Summary of the MM/GB-SA decomposition analysis (per-residue and per-group contributions) from the last 40 ns of the MD simulation of the A_2A-cisapride complex.

Figure 5

Normalized probability distribution function of the binding enthalpies for cisapride to WT-A2A and WT-hERG1 proteins.

Analysis of the per-residue and per-functional group contribution to binding enthalpies allows us to determine residues and drug fragments, which uniquely determines the binding mode of cisapride in both systems. (Figures 3 and4) The key residues that determine the thermodynamics of cisapride stabilization in the hERG1 cavity are S624 from the pore helix, Y652 and F656 from the distal S6 and, surprisingly, W568 from the S5 helix. The high-affinity intra-cavity blockade of hERG1 by cisapride depends on the hydrogen bonds formed by the pore-helix and on stabilizing aromatic interactions with residues in the S5 and S6 helices. The drug dimensions allow for inter-subunit interactions (see Figure 2) with matching hydrophobic moieties in the hERG1 cavity while the carbonyl backbone of cisapride interacts strongly with T623 and S624. The binding mode is similar to that for dofetilide and other molecules known to restrict hERG1 K⁺ currents by binding to the protein’s internal cavity[3, 26]. Most of the residues were mapped in previous experimental studies and are thought to be responsible for high-affinity blockades[34–38]. Importantly for both systems, docking studies were able to identify key amino acids and functional groups and thus results obtained from average protein structure and flexible ligand docking appear to be sufficiently predictive to allow for micromanagement of the drug.

The decomposition of binding enthalpies for cisapride binding to the A_2A adenosine receptor is illustrated in Figure 4. While there is a cluster of acidic residues involved in drug binding (E13, E169 and D170), the desolvation penalty for the negative charges located on the surface of the receptor counter balances favorable charge interactions and contributes unfavorably overall to drug stabilization in the pocket. Instead, the bound cisapride is stabilized to a large extent by a network of aromatic and amphipathic groups lining the binding pocket. The constellation of residues includes Y271, L267, I252, T256 and M270 that interact with the linker and the rings (R1, R3, and R6). (See Figure X for labeling of the regions of cisapride) The binding of cisapride is further stabilized by interactions between the residues S67, T68 and F168, and the moieties of the R2 and R3 regions. It is worth mentioning that all of these residues are involved in the formation of binding pockets for adenosine and other agonists, according to previous studies and recently published crystal structures[14]. F168 forms binding interactions with the central methoxy groups of the bound molecule, and other key residues (S67, T68 and Y271) interact with amphipathic functional groups on aromatic rings. The molecular fragments of the drug that form stable interactions with the residues in the A_2A receptor are as follows: a methoxy group on the central heterocyclic ring (R4), a terminal aromatic ring (R1), and a methoxy group (R7) on the benzamide ring. Interestingly, it has been shown that modification of the benzamide ring is essential for modulation of the activity of cisapride derivatives[39].

To better understand the effects of each fragment of cisapride on the A_2A adenosine receptor binding and blocking affinities of hERG1, we also performed per-group decomposition of the binding enthalpies for drug binding to hERG1. The cisapride fragments were defined according to the topology definition of the CHARMM general force field[25, 40] and are illustrated in Figure 2. The resulting decomposition is shown in Figures 3 and4, allowing us to focus on targeted modifications of the fragments. It is apparent that to inhibit binding to hERG1, one may need to modify one or both of the terminal rings (R1, R6) as well as the flexible linkers (R2, R3). Cisapride bound to the hERG1 cavity interacts with a ring of hydrophobic and aromatic residues (Y652, A653 and F656) and the Ser/Thr-rich apex of the pore helix (T623, S624)[2]. Experimental studies have also demonstrated that the matching interactions between cisapride and aromatic rings (F656 and Y652) are essential for the hERG blockade[41]. The decomposition analysis shows that one of the key interactions important for the hERG blockade but unrelated to the high-affinity binding of cisapride to the A_2A receptor is the van-der-Waals interactions between the long alkyl chain linker of the drug and a number of hydrophobic residues in the hERG cavity. The docking/MD simulations strategy adopted in this study appears to have yielded a potential route to rehabilitate cisapride. In our next step, we attempt to remove fragments responsible for drug stabilization in the inner-cavity of hERG1 channel, while retaining its binding to the A_2A receptor. Several parallel approaches were considered for de novo development of non-blocking cisapride collected in Table 1:

Table 1 Selected cisapride analogues from the de Novo drug design study

• (i) In the first approach, we simply targeted heteroatoms in the functional groups of cisapride molecule. They were exchanged with hydrogen (-H) to observe the effect of each of these fragments on the docking scores for both receptors.

• (ii) Next, the number of –CH₂ groups of the linker was modified to produce truncated and extended versions of the linker. This affects the distance between the terminal aromatic ring and the six-membered heterocyclic ring piperidine, thus affecting torsional energy. The length of the alkyl chain also modulates the interactions between the aromatic moieties of the drug and the Y652/F656 residues that are key components of the high-affinity blockade of hERG1.

• (iii) In the last approach, the functional groups were removed based on the per-fragment decomposition, as shown in Figures  1 and 3.

A library of 10,700 molecules from the small organic molecules fragment database of Maestro was used to create new analogues of cisapride guided by energy decomposition analysis as well as the per-residue interaction energy analysis described above. The summary of the modification strategy is shown in Additional file1: Tables S1 to S2. All of the developed molecules were initially docked using virtual high throughput screening (VHTS) followed by Glide extra precision (XP) docking for the derivatives showing low hERG blocking binding scores and high binding scores for the A_2A adenosine receptor (~500 compounds).

General note on molecular docking and De novo drug design

All of the derived cisapride analogues were docked to the A_2A (PDB ID: 3QAK) active site using the approach described in the methods section. They were also docked to the pore domain of hERG1 model, and the cisapride derivative results were compared to the original cisapride docking scores. (Figures 6 and7, and Additional file1: Table S1) Four different molecular docking programs (AutoDock, Glide, FlexX, and GOLD) were used to study the binding interactions of cisapride and the cisapride derivatives in the A_2A adenosine receptor and the hERG1 central cavity. Since the FlexX and AutoDock docking programs both underestimated the interaction energies of cisapride and its derivatives at both targets (i.e., hERG and A_2A), GOLD and Glide/XP are generally preferred for the evaluation of binding affinities of known compounds. Our previously developed model of the trans-membrane domains, S1-S6, of open-state hERG1[1, 26] was used to model the binding interactions of cisapride and its derivatives in the pore domain of the channel. Initially, cisapride was docked to A_2A receptor as well as hERG1 model, and these binding scores were used as threshold values for all four docking programs. Cisapride derivatives that displayed similar/higher binding scores (absolute values) in the A_2A receptor and lower binding scores in the hERG1 pore domain were considered for further analysis. For example, one of them (#11) labeled as Cisapride-D11 in the text had docking scores similar to original cisapride in the A_2A binding pocket and considerably lower docking scores in the hERG PD, thus, it was considered for further rehabilitation.

Figure 6

Ligand Binding to hERG1 pore domain. (i) The top-docking pose of cisapride in the hERG1 pore domain. (left) Binding interactions are zoomed and detailed (right). (ii) A 2D-ligand interaction map for cisapride binding to the hERG cavity.

Figure 7

Ligand Binding to A2A receptor. Top panel: The top-docking pose of cisapride in the human A_2A adenosine receptor. (left) Binding interactions are zoomed and detailed (right). Bottom panel: A 2D-ligand interaction map for cisapride binding to the A_2A receptor.

Effect of functional group substitution on drug binding to A_2A receptor and hERG1 channel

All of the candidates with hydrogen substitutions by heteroatoms (i.e., Cisapride-D1 and Cisapride-D8 in the Additional file1: Table S1) had docking scores similar to those of cisapride binding to the A_2A receptor. However all of these derivatives retained high docking scores for binding to the hERG PD and therefore were discarded. Some of the new derivatives, such as Cisapride-D2 (an exchange of a –Cl group with –H), led to a decrease in the docking score in the A_2A receptor (i.e., the Glide/XP docking score decreased from -8.48 to -5.07 kcal/mole for Cisapride-D2), thus these derivatives were also discarded. (Additional file1: Table S1) While the removal of some functional groups from cisapride did not significantly reduce the docking scores of compounds at A_2A, their hERG PD bindings were still too high (i.e., Cisapride-D4, Cisapride-D5).

Varying the number of –CH₂ groups in the linker, R3 between the terminal aromatic ring, R1, and the central heterocyclic ring, R4, affected the docking scores of compounds in the hERG PD. (See Figure2 for labeling) While decreasing the number of -CH₂ groups from three (original cisapride) to one reduced the docking score of the derivative in the hERG1 central cavity, (i.e., for Cisapride-D11, GOLD docking score went from -9.42 to -7.40 kcal/mole; the Glide/XP docking score went from -7.80 to -7.09 kcal/mole) increasing the number of –CH₂ groups from 3 to 4 did not significantly affect the docking score of the derivative in the hERG1 PD (i.e., for Cisapride-D9, the GOLD docking score went from -9.42 to -9.03 kcal/mole; the Glide/XP docking score went from -7.80 to -7.45 kcal/mole) (Additional file1: Table S1). Therefor for further rehabilitation studies, only a decreased number of –CH₂ groups in the cisapride backbone (i.e., Cisapride-D11) was considered. Derivatives of Cisapride-D11 were generated using the Combinatorial Library Enumeration and Screening module of the Maestro molecular modeling package and the corresponding docking results were compared with that of original Cisapride-D11 (Additional file1: Table S2). For this purpose, approximately 600 default fragments in Schrodinger’s Enumerate module were used, and 8600 new cisapride derivatives were generated. After the derivative generation, the compounds were prepared (protonation states were determined at a pH of 7) with the LigPrep module of Maestro and energy minimization was performed for the structures using the PRCG energy minimization method. Structures were then docked onto the A_2A adenosine receptor via the Glide High Throughput Virtual Screening method (Glide/HTVS). This was used to reduce the large number of derivatives to a manageable number for further analysis. The top 100 compounds according to the Glide/HTVS docking scores were used for Glide extra precision (XP) docking both in the A_2A active site and the hERG1 central cavity. The docking results of the selected derivatives of Cisapride-D11 are tabulated in the Additional file1: Table S2. One of the major outcomes of this in silico study is that simple truncation of the flexible linker appears to be sufficient to remediate unwanted off-target interactions of cisapride. Obviously, it is not sufficient just to claim that such modifications will lead to a potent drug; it requires verification. One of the reasons for choosing cisapride is the availability of a large number of previously evaluated cisapride derivatives. Accordingly, to test similarities and differences in predicted cardiotoxicity (hERG1 blockade) and efficacy of the cisapride derivatives we used the ZINC database and literature mining.

Comparison study with a database of available drug databases

94 cisapride derivatives (with >80% structural similarity with cisapride) taken from the ZINC databank were used to obtain the docking scores/poses for binding to the human adenosine A_2A receptor and the hERG1 pore domain using the Glide/XP docking program. The compounds that had docking scores (absolute values) greater than -7.00 kcal/mole in the A_2A receptor binding site (28 cisapride derivatives) were evaluated for their docking scores in the hERG PD (Additional file 1: Tables S3 and S4). Results showed that a few promising compounds (i.e., ZINC05998832, ZINC58529167, ZINC13834042, ZINC20621758, ZINC43023913) with docking scores less than (absolute values) -6.00 kcal/mole in the hERG1 PD, and greater than (absolute values) -7.00 kcal/mole in the binding pocket of the A_2A receptor. (Table  2)

Table 2 Cisapride derivatives from the ZINC database that show low binding scores in the hERG1 PD and high binding scores in the A _2A adenosine receptor

Interestingly, ZINC20621758, also known commercially as mosapride, was detected via blinded dual-target screening to be a safer 5HT-4 agonist than cisapride. Potet et al.[6] showed that mosapride does not block the hERG channel, which is in agreement with lower binding scores predicted in silico. It is also a well-known agonist of the 5HT-4 receptor. The IC₅₀ of mosapride binding to the hERG channel is ~16.5 μM. Mosapride carries a significant decrease of cardiovascular risks related to alterations in QT intervals according to preclinical studies[42, 39, 43]. Carlsson et al. used a rabbit model of the acquired long QT syndrome, and while cisapride prolonged the QT interval, mosapride did not[39, 44, 45]. One of the key differences between mosapride and cisapride is a shorter linker between aromatic rings that, as predicted in the dual-target de novo drug modeling, will have an immediate effect on the drug-hERG PD interactions. Figure 8 shows the top Induced Fit docking pose of mosapride in the A_2A receptor. Cisapride and mosapride have different binding modes in the A_2A receptor due to a different number of CH₂ groups in the linker. Additional file1: Figure S2 shows a superposition of the top-docking poses of cisapride and mosapride.

Figure 8

Mosapride bidning to hERG1 pore domain. Top panel: The top-docking pose of mosapride that shows a low hERG PD blocking affinity in the human A_2A adenosine receptor. (left) Binding interactions are zoomed and detailed (right). Bottom panel: 2D-ligand interaction map for mosapride binding to the A_2A receptor.

The other ZINC compounds in Table 2 also exhibit high docking scores in the A_2A receptor; however there are no experimental results available in the literature for these compounds. However, ZINC05998832 (aka fluoroclebopride) is a close analog of clebopride. The reported IC₅₀ of clebopride in the hERG PD is 0.62 μM. Clebopride is commercially available in Spain and Italy as a gastroprokinetic drug. Although its cardiotoxicity has not been reported in clinical studies, in vitro studies have shown that clebopride prolongs the cardiac action potential duration at 90% (but not 50%) repolarization at 10 μM[44]. Tack et al.[44] also found that no cardiovascular safety concerns were reported for the newer selective 5HT-4 agonists prucalopride, velusetrag, and naronapride or for the non-selective 5HT-4 agonists with no hERG or 5HT-1 affinity, such as renzapride, clebopride, and mosapride. Importantly, the experimental data available for several ZINC compounds indicated nM affinity against 5HT-4 and 5HT-2 receptors, corresponding to high affinity binding for several of the compounds studied (ZINC20621758, Ki ~113 nM; ZINC43023913, 79 nM; ZINC13834042, 2250 nM and ZINC05998832 ~ 283 nM, respectively), all with very low docking scores in the hERG PD. Although a few cisapride analogues (such as ZINC28087327, Ki 40.6 nM in 5HT-4; ZINC43024452, Ki 3.16 nM in 5HT-4) were correctly predicted by in silico screening to have high experimental binding affinities in the 5HT-4 receptor, they also had high docking scores in the hERG PD (Additional file1: Tables S2 and S3).

On the transferability of A_2A predictions to 5HT-4 receptors

There is a well-known challenge in relating the top binding poses produced by docking to the most stable ligand orientations in the binding pocket. Many organic molecules are flexible with a number of rotating bonds. The previous section on mapping the relevant binding sites in the receptor suggests that some of these residues are also responsible for selective substrate and agonist binding to human 5HT-4[15, 46]. It is worth mentioning that the main target of cisapride in gastrointestinal tissues is not the A_2A receptor but the 5HT-4 receptor. To investigate the similarities of the binding pockets of two GPCRs, namely, 5HT-4 and A_2A receptors, we performed homology modeling with subsequent drug docking. The homology modeling was performed using SWISS-MODEL for 5HT-4 (using the highest sequence identity scores of a multiple sequence alignment). β1 adrenergic receptors (with agonist bound) resulted in the highest sequence identity percentage from a sequence alignment (41%) and it was used as a template. Superposition of the 5HT-4 model protein and the A_2A crystal structure reconfirmed that there is high structural similarity between these two GPCRs (RMSD of 2.4 Å) (Additional file1: Figure S1). Furthermore, most of the residues identified by MD simulations as essential for drug stabilization are well-known determinants of high-affinity binding in the 5HT-4 system. We also used the ChEMBL database to search for compounds targeting the 5HT-4 receptor (CHEMBL1875). ChEMBL is a database for bioactive drug-like small molecules that contains abstracted experimental bioactivities from the scientific literature (i.e., binding constants and pharmacology and ADME/Tox data)[34, 35]. In total, 419 associated bioactivity results for the 5HT-4 receptor were found. The results were ranked based on Ki values, and compounds that had better Ki values than 1 nM in the 5HT-4 receptor were kept for consideration. Our aim was to test the docking score accuracies for compounds in the binding pockets of the A_2A and 5HT-4 receptors. The protonation states were determined at pH 7, and structures were optimized using the LigPrep module of Maestro. Subsequently, these compounds were docked to the A_2A and 5HT-4 receptors using Glide/XP Induced Fit Docking. Considering the simplicity of the docking approach, our results show a high correlation between the experimental and predicted binding energies (r² > 0.6 for both compounds in the A_2A and 5HT-4 receptors, docking results, Figure 9). Structures with bioactivity results in the 5HT-4 receptor were also docked in the hERG PD using Induced Fit Docking. Compounds with low predicted binding affinities in the hERG PD are listed in Additional file1: Table S5 and can guide future studies. Our findings reconfirm the similarity of the binding pockets of the A_2A and 5HT-4 receptors as well as the transferability of our results.

Figure 9

Comparison of experimental and calculated (IFD docking scores) binding energies of compounds. Top panel: The A2A receptor and Bottom panel: The 5HT-4 receptor. (Ki values are converted to kcal/mole for comparison, T = 300 K, R = 8.314 J/molK; If a compound has different Ki values for same target, the average of these values is used).

Conclusions

A cross analysis of the interactions between cisapride and its analogues with the human A_2A adenosine receptor and the hERG1 central cavity led us to formulate a computational approach to the rehabilitation of drugs withdrawn from the market due to pro-arrhythmic activity. Our SAR study of cisapride and its derivatives allows for the in silico evaluation of drug potency and the drug’s ability to block hERG. Molecular docking results using both de novo designed compounds and available cisapride derivatives in the ZINC drug database showed that the shorter alkyl chains in the cisapride analogues is key element to retaining their binding to the A_2A receptor and remediating the blockade of the hERG1 channel. We believe that a simple dual-target ‘rehabilitation’ strategy based on an optimization against two protein target structures may be applied to other drugs withdrawn from the market due to their side effects, and may lead to the reuse of these drugs.