Identification of berberine as a direct thrombin inhibitor from traditional Chinese medicine through structural, functional and binding studies

Thrombin acts as a key enzyme in the blood coagulation cascade and represents a potential drug target for the treatment of several cardiovascular diseases. The aim of this study was to identify small-molecule direct thrombin inhibitors from herbs used in traditional Chinese medicine (TCM). A pharmacophore model and molecular docking were utilized to virtually screen a library of chemicals contained in compositions of traditional Chinese herbs, and these analyses were followed by in vitro bioassay validation and binding studies. Berberine (BBR) was first confirmed as a thrombin inhibitor using an enzymatic assay. The BBR IC50 value for thrombin inhibition was 2.92 μM. Direct binding studies using surface plasmon resonance demonstrated that BBR directly interacted with thrombin with a KD value of 16.39 μM. Competitive binding assay indicated that BBR could bind to the same argartroban/thrombin interaction site. A platelet aggregation assay demonstrated that BBR had the ability to inhibit thrombin-induced platelet aggregation in washed platelets samples. This study proved that BBR is a direct thrombin inhibitor that has activity in inhibiting thrombin-induced platelet aggregation. BBR may be a potential candidate for the development of safe and effective thrombin-inhibiting drugs.

important role in disease treatment [11][12][13][14] . In this study, we describe a combination of in silico and in vitro experiments that identified a small-molecule direct thrombin inhibitor from TCM. A library of 23,033 natural compounds were screened in silico through pharmacophore modelling and molecular docking. The top 23 hits were evaluated for thrombin inhibition with an enzymatic assay, and berberine (BBR) showed direct thrombin inhibitory activity. Additionally, a surface plasmon resonance (SPR)-based binding study and molecular docking were carried out to characterize the interaction between BBR and thrombin. A thrombin-induced platelet aggregation assay was conducted to evaluate the in vitro bioactivity of BBR. The strategy used in this work provided an effective and feasible approach for identifying direct thrombin inhibitors from natural products and could promote the development of safe and effective thrombin-inhibiting drugs.

