Computer-Aided Lead Optimization : Improved Small-Molecule Inhibitor of the Zinc Endopeptidase of Botulinum Neurotoxin Serotype A

Optimization of a serotype-selective, small-molecule inhibitor of botulinum neurotoxin serotype A (BoNTA) endopeptidase is a formidable challenge because the enzyme-substrate interface is unusually large and the endopeptidase itself is a large, zincbinding protein with a complex fold that is difficult to simulate computationally. We conducted multiple molecular dynamics simulations of the endopeptidase in complex with a previously described inhibitor (Ki app of 762.4 mM) using the cationic dummy atom approach. Based on our computational results, we hypothesized that introducing a hydroxyl group to the inhibitor could improve its potency. Synthesis and testing of the hydroxyl-containing analog as a BoNTA endopeptidase inhibitor showed a twofold improvement in inhibitory potency (Ki app of 3.860.8 mM) with a relatively small increase in molecular weight (16 Da). The results offer an improved template for further optimization of BoNTA endopeptidase inhibitors and demonstrate the effectiveness of the cationic dummy atom approach in the design and optimization of zinc protease inhibitors.


INTRODUCTION
Botulinum neurotoxin serotype A (BoNTA) is a protein produced by the spore-forming anaerobic bacterium, Clostridium botulinum. The neurotoxin consists of a light chain (M r ,50,000) and a heavy chain (M r ,100,000) linked covalently by a disulfide bond. BoNTA poisoning selectively inhibits the release of acetylcholine from presynaptic nerve terminals at neuromuscular junctions, thus causing flaccid paralysis and leading to prolonged mechanical ventilation with serious medical sequelae or death following respiratory arrest [1]. In small doses BoNTA is an effective medical treatment for a variety of cholinergic nerve and muscle dysfunctions [2,3]. In addition to its medical applications, BoNTA has gained notoriety as a bioterror agent [4]. Because there is no antidote to BoNTA, effective small-molecule inhibitors of BoNTA are highly sought after as antidotes and as potential medical tools for modulating clinical uses of the toxin.
The light-chain domain of BoNTA is a zinc endopeptidase that specifically cleaves SNAP-25, a neuronal protein required for acetylcholine release [5]. The endopeptidase has been used as a target for developing small-molecule inhibitors of BoNTA [6][7][8][9][10][11][12][13][14][15]. Recently, we reported the development of a serotypeselective, small-molecule inhibitor of BoNTA with a K i of 1262.6 mM (1, Figure 1) [7]. Developing or optimizing smallmolecule inhibitors of BoNTA endopeptidase is, however, as challenging as developing small-molecule inhibitors of proteinprotein complexes, a challenge that has been known for decades [16]. The challenge in the case of BoNTA endopeptidase inhibitors can be partly attributed to the substrate that wraps around the circumference of BoNTA endopeptidase upon binding and creates a substrate-enzyme interface of 4840 Å 2 [17]. This large interface requires a high-affinity small molecule to block it. The extent of this challenge can be appreciated by comparing the interface area of BoNTA endopeptidase (4840 Å 2 ) to the typical protein-protein interface area of 750-1500 Å 2 [16]. Another part of the challenge in developing and optimizing these inhibitors is the endopeptidase itself, which is a large zinc-binding protein with a complex fold that is capable of large-scale conformational changes in solution and that is therefore difficult to simulate computationally [17]. Consequently, computer-aided optimization of BoNTA endopeptidase inhibitor leads has not been reported until now.
We report computer-aided optimization of 1 to arrive at an improved, serotype-selective, small-molecule BoNTA endopeptidase inhibitor with an ''apparent'' K i (K i app ) of 3.860.8 mM. This optimization was guided by multiple molecular dynamics simulations (MMDSs) of the zinc-containing endopeptidase in complex with 1 using the cationic dummy atom (CaDA) approach. This approach introduces four identical dummy atoms in a tetrahedral geometry around the zinc ion and transfers all the atomic charge of the zinc divalent cation evenly to the dummy atoms. The four peripheral atoms are ''dummy'' in that they interact with other atoms electrostatically but not sterically, thus mimicking the 4s4p 3 vacant orbitals of the zinc divalent cation that accommodate the lone-pair electrons of zinc coordinates [18][19][20][21][22][23]. The results offer an improved template for further optimization of BoNTA endopeptidase inhibitors and demonstrate that the CaDA approach is useful for both design and optimization of zinc protease inhibitors.

Design
Inhibitor 1 was designed to coordinate the zinc divalent cation embedded in the active site of BoNTA endopeptidase for affinity and simultaneously to interact with serotype-specific residues in the active site for selectivity [7]. This design was based on previous MMDSs using the CaDA approach. The MMDSs (20 simulations) of the endopeptidase in complex with 1 showed that (1) the hydroxamate group coordinated the active-site zinc ion; (2) the phenyl group substituted at the thiophene ring had a p-p interaction with Phe193 and a cation-p interaction with Arg362; (3) the indole ring was engaged in a cation-p interaction with Lys165; (4) the phenyl group attached to the indole ring had a van der Waals interaction with the side chain of Leu527 and a cation-p interaction with Lys165; (5) the ammonium group interacted with the carboxylates of Glu54 and Glu55 [7]. The absolute free energy binding between 1 and the endopeptidase was estimated to be 27.5 kcal/mol according to a free energy perturbation calculation of the MMDS-derived model of the 1-bound endopeptidase using a published method [24] with modifications described in MATERIALS AND METHODS. These computational observations were consistent with the experimentally determined K i app of 762.4 mM (DG ranging from -6.9 to -7.3 kcal/mol) for 1 (see below). In this context, the MMDS-derived model of the 1-bound endopeptidase was used for the following lead optimization. An analysis of the MMDS trajectories of the endopeptidase in complex with 1 identified a vacancy around the meta position of the phenyl group substituted on the thiophene of 1. Synthetically, this void can be filled by a hydroxyl group substituted at the phenyl ring. This hydroxyl group can form hydrogen bonds with active-side residues to improve the affinity for the endopeptidase and the introduction of this hydroxyl group can also increase the hydrophilicity of 1 because dimethyl sulfoxide is needed to dissolve 1 in water. These considerations led to the design of inhibitor 2 ( Figure 1). MMDSs (20 simulations) of the endopeptidase in complex with 2 were carried out to confirm the anticipated hydrogen bonds. The result of these simulations suggested that 2 binds at the active site of BoNTA endopeptidase in a manner similar to that of 1 and that the hydroxyl group attached to the phenyl group of 2 indeed has hydrogen bonds with Arg362 and Asp369 of the endopeptidase ( Figure 2). In the average structure of the endopeptidase complex obtained from 10,000 instantaneous structures at 1.0-ps intervals during the last 0.5-ns period of the 20 different simulations using an explicit water model [25], the hydrogen bond of 2 to Arg362 is bridged by a water molecule; the average distances from the phenolic oxygen atom to the carboxylate oxygen atom of Asp369 and the water oxygen atom are 2.9 Å and 2.3 Å , respectively; the average distance between the water oxygen atom and the closest guanidinium nitrogen atom of Arg362 is 3.3 Å .