Results
In silico screening for potential thrombin inhibitors. Ten pharmacophore models (Table S1) were generated based on the common features of six known direct thrombin inhibitors. Model assessment studies (Table 1) indicated that Model_10 (Fig. 1A) had the highest comprehensive appraisal index (CAI) and identified effective index (N), indicating that this model had the best ability to identify active compounds and exclude inactive compounds comprehensively 15 . Model_10 contained one H-bond acceptor (HBA, marked with green), one aromatic ring (AR, marked with yellow) and one hydrophobic group (HY, marked with cyan). The best active compound (CHEMBL377303) could map all features of Model_10 with a fit value of 3.00 (Fig. 1B). Model_10 was used to screen traditional Chinese medicine database 2009 (TCMD2009, Chinese Academy of Sciences), resulting in a hit list of 93 compounds (Table S2).
Molecular docking was performed on the basis of the crystal structure of thrombin (PDB ID 4UFD). The known direct thrombin inhibitor ligand S49 (molecular formula: C 29 H 33 N 5 O 4 S), the structure of which was solved as a co-crystallized thrombin complex, was re-docked into the active site of thrombin to validate the reliability of the docking protocol. The result showed that S49 could bind to thrombin via H-bond interaction with Asp189, Ser214, Gly216 and Gly219 (Fig. 2), which was consistent with previous studies 16 . The low root mean-square deviation (RMSD) of the re-docked and co-crystallized conformation of S49 was calculated to be 1.48 Å (Fig. 2 Table 1. Assessment results for each pharmacophore model. a Ht is the number of hits. b Ha is the number of active hits. c A% represents the ability to identify active compounds from the external database (A% = Ha/A, while A is the number of active compounds in the external database). d Y% represents the proportion of active hits in total hits (Y% = Ha/Ht). e N represents the ability to distinguish active compounds from non-active compounds (N = (Ha/Ht)/(A/D), while D is the number of compounds in the external database). f CAI was utilized to evaluate the models comprehensively (CAI = N × A%).

Figure 1.
The pharmacophore model_10 of thrombin inhibitors (A) and the matching pattern between pharmacophore model_10 and CHEMBL377303 (B). The numbers in (A) represent the distance between the two pharmacophore features. In (A) and (B), the arrows represent the direction of the hydrogen bond groups. Grey, red, blue and yellow atoms represent carbon, oxygen, nitrogen and sulfur atoms, respectively.
which indicated that the docking protocol established in this study could reasonably predict the binding mode of a known thrombin inhibitor. Ninety-three compounds identified by the pharmacophore model were docked into the active site of thrombin, resulting in a hit list of 39 compounds with docking scores above 5.0 (Table S3). Compared with pharmacophore-based virtual screening, 23 compounds (Fig. 3) both fit the pharmacophore model well and provided reasonable docking results.
In vitro screening for direct thrombin inhibitors. The

SPR-based binding studies. For a deeper investigation on the interaction between BBR and thrombin,
SPR-based binding analysis was performed. As shown in Fig. 6 and Table 2, BBR bound to thrombin (K D = 16.39 μ M) much tighter than argatroban (K D = 53.78 μ M) and hydroxyl alizarin (K D = 1,117 μ M), despite of their similar binding profile. We next tested whether BBR occupies the same binding site in thrombin as argartroban. A competitive binding assay was carried out and sensorgrams from binding of mixture of argartroban and BBR were compared with those from argartroban alone. The comparison showed that the binding of argartroban to thrombin was decreased in the presence of 300 μ M BBR. The inhibitory effect could be observed even with up to 120 μ M argartroban (Fig. 7). Altogether, the data indicated that BBR could bind to the same argartroban/thrombin interaction site with enhanced affinity.
Binding site of BBR in thrombin model. A ligand-thrombin docking model was used to reveal the mechanism underlying the molecular recognition between BBR and thrombin. The results showed that the C10 methoxy group of thrombin inserted into the catalytic centre of thrombin to form two hydrogen bond interactions with the side chain containing Phe227 and Trp215. The aromatic ring A of BBR formed π -π interactions with the side chain of Trp-60D (Fig. 8). The functional groups of BBR predicted via molecular docking were consistent with the pharmacophoric features generated with the pharmacophore model ( Fig. 9).
Antiplatelet activity of BBR. To evaluate the inhibitory effects of BBR on thrombin-induced platelet aggregation in vitro, washed platelets (WP) samples were pretreated with different concentrations of BBR. Then, 0.5 U/ml thrombin-induced platelet aggregation was monitored after the addition of test compounds for 5 minutes. As shown in Fig. 10, there were significant differences between the BBR groups and the negative control. BBR can significantly prevented platelet aggregation in WP samples.

Cytotoxicity evaluation.
To assess the safety of BBR for its biological applicability and therapeutic uses, the cytotoxicity of BBR was evaluated with a luciferase-coupled adenosine triphosphate (ATP) quantitative assay. Staurosporine, a prototypical ATP-competitive kinase inhibitor that can induce cell apoptosis 17 , was used as a positive control in this assay. Compounds were incubated with HEK293 cells for 24, 48 or 72 hours before luminescence signal measurement. Compared to the control group, BBR showed no significant cytotoxicity on HEK293 cells after incubation for 72 hours (Fig. 11). The values represent the means and standard error of replicates from three independent experiments.

Discussion
The combination of virtual screening, enzymatic bioassay and binding studies performed in this study showed remarkable advantages in identifying direct thrombin inhibitors from traditional Chinese herbs. Virtual screening focused on quick identification of potential active compounds from a large molecular dataset. The enzymatic bioassay focused on functional verification of direct thrombin inhibitors from the perspective of molecular biology. SPR-based binding studies focused on the physical binding kinetics of ligands and target. The platelet aggregation assay focused on in vitro biological efficacy. The potency of BBR was highly consistent with its ability to inhibit the thrombin target. Finally, BBR was confirmed as a direct thrombin inhibitor without cytotoxicity in HEK293 cells. Several methods were used to confirm that BBR is a direct thrombin inhibitor from multiple perspectives, providing a reference to reveal the structural and activity characteristics of a direct thrombin inhibitor.
BBR, an isoquinoline alkaloid, is a natural compound that exists in many natural herbal products, such as Coptis chinensis, Phellodendron amurense, Chelidonium majus, Stephania cepharantha, Hydrastis canadensis, etc. BBR is popular in TCM for its beneficial effects in human type 2 diabetes 18 , stroke 19 , hypertension 20 , inflammation 21 , microbial infection and cancer 22 , although its mechanism of action has not been studied clearly. This was the first work to identify BBR as a thrombin inhibitor, which was helpful in elucidating the mechanism underlying the anti-platelet aggregation pharmacological effect of BBR. This could be considered as a starting point for optimization to generate a dependable and efficient direct thrombin inhibitor in a follow-up study.
The essential characteristics of BBR were revealed by investigating the overlap between the pharmacophore model and docking studies. The C10 methoxy group of BBR acted as a key hydrogen bond acceptor by interacting with Phe227 and Trp215 of thrombin, which is consistent with the key amino acid residues in the active-site region of thrombin 23 . The aromatic ring A of BBR could interact with Trp-60D of thrombin through π -π interactions, suggesting that Trp-60D might be a critical amino acid residue that binds to BBR in the interaction process.
It took less than one month to carry out all the computational studies on a single computer. Together with the enzymatic bioassay, binding studies and platelet aggregation assay, the entire process took no more than two months and required significantly fewer resources than similar lead compound discovery steps in drug design programs. Therefore, this is a time-and resource-saving method for discovering direct thrombin inhibitors from natural products.

Methods
Pharmacophore-based virtual screening. A diverse dataset of 116 experimentally known thrombin inhibitors with IC 50 values below 1 μ M were obtained from published literature [24][25][26][27][28][29][30] . All the compounds were sketched and converted to 3D structures with all proton and MMFF94 charges added using SYBYL-X 1.2. To generate the pharmacophore model, six compounds that met the following criteria were selected as a training set  The Common Feature Pharmacophore Generation protocol in Discovery Studio v3.5 (Accelrys, San Diego, CA, USA) was used to generate the ligand-based pharmacophore models. A principal value of 2 and a maximum omit feature value of 0 was assigned to the six compounds in the training set. Energy minimization was performed using the CHARMM force field for all the compounds. Poling algorithm was used to generate a maximum of 255 diverse conformations with the threshold of 20 kcal·mol −1 above the calculated lowest energy for each compound in the training set. Conformers were generated using the diverse conformer generation protocol running with the best conformer generation option as available in DS. All the training set compounds associated with their conformations were used to generate the pharmacophore using the common feature pharmacophore generation module in DS. The feature mapping protocol was used to identify common features shared by the training set. As predictors, hydrogen-bond acceptor (HBA), hydrogen-bond donor (HBD), hydrophobic (HY) and aromatic ring

Table 2. The binding parameters of the interaction between immobilized thrombin and test compounds.
(AR) features were selected during pharmacophore generation. Ten possible pharmacophore models having one more arrangement of constituent features were generated in each pharmacophore running. The pharmacophore models were classified according to the ranking scores. Redundant hypotheses having the exact same chemical characteristics and nearly equal distances between these functions were deleted 31 .    D is the number of compounds in the external database. A is the number of active compounds. Ht is the number of hits and Ha is the number of active hits. A% represents the ability to identify active compounds from the external database. Y% represents the proportion of active hits in total hits. N represents the ability to distinguish active compounds from non-active compounds. CAI was utilized to evaluate the models comprehensively. The model with the highest CAI was used to screen TCMD2009.
A total of 23,033 compounds from TCMD 2009 were extracted and converted into 3D conformers using the CONCORD module in SYBYL-X 1.2 software (Tripos, Inc., St. Louis, MO, USA). The structures of all 23,033 compounds were energy minimized using the Tripos force field (Powell method and 0.05 kcal/mol Å energy gradient convergence criteria) and assigned an electrostatic charge with the Gasteiger-Hückel method. All the energy-minimized structures were stored as a 3D database in TCMD2009.
Virtual screening was performed using the search 3D database module in DS, with the minimum interference distance set to 1 Å and the search method set to best. All other protocol parameters were kept as the default settings. Fit value was calculated to indicate the matching degree of each ligand on the pharmacophoric features. A higher fit value suggests a better alignment between ligand's conformer and pharmacophore model.

Molecular docking-based virtual screening.
Surflex-Dock, a well-recognized method in the field of molecular docking 33,34 , was used to perform virtual screening and calculate the ligand-receptor interaction. The X-ray crystal structure of thrombin (PDB code: 4UFD. Resolution: 1.43 Å) resolved recently was designated to be the docking template 35 . Crystallographic water molecules in the structure were deleted and hydrogen atoms were added. Using the standard parameters employed in SYBYL-X 1.2 software package, the structure was energy minimized using AMBER7 F99 force field.
To check the accuracy of the docking program, ligand S49 was flexibly re-docked into the active site of thrombin. The RMSD of the re-docked and co-crystallized conformation of ligand S49 was calculated according to the following formula: where d is the distance between N pairs of equivalent atoms excluding hydrogens. A lower RMSD indicates a higher overlap between the re-docked and co-crystallized conformation of ligand S49. The active site of thrombin was defined as the protomol generated using the steric hydrophobic group (CH 4 ), the hydrogen bond group (C = O), and the hydrogen acceptor (N-H) within 0.5 Å of the ligand S49 binding site. After validation, the compounds identified by pharmacophore model were docked into the active site of thrombin using the default parameters in SYBYL-X 1.2 software. After each docking run, the best ten docked conformers were sorted in a molecular spreadsheet. They represent binding affinities in − log 10(Kd) based on surflex-dock scoring function (crash score (also pKd units), polar score, D-score, PMF-score, G-score, ChemSco and CScore) 36 . The common hits through ligand-and structure-based virtual screening were further evaluated using thrombin inhibition assay.

Thrombin inhibition assay.
To evaluate the thrombin-inhibitory activity of the hits, a fluorescence assay utilizing a SensoLyte 520 Thrombin Activity Assay Kit (AnaSpec, Inc., San Jose, CA) was performed following the manufacturer's instructions. This kit contained a novel 5-FAM/QXL 520 (fluorophore/quencher pair) FRET modified substrate (approximately 6 amino acids) that can be cleaved by thrombin at the Arg-Gly (R-G) site. The FRET substrate can be cleaved into two separate fragments, resulting in the release of 5-FAM fluorescence, which can be monitored at an excitation/emission = 490 nm/520 nm 37 . Argatroban (CAS No. 74863-84-6), a known direct thrombin inhibitor 38 , was used as a reference positive control for thrombin activity assay. All the test compounds were acquired from the National Institutes for Food and Drug Control (China) for thrombin inhibition assay. The purity of all the compounds was over 98% on the basis of HPLC analysis. . Luminescence was determined with an Envision 2100 multilabel reader to detect viability following incubation with the CellTiter-Glo reagent for 15 min. Compared to the control group, BBR showed no significant cytotoxicity on HEK293 cells after incubation for 72 hours. All error bars indicate the SE of three replicates. ns, no significant difference between the control group and the BBR groups. ***P < 0.001. Ten microlitres of compound solution (300 μ M in assay buffer) and 40 μ L thrombin solution (0.33 μ g/ml in assay buffer) were sequentially added to each well of a black, flat-bottom, non-binding 96-well plate (Corning #3650, Corning, NY). The plate was incubated for 10 min at 37 °C. Fifty microlitres of thrombin-substrate solution (6.67 μ M in assay buffer) was then added to each well. The reaction mixture was incubated for 1 hour at 37 °C. After incubation, the fluorescence (excitation 490 nm, emission 520 nm) of each mixture was detected with a 2104 EnVision Multilabel Reader (PerkinElmer).
For enzymatic reaction, the thrombin solution treated with vehicle buffer alone was used as negative control. The thrombin solution incubated with argatroban (60 μ M in assay buffer) was used as positive control. Furthermore, 100 μ L hits solutions (300 μ M in assay buffer) were served as test compound control to avoid false results due to compound autofluorescence. All fluorescence readings are taken into account in relative fluorescence units (RFU). The experimental data were analyzed by GraphPad Prism software (version 5, Graphpad Software Incorporation, USA).
Binding affinity measurement using SPR. SPR binding studies were performed using a Biacore T200 biosensor system (GE Healthcare, Uppsala, Sweden). All the SPR-based materials were acquired from GE Healthcare. Thrombin from human plasma (T6884), dimethyl sulfoxide (DMSO, W387520) were purchased from Sigma-Aldrich Chemical Co. (St. Louis, MO, USA).
Thrombin was diluted in 10 mM sodium acetate buffer at pH 5.0 and immobilized on a CM5 chip using an amine coupling kit. The immobilization of thrombin on a CM5 chip was performed according to the BIA applications handbook 39 . A total of 8419 RU of immobilized proteins was obtained. Sensor preparation and interaction analyses were performed at 25 °C in a PBS-P [10 mM phosphate buffer, 137 mM NaCl, 2.7 mM KCl, 0.05% P20 (pH 7.4)] running buffer containing 5% DMSO.
Tested compounds were prepared in a two-fold dilution concentration series by running buffer containing 5% DMSO (0.35-11.25 μ M for berberin, 4.7-75 μ M for argatroban and 1.17-37.50 μ M for hydroxy alizarin). Argatroban and hydroxy alizarin were set as positive control and negative control, respectively. Solvent-correction procedures were included to compensate for any DMSO-related bulk refractive index variations. Reference flow cell without immobilized thrombin served as a non-specific binding control. Biacore traces were baseline subtracted and the signal was presented in sensorgrams and measured in RU. Empirically in the BIAcore technology, 1 ng of analyte bound at the surface gives a response of 1000RU 2,40 . Equilibrium constants (K D ) were calculated using the 'affinity' model in Biacore T200 evaluation software version 2.0.
Using the same thrombin-immobilized CM5 chip, SPR competitive binding assays were performed following "ABA" method in a Biacore S200 biosensor system (GE Healthcare). Argartroban (30,60 or 120 μ M) were injected with or without BBR (300 μ M) over the surface of the chip for 60 seconds at 30 μ L/min. For test of mixture binding, PBS-P containing BBR (300 μ M) was used as running buffer instead of PBS-P alone. The response units were measured before the end of injection.

Inhibition of thrombin-induced platelet aggregation assay.
To further evaluate the anti-thrombin activities of BBR, inhibition of thrombin-induced platelet aggregation was measured using the turbidimetric method 41 . Male rabbits (weight 2.88 kg) were anesthetized with 10% chloral hydrate and blood was collected from the aortic neck and anticoagulated with sodium citrate (109 mM; 9:1, v/v). PRP was obtained by centrifugation of collected blood at 190 g for 10 min. The platelets were isolated from PRP by centrifugation for 10 min at 650 g, resuspended and washed twice in Ca 2+ -free Tyrode's buffer (138 mM NaCl, 2.68 mM KCI, 11.9 mM NaHCO 3 , 0.36 mM NaH 2 PO 4 , 0.49 mM MgSO 4 , 5.5 mM glucose, pH7.5) 42 . Washed platelets were then resuspended in Tyrode's buffer including 1.8 mM CaCl 2 at a density of 2 × 10 7 platelets/ml. The trial protocol was approved by the Ethical Committee of Capital Medical University, China.
In brief, 300 μ L WP was preincubated with BBR (final concentrations of 0.83, 2.50 and 7.50 μ M), argatroban (final concentrations of 0.62 μ M) or normal saline (control group) for 5 min at 37 °C. Aggregation was initiated by the addition of 10 μ L of thrombin (final concentrations of 0.5 U/ml). Aggregation was recorded at 37 °C for 5 min using a four-channel platelet aggregation analyser (LBY-NJ4, Beijing Precil Instrument Co., Ltd., China). Platelet aggregation inhibitory percentage was calculated using the following equation: Inhibitory percentage (%) = (1 − A t /A 0 ) × 100%. A t is the platelet aggregation percentage of the test sample, and A 0 is the platelet aggregation percentage of the control sample.
Luciferase-coupled ATP quantitative assay. HEK293 cells were seeded at 5.0 × 10 3 per well into 96-well clear-bottom black plates and incubated in 5% CO 2 at 37 °C overnight. Different concentrations of BBR or 1.0 μ M STSP were added to the 96-well plates and incubated with the cells in 5% CO 2 at 37 °C for 24 hours, 48 hours or 72 hours. Luminescence was determined with an Envision 2100 multilabel reader to detect viability following incubation with CellTiter-Glo reagent for 15 min.

Statistics.
For platelet aggregation and cytotoxicity evaluation experiment, statistical analyses were performed by one-way analysis of variance (ANOVA) followed by Dunnett's multiple comparison test. A value of p < 0.05 was considered as statistically significant.