Synthesis
The initial synthesis of 2 followed a published scheme [7] that was devised to synthesize 1. The starting material methyl 2-(2-(3hydroxyphenyl)thiophen-3-yl)acetate (4) was prepared using Suzuki coupling [26][27][28] (Figure 3). However, the yield of Friedel-Crafts acylation [29,30] for preparing 5 (Figure 3) was reduced to  10%, presumably owing to the hydroxyl group substituted at the phenyl ring. To increase the yield, a new scheme was devised to perform Friedel-Crafts acylation first and then Suzuki coupling ( Figure 4); this scheme enables facile derivatization of the phenyl group substituted at the thiophene ring via a traditional or combinatorial chemistry approach. As shown in Figure 4, Heck alkynylation [31] of 7, which carries both bromo and iodo atoms, was selectively achieved to afford 8 by using PhCCH/ Pd(PPh 3 ) 2 Cl 2 , K 2 CO 3 , and Et 3 N in DMF. A catalytic amount of InBr 3 [32] was used for the indole formation to obtain 9 in a high yield. To obtain 10, N-alkylation of the indole ring was carried out under the same conditions reported for the synthesis of 1 [7]. Refluxing 10 with 3-hydroxyphenylboronic acid in toluene:ethanol containing PdCl 2 (PPh 3 ) 2 and Na 2 CO 3 gave 11 in a higher yield than a known procedure [33]. The final product 2 was obtained by a reaction of 11 with excess hydroxylamine [34] that simultaneously converted methyl ester and phthalimide to hydroxamic acid and amine, respectively [35,36]. High-performance liquid chromatography (HPLC) purification using solvent containing trifluoroacetic acid yielded pure 2 in its trifluoroacetic acid salt form that was used for the following testing.

Testing
HPLC-based kinetics assays [37] were used to measure the inhibition of BoNTA endopeptidase and a related endopeptidase from botulinum neurotoxin serotype B by 1 and 2. A K i value of 1262.6 mM was previously obtained for 1 that was assumed to be a competitive inhibitor [7]. A further investigation on inhibitor effectiveness in this study revealed nonlinear Dixon plots for both 1 and 2 with data obtained at the condition of [S],,K M ( Figure 5). The curvature of these plots indicates a partially competitive inhibition mechanism for 1 and 2 [38]; it precludes determination of K i values using standard kinetic methods and invalidates the previously reported K i value for 1. To compare the relative inhibitory potencies of the two inhibitors, K i app values of 1 and 2 were therefore obtained from the slope of 1/n versus [I] at the conditions of [S],,K M and inhibitor concentrations less than 15.0 mM and 5.4 mM, respectively, according to a literature procedure [38,39]. The Dixon plots for 1 and 2 are approximately linear at relatively low concentrations ( Figure 5). Consequently, the kinetics assays showed that (1)    endopeptidase by 1 or 2 was due to nonspecific zinc chelation. All the observations demonstrate that 2 is a more potent BoNTA endopeptidase inhibitor than 1.

DISCUSSION Improved Inhibitor of BoNTA Endopeptidase
Our results show that 1 and 2 are competitive inhibitors of BoNTA endopeptidase at inhibitor concentrations less than 15.0 mM and 5.4 mM, respectively, and that the K i app values of 1 and 2 are 762.4 mM and 3.860.8 mM, respectively. However, both inhibitors are partially competitive inhibitors at relatively high inhibitor concentrations. These observations are not unprecedented. It has been reported that a Dixon plot of a dipeptide inhibitor (Gly-Trp) of angiotensin converting enzyme showed a competitive inhibition mechanism at relatively low inhibitor concentrations and a complex inhibition mechanism at relatively high concentrations, whereas the same inhibitor showed a competitive inhibition mechanism in a Lineweaver-Burk plot [39]. Interestingly, both BoNTA endopeptidase and angiotensin converting enzyme are zinc-bound enzymes. Further studies are needed to investigate the detailed inhibition mechanisms of 1 and 2 at relatively high concentrations. Nevertheless, the results of this study show clearly that 2 is a more potent BoNTA endopeptidase inhibitor than 1, as evident from the Dixon plots shown in Figure 5 and from the half maximal inhibitory concentrations determined for the inhibition of BoNTA endopeptidase by 1 (1561.5 mM) and 2 (7.960.8 mM).

Efficiency of computer-aided optimization
While the optimization of 1 reported here was guided by MMDSs using the CaDA approach, 12 other analogs of 1 were made according to the scheme shown in Figure 4 using commercially available phenylboronic acids without the MMDSs guidance. Interestingly, none of these brute-force-approach-based derivatives was found to be more potent than 2 (data not shown). In terms of lead optimization for in vitro potency, which is an important step in drug development, the present work offers another demonstration of the relative efficiency of using computer simulations instead of the brute force approach for lead optimization. The results suggest that the CaDA approach is useful for both design and optimization of zinc protease inhibitors. Strategy for endopeptidase inhibitor optimization As described above, the interface between BoNTA endopeptidase and its substrate is more than twice as large as typical proteinprotein interfaces [16,17]. This large interface causes a conflict in the design and optimization of medically useful BoNTA endopeptidase inhibitors in that, on the one hand, a large molecule is needed to gain more interactions at the unusually large interface, and on the other hand, a small molecule is required to ensure good cell permeability. To avoid the conflict, optimization should be directed towards improving inhibitor affinity with a small increase in molecular weight. This strategy is necessary in early stages of optimization of a micromolar lead because a bivalence or fragment-based approach will be used in later stages of the optimization to take advantage of peripheral sites or exosites on the endopeptidase to achieve high affinity [12,13,17,40,41]. The bivalence or fragment-based approach will be inapplicable to the optimization if the molecular weight of a lead is already too high. At present, the molecular weights of 1 and 2 are 524 and 540 Da, respectively. This study shows that it is practical to achieve a twofold improvement in inhibitory potency with a molecular weight increase of only 16 Da by adding one oxygen atom, given the guidance of the MMDSs of BoNTA endopeptidase using the CaDA approach. Further optimization of 2 following the example reported herein is underway and will be reported in due course.

Computational Methods
Preparation of the zinc-and inhibitor-bound endopeptidase The initial structure of the zinc-containing BoNTA endopeptidase was taken from an available crystal structure of a BoNTA endopeptidase mutant (Glu224Gln and Tyr366Phe; Protein Data Bank code: 1XTG, residue numbers of 1XTG described in this section deviate by one from the residue numbers described in RESULTS and in Figure 2) [17]. To change the mutant back to the wildtype, residues 224 and 366 were mutated to Glu and Tyr, respectively, and the dihedral N-CA-CB-CG of Tyr366 was changed from 65.66 to -171.84u of arc. Hydrogen atoms of the endopeptidase were added by using the QUANTA97 program (Accelrys Software, Inc, San Diego, California). The zinc divalent cation in the crystal structure was replaced by the tetrahedron-shaped zinc divalent cation that has four cationic dummy atoms surrounding the central zinc ion [19] (for more information on zinc protein simulations see ''http:// mayoresearch.mayo.edu/mayo/research/camdl/zinc_protein. cfm''). His223, His227, and Glu262, which are the first-shell zinc coordinates, were treated as histidinate (HIN) and glutamate (GLU), respectively [19,20,[42][43][44]. Glu261 and Glu351, which form a hydrogen bond with His223 and His227 respectively, were treated as glutamic acid (GLH) [19,20,[42][43][44]. His170 and His269 were treated as the histidine whose epsilon nitrogen atom is protonated (HIE); His39 and His230 were treated as the histidine whose delta nitrogen atom is protonated (HID) and histidinium (HIP), respectively. Inhibitor 1 or 2 was manually docked into the active site of the endopeptidase according to a previously reported 3D model of inhibitor 1-bound BoNTA endopeptidase [7]. One chloride ion was added next to the ammonium group of Lys381 located on the surface of the protein to neutralize the protein. The published force field parameters of the tetrahedron-shaped zinc divalent cation [7] were used in energy minimizations and in the following described molecular dynamics simulations. The RESP charges and the AMBER force field parameters of inhibitors 1 and 2 were generated by using the ANTECHAMBER module of the AMBER 7 program [45] and the structures of 1 and 2 that were optimized at the HF/6-31G* level with the Gaussian 98 Program [46] (Figure S1 and Tables S1 and S2). A 10,000-step energy minimization was first performed on inhibitor 2 with a positional constraint applied to the rest of the complex. A 50-step energy minimization was then performed on the tetrahedron-shaped zinc divalent cation and 2 with a positional constraint applied to the endopeptidase only. A 200-step energy minimization was lastly performed on the entire complex without any constraint.
Multiple molecular dynamics simulations All MMDSs (2.0 ns for each simulation with a unique seed number for starting velocities of the system) were performed according to a published protocol [47] using the PMEMD module of the AMBER 8 program [45] with the Cornell et al. force field (parm99.dat) [48]. The topology and coordinate files used for the MMDSs were generated by the PREP, LINK, EDIT, and PARM modules of the AMBER 5 program [45]. All simulations used (1)  The solvated endopeptidase complex system was first energyminimized for 200 steps to remove close van der Waals contacts in the solvated system, slowly heated to 300 K (10 K/ps), and then equilibrated for 1.5 ns. The CARNAL module was used for geometric analysis and for obtaining the time-average structure. All MMDSs were performed on 40 Apple G5 processors. The energy minimization of the time-average structure used ''NCYC = 50'' and other default inputs of the SANDER module of AMBER 5 [45], respectively.
Free energy perturbation calculation [24] The absolute free energy of binding between 1 and BoNTA endopeptidase was obtained by perturbing 1 to nothing over 24 windows of perturbation using the thermodynamics integration method implemented in the SANDER module of AMBER 8 [45]. Each window was equilibrated and sampled with MMDSs (10 simulations) using the procedure described above. Each of the MMDSs used a 1.0-fs time step and lasted 0.5 ns and 1.0 ns for equilibration and sampling, respectively. The perturbation used 12/24-point Gaussian quadrature and the mixing rule shown in Equation 1. This exceptionally computing intensive calculation was performed on 480 Apple G5 processors for two months. To a solution of the above ester (10.00 g, 64.02 mmol) in THF (100 mL) was added NBS (11.40 g, 64.02 mmol). The resulting mixture was refluxed for two and a half hours. The solvent was removed in vacuo. MPLC purification (Hex:EtOAc/9:1) of the residue gave (2-bromothiophen-3-yl)acetic acid as a colorless oil (14.55 g, 97%). 1  To a stirred solution of methyl (2-bromothiophen-3-yl)acetate (120 mg, 0.51 mmol) and 4-iodo-3-nitrobenzoyl chloride (150 mg, 0.48 mmol) in anhydrous CH 2 Cl 2 (10 mL) was added AlCl 3 (260 mg, 1.95 mmol) in four portions at 10-minute intervals at room temperature. The resulting mixture was stirred overnight. The reaction mixture was slowly poured onto 5 g of ice and allowed to warm to room temperature. The aqueous phase was extracted with CH 2 Cl 2 (3615 mL). The combined organic layer was dried over MgSO 4 , filtered, and then concentrated in vacuo. MPLC purification (Hex:EtOAc/5:1) of the residue gave 7 as a light yellow solid (158 mg, 61%). 1  Methyl (2-bromo-5-(2-phenyl-1H-indole-6-carbonyl)thiophen-3-yl)acetate (9) To a solution of 8 (25 mg, 0.05 mmol) in EtOAc (5 mL) was added stannous chloride dihydrate (58 mg, 0.26 mmol). The resulting mixture was refluxed for one hour under N 2 . The reaction mixture was poured onto ice (5 g), and basified with saturated NaHCO 3 solution to pH 8. The white milky mixture was filtered through a Celite pad to remove tin oxides. The organic layer from the filtrate was dried over MgSO 4 , filtered, and then concentrated in vacuo. MPLC purification (Hex:EtOAc/4:1) of the crude product gave the desired intermediate methyl (2-bromo-5-(3-amino-4-(phenylethynyl)benzoyl)thiophen-3-yl)acetate as a yellow foam (21 mg, 90%). 1  To a 250 mL flask containing freshly distilled toluene (50 mL) were added methyl (2-bromo-5-(3-amino-4-(phenylethynyl)benzoyl)thiophen-3-yl)acetate (2.93 g, 6.45 mmol) and indium tribromide (1.14 g, 3.22 mmol) under N 2 . The resulting mixture was refluxed for one hour. The solvent was removed in vacuo. MPLC purification (Hex:EtOAc/4:1) of the crude product gave 9 as a yellow solid (2.50 g, 85%). 1 To a stirred solution of 9 (45 mg, 0.10 mmol) in anhydrous CH 3 CN (5 mL) was added CsF-Celite (125 mg) under N 2 , followed by adding N-(4-bromobutyl)-phthalimide (28 mg, 0.10 mmol). The resulting mixture was refluxed for five hours, cooled to room temperature, and filtered. The filtrate was evaporated in vacuo. MPLC purification (Hex:EtOAc/4:1) of the residue gave 10 as a yellow solid (28 mg, 43%). 1 (5 mL 2-(5-1-(4-Aminobutyl)-2-phenyl-1H-indole-6-carbonyl)-2-(3-hydroxyphenyl)-thiophen-3-yl)-N-hydroxyacetamide trifluoroacetic acid salt (2NTFA) To a stirred solution of 11 (50 mg, 0.07 mmol) in THF/MeOH (1 mL each), 1 mL of 50% aqueous NH 2 OH was added, followed by adding a catalytic amount of KCN (CAUTION: KCN is highly toxic and must be handled with extreme care by trained personnel). The resulting mixture was stirred overnight at room temperature. After the solvent was removed in vacuo, the residue was washed with water (365 mL). HPLC purification of the residue gave 2NTFA as a yellow powder (32 mg, 65%). HPLC purification condition: Phenomenex Gemini 5 mm, C18, 4.66250 mm, eluting with linear gradient of 20% of solution A (1000 mL of H 2 O and 45 mL of TFA) to 100% of solution B (100 mL of H 2 O, 900 mL of CH 3 CN, and 45 mL of TFA) over 20 minutes, and flow rate of 1.5 mL/min with a retention time of 12.20 minutes for 2NTFA (see Figure S2 for chromatograms of 2 before and after the HPLC purification). 1   Botulinum neurotoxin inhibition assays Assays of protease activities of BoNTA and botulinum neurotoxin serotype B (BoNTB) were done at 37uC and contained 0.5 mM substrate, 0.5-1.5 mg/mL recombinant BoNTA/BoNTB endopeptidase, 40 mM HEPES, and 0.05% tween at pH 7.3. BoNTA endopeptidase assays also contained 1 mM dithiothreitol, 25 mM ZnCl 2 , and 0.5 mg/mL bovine serum albumin, while BoNTB endopeptidase assays were supplemented with 1 mM dithiothreitol only. Substrate for BoNTA was a peptide containing residues 187-203 of SNAP-25 [37], while that for BoNTB was residues 60-94 of VAMP [51]. Inhibitors were dissolved in dimethyl sulfoxide at 10 times the final assay concentration, then diluted into the assay mixture containing substrate, followed by addition of endopeptidase (i.e., inhibitor and endopeptidase were not preincubated). Assay times and endopeptidase concentrations were adjusted so that less than 10% of the substrate was hydrolyzed. Assays were stopped by acidification with trifluoroacetic acid and analyzed by reverse-phase HPLC as described [37]. Inhibition of BoNTA endopeptidase activity by 1 or 2 was determined in three independent experiments using eight and nine concentrations of 1 and 2 in each, respectively.

Determination of K i app
For both inhibitors, Dixon plots were nonlinear when the entire range of tested concentrations was included, suggesting a partially competitive inhibition mechanism. Nonetheless, replots of the linear portion of the data, where [S],,K M and inhibitor concentrations of 1 and 2 were less than 15.0 mM and 5.4 mM respectively, permitted calculations of K i app values for 1 and 2 using the previously established K M for the substrate measured in the absence of inhibitor [52] according to a literature procedure [38,39].