Virtual Screening of Chemical Compounds for Discovery of Complement C3 Ligands.

The complement system is our first line of defense against foreign pathogens, but when it is not properly regulated, complement is implicated in the pathology of several autoimmune and inflammatory disorders. Compstatin is a peptidic complement inhibitor that acts by blocking the cleavage of complement protein C3 to the proinflammatory fragment C3a and opsonin fragment C3b. In this study, we aim to identify druglike small-molecule complement inhibitors with physicochemical, geometric, and binding properties similar to those of compstatin. We employed two approaches using various high-throughput virtual screening methods, which incorporate molecular dynamics (MD) simulations, pharmacophore model design, energy calculations, and molecular docking and scoring. We have generated a library of 274 chemical compounds with computationally predicted binding affinities for the compstatin binding site of C3. We have tested subsets of these chemical compounds experimentally for complement inhibitory activity, using hemolytic assays, and for binding affinity, using microscale thermophoresis. As a result, although none of the compounds showed inhibitory activity, compound 29 was identified to exhibit weak competitive binding against a potent compstatin analogue, therefore validating our computational approaches. Additional docking and MD simulation studies suggest that compound 29 interacts with C3 residues, which have been shown to be important in binding of compstatin to the C3c fragment of C3. Compound 29 is amenable to physicochemical optimization to acquire inhibitory properties. Additionally, it is possible that some of the untested compounds will demonstrate binding and inhibition in future experimental studies.

Introduction

The complement system, consisting of over 40 soluble and cell-bound proteins, is integral to innate immunity.[1−5] In the event of pathogen exposure or injury, cascading complement response occurs resulting in opsonization, chemotaxis, phagocytosis, and lysis.[6,7] Complement activity is also double-edged as a lack of regulation or balance in its response can be observed in numerous autoimmune and inflammatory diseases, including age-related macular degeneration, lupus, rheumatoid arthritis, multiple sclerosis, Sjögren syndrome, scleroderma, chronic obstructive pulmonary disease, ischemia reperfusion injuries, and rare diseases such as paroxysmal nocturnal hemoglobinuria, atypical hemolytic uremic syndrome, and C3 glomerulopathy, among others.[8−10] Currently, there are clinically available drugs for only two targets within the complement cascade, variations of the natural protein inhibitor C1-INH and a C5 inhibiting monoclonal antibody eculizumab, both of them being protein-based therapeutics.[9,10]

Compstatin[11] is a cyclic peptide capable of inhibiting complement response through C3,[12] originally discovered using a phage-displayed peptide library screening,[13] subsequently reaching clinical trials for age-related macular degeneration and other complement-related diseases.[14,15] We have been involved in structure- and computation-based optimization of compstatin family peptides, originally using a major structural conformer of free compstatin from solution NMR studies (reviewed in refs (16−21)) and subsequently using bound structures from computational de novo design studies and molecular dynamics (MD) simulations, on the basis of the crystal structure of a compstatin analogue bound to C3c[22] (e.g., see refs (23−28)). Although our most recent design has led to overcoming solubility/aggregation issues of the previously most potent compstatin analogues,[26,28] peptides in general suffer from low stability and bioavailability in vivo and often require intravenous administration. Chemical compounds are typically orally administered and are more cost effective for scaled-up industrial production. Therefore, there is a need for the development of nonpeptidic low-molecular-mass inhibitors. Toward this goal, we launched a pharmacophore-based virtual screening study, described here, to identify druglike chemical compounds with the geometric and physicochemical characteristics of compstatin, which are capable of binding to the compstatin binding site of target protein C3.

Virtual screening has proven to be a valuable methodology for identifying potential therapeutic candidates[29] and an alternative to fragment-based chemical compound design[30] or rational peptide design.[31−33] Virtual screening has the benefit of being a high-throughput method and can be used to screen millions of chemical compounds, while being more time and resource efficient than experimental high-throughput screening methods, thanks to the advances in computer hardware architecture and drug design-related algorithms.[34−36] Small druglike compounds are desirable because they typically exhibit better pharmacological properties, but at the expense of lower specificity, compared with peptide- or protein-based therapeutics. In this study, we utilize virtual screening with the objective to identify novel druglike compounds capable of binding in the compstatin binding site of C3 and potentially inhibiting complement response. We use a molecular dynamics structure of a potent compstatin analogue bound to C3c as the basis for the development of pharmacophore models and docking of molecules. Our virtual screening framework is similar to that used in a recent identification of 11 druglike ligands of complement fragment C3d, 10 of which are fluorescent markers of complement activation.[37]

Methods

Primary Approach

Pharmacophore Models

A pharmacophore model is represented by a framework of features corresponding to the spatial distribution of physicochemical properties (aromaticity, hydrophobicity, hydrogen bond donor/acceptor capability, and positive/negative charge) of an active ligand. During screening, a database of molecules is compared against the pharmacophore model and if the molecule contains matching pharmacophore features, then it is considered a positive hit.

The workflow of the primary approach procedure is shown in Figure . Pharmacophore models were developed using a molecular dynamics (MD) trajectory of the complex between C3c and the RSI-compstatin analogue with sequence Ac-RSI[CVWQDWGAHRC]T-NH2 (brackets denote cyclization through a disulfide bridge).[26] Pharmacophore features were primarily selected from earlier molecular dynamics data of C3c-bound compstatin analogues,[27,38] with the aid of prior knowledge from MD studies[24,39] and optimization studies[23−25,27] of C3c-bound compstatin analogues and optimization studies of free compstatin analogues.[16] All MD simulations were based on the crystal structure of C3c with the W4A9 analogue of compstatin.[22]

Figure 1

Flowchart of the primary virtual screening approach.

Mechanistic binding analysis of compstatin structure throughout the MD trajectory was performed using the R package Bio3D, the Python library MDTraj, and Chimera[40−43] to identify significant physicochemical properties and nonpolar contacts at the binding site. In addition, free-energy contributions of individual amino acids and hydrogen bond occupancies from previous studies[27,38] informed the selection of pharmacophore features. Mean positions of centers of mass were calculated for the atoms identified in each selected feature. The tolerance radii of each pharmacophore feature were defined by calculating the conformational flexibility of the feature from the MD trajectory. In the first round of screening, 473 pharmacophore models were developed using subsets of 3–5 features identified to be of interest. The purchasable subset of the ZINC 12 database,[44] consisting at the time of the study of ∼19 million molecules (190 million conformers) in stock or to be made-to-order, was screened using ZINCPharmer[45] with each of the 473 pharmacophore models. A molecule was a positive hit during screening if its chemical moieties were spatially distributed such that there was overlap with the tolerance radii of the specific features of the pharmacophore model.

In the second round of screening, 40 new pharmacophore models consisting of subsets of 3–6 features were developed with iterative improvements over the initial 473 models. A set of ∼7 million molecules with ∼1.1 billion conformers that was used in one of our previous virtual screening studies[37] was used for screening of these 40 pharmacophore models using Phase.[46,47] This set of molecules was from the Drugs Now subset of the ZINC 12 database, consisting of molecules that were in stock at the time of the study, and the conformers were generated using Phase.

Docking

Molecules identified to fulfill the selection criteria of the pharmacophore models as positive hits during the pharmacophore screening were docked to C3c in the binding region of RSI-compstatin. In the first round of pharmacophore screening and docking, the molecules were docked to a single conformation of C3c (acquired from the final frame of the MD trajectory of the C3c/RSI-compstatin complex). In the second round of screening, representative conformational states of the binding site of C3c were extracted from the MD trajectory of the C3c/RSI-compstatin complex to accurately capture the conformational variations of the binding site. The conformational states of C3c observed in each frame of the MD trajectory were superimposed on the basis of Cα atoms and hierarchical clustering was performed on the basis of the root-mean-square deviation (RMSD) of C3c amino acids identified to be involved in the interaction with RSI-compstatin.[27] Five clusters were calculated, corresponding to five representative structures of C3c. Molecular docking was performed using AutoDock Vina[48] with preprocessing of structures using AutoDock Tools. Each molecule was docked to the representative structures of C3c within a box (with dimensions of 28 Å × 28 Å × 28 Å), encompassing the entire RSI-compstatin binding site of C3c. The exhaustiveness parameter of AutoDock Vina was set to 20 to improve docking accuracy, and the top 20 docked poses of each molecule, based on predicted binding energies, were returned.

Scoring

Each of the docked poses were scored using the Vina scoring function, and the predicted binding energies were reported. Mean predicted binding energies were calculated for docked poses to each of the five representative structures of C3c. The predicted solubility of the molecules was calculated using the partition coefficient (log P) using the ChemmineR package in R.[49] Molecules were also evaluated using ChemmineR, visually inspected with Chimera, and verified through the ZINC 12 database for adherence to Lipinski’s Rule of Five.[50] A combination of the above scoring methods, in addition to visual inspection of the molecules for geometric properties and occurrence of specific pharmacophore features, were utilized to select 58 compounds for ordering and experimental testing.

Alternative Approach

In an alternative approach, two stages of additional searches were performed. In stage 1, two distinct pharmacophore searches (search A and search B) were performed. In stage 2, an additional search was performed on the basis of the structures of the most promising molecules identified in stage 1. The pharmacophore features and the targeted C3c residues used in this approach are shown in Appendix S8, Table S2. All searches were performed using ZINCPharmer.[45] The workflow of the alternative approach procedure is shown in Figure .

Figure 2

Flowchart of the alternate virtual screening approach.

Stage 1

Two separate searches were performed in the first stage. The two separate searches differ in constraints to the allowable molecular weight of the selected molecules. In search A, the results were limited to molecules with a molecular weight less than 500. In search B, the results were limited to molecules with a molecular weight greater than 500. It has been shown that molecules with a molecular weight less than 500 have a greater bioavailability,[50] but due to the size and complexity of the compstatin binding site, molecules with a molecular weight greater than 500 were also considered in this study.

For both search A and search B, pharmacophore models were generated on the basis of the resolved crystal structure of compstatin bound to human C3c (PDB code 2QKI(22)) as well as docked structures we produced from C3c in complex to cmp-5 (ZINC61197239), a compound recently shown to significantly inhibit the cleavage of C5,[51] a downstream process of the complement system. The pharmacophore models based on the binding of compstatin to C3c were developed using subsets of four features, each of which target an amino acid within one of four amino acid sectors of human C3c (sector I: 344–349, sector II: 388–393, sector III: 454–462, and sector IV: 488–492). In this way, molecules encompassing the pharmacophore features were expected to bind to each of the four aforementioned C3c sectors. Amino acids within these sectors have been shown to be in direct contact with compstatin,[1−3,22] and previous studies identified individual amino acids within these sectors to be key to compstatin binding.[24,38] To create the additional pharmacophore models, cmp-5 was first docked to the compstatin binding site (PDB: 2QKI(22)) using AutoDock Vina. The generated docked poses of cmp-5 binding to C3c were then clustered on the basis of RMSD. The docked poses acquiring the most favorable binding energies within separate clusters were selected to generate the pharmacophore models based on cmp-5 binding to C3c. Only molecules containing less than 10 rotatable bonds were considered for investigation. Molecules that fit the tolerance radii of each pharmacophore feature were analyzed further on the basis of charge complementarity, RMSD, and a visual inspection of each molecules’ structure-based interactions with C3c to determine the potential of the molecule as a possible C3 inhibitor. In the case that a pharmacophore model had over 500 matches, the first 500 molecules, ranked by the number of rotatable bonds from least to greatest, were considered. Molecules with a negative net charge were excluded from further investigation because the most potent compstatin analogues, including those used in this study, have positive or neutral net charge. In the case that a molecule had multiple conformations fitting a pharmacophore model, the conformation with the lowest RMSD to the pharmacophore model was selected. Molecules were visually inspected for structural interactions with the compstatin binding site. Only molecules within the purchasable subset of the ZINC 12 database[44] were considered for further analysis.

Single MD simulations were conducted for the molecules identified by each of the two searches and the computationally derived binding affinity of the molecules was assessed on the basis of average binding energies estimated by the Vina scoring function.[48] The top five computationally derived molecules with the lowest average binding energies were used to develop an additional set of pharmacophore models for the second stage of searches (described below).

Stage 2

In the second stage of the alternative approach, the top five molecules with the lowest average binding energies estimated by the Vina scoring function[48] out of all of the molecules from the two searches in stage 1 were selected to develop an additional set of pharmacophore models for the second stage of searches. From the simulations of the selected five molecules in complex with C3c, the lowest binding energy snapshot was extracted and used to generate additional pharmacophore models. The additional pharmacophore models consisted of at least three features and target at least three of the four compstatin binding sectors. Additional single molecular dynamics simulations were conducted for the molecules identified in this stage, and the computational binding affinity of the molecules was assessed on the basis of the average binding energies from the simulation snapshots (described below).

Molecular Dynamics Simulations

Molecular dynamics (MD) simulations were conducted for the selected molecules fulfilling the selection criteria described above in complex with C3c independently. The MD simulations were employed to refine receptor ligand conformations, assess the energetic favorability of the molecules binding to C3c, and determine structural stability of the molecules within the compstatin binding site of human C3c.

The initial coordinates for C3c were obtained from the PDB 2QKI.[22] To decrease computational demand, the protein was truncated for the simulations.[38] The simulated system included the entire MG4 and MG5 domains of the compstatin binding site (amino acids 329–534) and segment 607–620, which is proximal to MG4 and MG5. The initial docked position of the molecules in complex with C3c corresponded to the coordinates of the molecule oriented by ZINCPharmer[45] during the pharmacophore search. CHARMM version c39b2[52] was used to set up and perform all simulations.

The MD simulations were performed using the CHARMM36 topology and parameters.[53] Initial energy minimization steps were performed on each individual molecule/C3c complex by sequentially performing 200 steps of steepest decent minimization, 200 steps of adopted basis Newton–Raphson minimization, and 200 steps of steepest decent minimization. During the two energy minimization steps, the ligand, excluding hydrogens, was subjected to 0.5 kcal/(mol Å2) harmonic constraints and the protein, excluding hydrogens, was subjected to 1.5 kcal/(mol Å2) harmonic constraints.

For each of the explicit-solvent MD simulations, the molecule/C3c complex was solvated in a truncated octahedral water box with a length of 100 Å and a volume of approximately 769 800 Å3. For example, as a representative of all simulation systems, in the explicit-solvent MD simulations of compound 29, discussed later on, in complex with C3c, 22 895 water molecules were used to solvate the complex and the simulation system contained a total of 72 357 atoms. Because of the size of the simulation systems, the duration of all simulations was limited to 20 ns except for the best binding mode (see Refined Docking and MD Simulation Investigation of the Compound 29/C3c Complex Structure) of compound 29, which was simulated for 50 ns (see Results and Discussion). The potassium chloride concentration in the water box was set to 0.15 M, and additional potassium and chloride ions were added to neutralize the charge of the system. The ions were placed through 2000 steps of Monte Carlo simulations.[54,55] An energy minimization of 100 steps of steepest decent minimization, 100 steps of adopted basis Newton–Raphson minimization, and 100 steps of steepest decent minimization were sequentially performed on the solvent molecules.

The system was subsequently equilibrated for 1 ns with 1 kcal/(mol Å2) harmonic constraints on the protein backbone and ligand heavy atoms and 0.2 kcal/(mol Å2) harmonic constraints on the protein side-chain heavy atoms. After the equilibration stage, all constraints were released and amino acids greater than 20 Å away from any atom in compstatin according to the crystal structure[22] as well as truncated termini were subjected to 0.5 kcal/(mol Å2) of harmonic constraints for backbone atoms and 0.1 kcal/(mol Å2) for side-chain heavy atoms. At this stage, the system was simulated for 20 ns and frames were extracted every 20 ps. The first 10 ns of this stage were treated as further equilibration to allow the energy and ligand conformation to converge. The last 10 ns of the final stage were treated as the production stage. Upon completion of the simulations, the water molecules and ions were removed and the trajectory was analyzed as follows. Triplicate simulations were performed for the five molecules with the lowest average binding energy values across both stages of the alternative approach to ensure reproducibility. All simulations were performed using the Leapfrog Verlet integrator at isothermal and isobaric conditions. The temperature of the simulations was held at 300 K using the Hoover thermostat. Fast table lookup routines were applied to nonbonded interactions, and the SHAKE algorithm was implemented to constrain the bond lengths to hydrogen atoms.[56,57] Periodic boundary conditions were applied in each simulation.

Upon completion of the MD simulations, the binding of the candidate molecules from both stages of the alternative approach to human C3c was assessed using the average binding energy estimated by the Vina scoring function.[48] Each simulation snapshot extracted in increments of 20 ps from the last 10 ns of the MD simulation for each molecule/C3c complex was scored using the Vina scoring function,[48] and the binding energies were averaged over the total number of simulation snapshots analyzed. The average RMSD of the simulated molecules from both stages of the alternative approach in complex with C3c was also calculated to estimate their stability within the binding site. The RMSD calculations were performed using the heavy atoms of the molecules and with respect to the average structure of the last 10 ns as a basis through Wordom.[58] The RMSD of the molecules within their respective simulations was calculated for each simulation snapshot in increments of 20 ps for the final 10 ns of the MD simulation and averaged over the total number of snapshots analyzed. Molecules with an average RMSD value larger than 3.0 Å were discarded and omitted from further investigation.

The six compounds (see Results and Discussion) with the lowest average binding energies that exhibit stable binding throughout their simulations (having average RMSD values with respect to their average simulation structure of less than 3.0 Å) that were also readily available for purchase were evaluated using hemolytic assays and microscale thermophoresis (MST) described in the following sections.

Experimental Validation

Hemolytic Assays

Selected compounds were obtained from ChemBridge (compounds 1–55), Asinex (compounds 57 and 58), MolPort (compounds 56, A–C, E, and F), and Specs (compound D). Stock solutions were prepared by dissolving each compound in dimethyl sulfoxide (DMSO) to concentrations of 4 mM.

Rabbit erythrocytes (Complement Technology, Inc.) were washed in phosphate-buffered saline and resuspended in a veronal-buffered saline solution containing 5 mM MgCl2 and 10 mM ethylene glycol tetraacetic acid (EGTA) (VBS-MgEGTA). Each compound was diluted in VBS-MgEGTA to end up with a final concentration of 1% DMSO and was added to round-bottom 96-well plates. Normal human serum (NHS) diluted in VBS-MgEGTA was added to each well, and the plates were incubated at room temperature for 15 min. Thirty microliters of rabbit erythrocytes at a concentration of 1.25 × 108 cells/mL were then added to each well. Positive controls for lysis included rabbit erythrocytes in deionized water and rabbit erythrocytes in NHS diluted in VBS-MgEGTA, whereas negative controls included rabbit erythrocytes in VBS-MgEGTA and VBS-EDTA (20 mM EDTA). Following the addition of rabbit erythrocytes, the plates were incubated at 37 °C for 20 min and then quenched with ice-cold VBS with 50 mM EDTA. After centrifugation at 1000g for 5 min, the supernatant from the plates was diluted 1:1 with deionized water in flat-bottom 96-well plates and the lysis was quantified spectrophotometrically at 405 nm. The assays were performed in triplicate to ensure reproducibility.

Microscale Thermophoresis

The binding affinity for C3c of the top 10 compounds identified in the primary approach and the top 6 compounds identified in the secondary approach were evaluated in a competitive microscale thermophoresis (MST) assay, using a Monolith NT.115 instrument (NanoTemper Technologies GmbH, Munich, Germany). Competition was performed against a competition peptide, synthesized by ELIM Biopharm (Hayward, CA). The competition peptide was labeled with the cyanine fluorophore CY5, which was attached at the side chain of the preceding lysine, and had sequence Ac-I[CVWQDWGAHRC]TAGK-(CY5)-NH2. The peptide was cyclized by a disulfide bridge between the two cysteine amino acids and acetylated at the N-terminus and amidated at the C-terminus. A 1:1 serial dilution series of each of the selected compounds was performed in MST buffer (50 mM Tris–HCl, 150 mM NaCl, 10 mM MgCl2, and 0.05% Tween 20) with 5% DMSO. Each dilution series started with a final concentration of 500 μM and ended in a final concentration of 15.3 nM. To each dilution series of a compound, purified C3c (Complement Technology) and the competition peptide (labeled with Cy5) were dissolved to a final concentration of 20.4 nM (C3c) and 50 nM (competition peptide). The resulting mixture was incubated for 15 min in the dark at room temperature. Following incubation, the samples were loaded into standard capillary tubes and the thermophoretic response of the fluorescently labeled marker was measured. Each dilution series was performed in triplicate and estimation of the IC50 was performed through nonlinear regression.

The fragment C3c that was chosen for the MST assay binds compstatin but it does not contain the thioester domain. The choice of C3c for the MST assay is because the only cocrystal structure of a compstatin family peptide bound to C3/C3 fragment reported in the literature is with C3c,[22] and the sequence of the peptide is Ac-I[CVWQDWGAHRC]T-NH2. The bound compstatin analogue of the crystal structure is the parent peptide of the competition peptide.

We utilized MST to determine protein–ligand binding affinities.[59,60] MST measures binding that occurs in a bulk solution, avoiding artifacts encountered in methods similar to surface plasmon resonance where one binding partner must be immobilized onto a chip. Although other methods such as isothermal titration calorimetry (ITC) and fluorescent polarization (FP) can also be used to measure binding in solution, these methods are not as suitable for our application as MST. ITC is a low-throughput method that compensates for low sensitivity by increasing the amount of sample measured. This strategy is prohibitive in cost for evaluating multiple protein–ligand pairs. FP, on the other hand, can be performed in a high-throughput manner but would require fluorescent labeling of small molecules to detect changes in polarization of emitted light when binding occurs. As our molecules have small molecular weight, less than 500 Da, labeling is likely to perturb protein–ligand interactions. In contrast to ITC and FP, MST is more sensitive, requiring minimal sample volume (20 μL for a dilution series), and detects changes in surface area, hydration entropy, and net charge. Our MST instrument requires fluorescent labeling of one binding partner; labeling is less likely to perturb binding of a small molecule as the target protein with a large molecular weight can be labeled. Although not strictly a high-throughput device, our MST device is relatively fast and can determine a dissociation constant for a single protein–ligand interaction within the period of an hour.

Results and Discussion

Primary Approach

Virtual High-Throughput Screening

The objective for our study was to identify low-molecular mass molecules that are capable of binding complement protein C3 or its activation fragments, C3b/C3c, in the binding site of compstatin. Our study was based on previous knowledge of the structure of free and bound compstatin and many computational and experimental studies, which pointed to key structural and physicochemical features of compstatin that are important for binding to C3/C3b/C3c and for inhibiting the complement system. Our approach involved the development of pharmacophore models, pharmacophore-based virtual screening of conformationally flexible molecules, docking of pharmacophore-matched molecules to multiple conformations of the C3/C3b/C3c binding site, and scoring of docking poses using energetics, lipophilicity, and Lipinski’s rule of five criteria. Figure shows a schematic flowchart of our approach.

Figure 3

Diagram outlining an example of the pharmacophore model generation and docking output. (A) Molecular graphics showing the binding interaction between compstatin and C3c. (B) A pharmacophore model generated through selection of features identified from the C3c–compstatin interaction. (C) One of the molecules identified through the pharmacophore model superimposed on the structure of compstatin. (D) The same molecule as in (C) docked on C3c.

In the first round of virtual screening with the initial 473 pharmacophore models (Appendices S1 and S2), we were able to identify specific combinations of features that were likely to yield molecules of interest. The screening of the purchasable subset of the ZINC 12 database with the pharmacophore models suggested that features corresponding to R(−1), V3, W4, Q5, W7, A9, and H10 on RSI-compstatin were necessary as those features resulted in positive hits that had favorable predicted binding affinities as well as a proper docked fit when visually inspected. A positive charge pharmacophore feature was chosen at the position of R(−1). The R(−1)-S0 N-terminal modification was chosen to improve solubility and was found to retain inhibitory activity.[26,27] The addition of the R(−1) side chain introduced an ionic interaction with E372 of C3c (Figure A), which contributed in the binding affinity and stability of the C3/RSI-compstatin complex.[26,27] The hydrophobic character of V3 was included as a pharmacophore feature, as V3 inserts into a hydrophobic subcavity in C3c (Figure A). The amino acids W4 and W7 were observed in the MD trajectory to participate in highly conserved hydrogen bonds with C3c (Appendix S3) and as a result, hydrogen bond donor and acceptor features at their corresponding positions were included in pharmacophore models. Additionally, the aromatic and hydrophobic properties of W4 and W7 were utilized as pharmacophore features due to the pervasiveness of these features in druglike molecules and their role in favorable interactions with C3c (Figure A). Amino acids Q5, A9, and H10 are participating in hydrogen bonds (Appendix S3) and were included as pharmacophore features. The specific locations of A9 and H10 elongate the pharmacophore models, therefore increasing the screening diversity.

Figure 4

Molecular graphics of (A) compstatin bound to C3c and (B–K) the top 10 compounds (of the selected 58 from the primary approach) docked to the corresponding representative conformation of C3c. In (A), three main features of compstatin, V3, W4, and W7, important to the interaction with C3c are identified with dashed circles.

On the other hand, specific features or combinations of features were identified to be either too lenient as screening parameters or unlikely to exist in drug molecules. For example, pharmacophores that included more than one hydrophobic feature resulted in too many molecules matched. Features such as the hydrogen bond donor/acceptor capability of D6 in compstatin were found to result in positive hits (∼10 000) during the pharmacophore screen but did not result in viable predicted binding energies (>−6.5 kcal/mol) during docking.

After screening, docking and scoring, features identified to result in potentially viable molecules were used to iteratively improve upon the initial pharmacophore models, resulting in 40 new pharmacophore models. Upon screening the conformers generated from the Drugs Now subset of ZINC and docking the resulted hits to C3c, we identified ∼81 000 docked conformer poses. We refined the list of molecules by applying a threshold value of <−7 kcal/mol in predicted binding energies and a log P threshold of 5. In a few cases, molecules exhibiting log P values >5 but significantly favorable predicted binding energies were considered as well. Further filtering was performed by evaluating the molecules using Lipinski’s rule of five and visual inspection for physicochemical and geometric properties. A final list of 167 molecules were identified (Appendix S4), out of which 58 (Appendices S5 and S6) were purchased and tested experimentally. Figure B–K shows the top 10 compounds out of the selected 58 and displays the variety in physicochemical properties and spatial placement when docked. The top 10 compounds are ZINC72382898 (compound 58), ZINC72382894 (compound 57), ZINC67742743 (compound 29), ZINC29862046 (compound 56), ZINC14995377 (compound 6), ZINC67881194 (compound 41), ZINC12000754 (compound 1), ZINC67605047 (compound 14), ZINC12079160 (compound 2), and ZINC67974289 (compound 51). Compstatin has three main features of importance to the interaction with C3c: the V3 interaction with a smaller hydrophobic subcavity, the W4 interaction with a steric wall-like structural feature, and the W7 interaction with a larger hydrophobic subcavity (Figure A). A few of the top 10 compounds (compounds 56, 6, and 51) match all three of these features, whereas the other compounds match at least one of these features. The other compounds in the selected 58 (Appendix S7) also match at least one of the main features of compstatin, except for two compounds (compounds 28 and 16) that were docked to another region.

Initial Screening Based on Pharmacophore Models

A total of 107 molecules were selected from both stages of the alternative approach and further investigated using MD simulations and energy calculations (Appendix S4). From the first stage of the alternative approach, a total of 84 molecules were selected for MD simulations and binding energy calculations (Appendix S4, stage 1), of which 50 originated from search A and 34 originated from search B. The targeted C3c amino acids and the pharmacophore features used to target them are listed in Appendix S8, Table S2. Of these 84 molecules, the 5 molecules acquiring the lowest average binding energy were selected to generate additional sets of pharmacophore models. The selected molecules are ZINC64623185, ZINC1060630, ZINC39961116, ZINC08437931, and ZINC12558945, of which the first three listed molecules originated from search A and the last two listed molecules originated from search B. From each of the MD simulations of these five molecules binding to C3c, the lowest binding energy snapshot was extracted and used to develop additional pharmacophore models used for the second stage of the alternative approach. From the second stage of the alternative approach, a total of 23 additional molecules were selected for MD simulations and energy calculations (Appendix S4, stage 2).

The six compounds that were readily available for purchase, with the lowest average binding energies across all other molecules from both stages in complex with C3c and average RMSD values less than 3.0 Å, were selected for experimental testing. The six selected compounds are ZINC64623185 (compound A), ZINC63743940 (compound B), ZINC08437931 (compound C), ZINC01060630 (compound D), ZINC13688614 (compound E), and ZINC13628667 (compound F). Compounds A–C were selected from the first stage of the alternative approach, and compounds D–F were selected from the second stage of the alternative approach. Molecular graphics images of the six purchased compounds, in complex with C3c, and their simulation-based conformations within the compstatin binding pocket of C3 are shown in Figure .

Figure 5

Molecular graphics images from the MD simulations of the six compounds (compounds A–F, shown respectively in panels A–F) discovered in the alternative approach, which were selected for experimental testing. The images show representative snapshots of the compounds’ simulation-based conformations during the final 10 ns of the triplicate MD simulations. The compounds are shown in licorice representation. C3c is shown in tube representation, with key interacting amino acids shown in thin licorice representation. The colored amino acids represent the four compstatin binding sectors. From left to right the sectors are: VI (gray), III (yellow), I (blue), and II (orange). The dotted black lines represent hydrogen bonds present throughout the MD simulations.

Compound A was selected using pharmacophore features targeting C3c amino acids R459, H392, M346, and D491. These pharmacophore features used for the selection and initial orientation of compound A in the compstatin binding site were a hydrophobic feature targeting M346, an aromatic feature targeting H392, a negative charge targeting D491, and a hydrogen bond donor targeting R459. During all three MD simulations, interactions of compound A to M346, R459, and D491 are maintained, whereas the interaction to H392 is not. During the simulations, the aromatic ring of compound A, originally targeting H392, moves to form a cation−π interaction with R456; thus, compound A cannot maintain interactions to any amino acid of sector I (Figure A).

Compound B was selected using pharmacophore features targeting C3c amino acids R456, N390, M346, and D491. These pharmacophore features used for the selection and initial orientation of compound B in the compstatin binding site were a hydrogen bond acceptor targeting R456, a hydrogen bond acceptor targeting N390, a hydrophobic feature targeting M346, and a hydrogen bond donor targeting D491. Within all three MD simulations, interactions between compound B and R456, N390, and M346 are not maintained while a hydrogen bond to D491 is formed and maintained. The aromatic ring of compound B originally positioned to target N390 shifts to form cation−π interactions with both R459 and R456 as the simulations progress; thus, interactions to M346 and N390, as well as to any amino acid of sectors I or II are lost (Figure B).

Compound C was selected using pharmacophore features targeting C3c amino acids R459, N390, H392, M346, and D491. The pharmacophore features used for the selection and initial orientation of compound D in the compstatin binding site were a hydrophobic feature targeting M346, a hydrogen bond donor targeting N390, an aromatic feature targeting H392, a hydrogen bond acceptor targeting D491, and an aromatic interaction targeting R459. During all three MD simulations, all interactions were maintained except for the interaction targeting D491 (Figure C).

Compound D was selected using pharmacophore features targeting C3c amino acids M346, R456, and N390. These pharmacophore features used for the selection and initial orientation of compound C in the compstatin binding site were a hydrophobic feature targeting M346, a hydrogen bond acceptor targeting R456, and a hydrogen bond acceptor targeting N390. During the MD simulation, the compound maintains a hydrogen bond to N390, a hydrogen bond to R456, and hydrophobic interactions targeting M346. All interactions from the initial docking were conserved; however, the compound does not strongly interact with any amino acid of sector VI (Figure D).

Compound E was selected using pharmacophore features targeting C3c amino acids R456, N390, M346, and D491. The pharmacophore features used for the selection and initial orientation of compound E in the compstatin binding site were a hydrophobic feature targeting M346, a hydrogen bond donor targeting N390, a hydrogen bond donor targeting D491, and a hydrogen bond acceptor targeting R456. During all three MD simulations, interactions of compound E to M346, R456, and N390 are maintained, whereas the interaction to D491 is not. In addition, the compound forms a weak cation−π interaction to R459 (Figure E).

Compound F was selected using pharmacophore features targeting C3c amino acids R459, N390, M346, and L492. The pharmacophore features used for the selection and initial orientation of compound F in the compstatin binding site were a hydrophobic feature targeting M346, an aromatic feature targeting N390, a hydrophobic feature targeting L492, and a hydrogen bond acceptor targeting R459. During all three MD simulations, interactions from compound F to M346 and R459 are maintained, whereas the interaction to L492 and N390 are not. During the simulations, the aromatic ring of compound F, originally targeting L492, moves to form a cation−π interaction with R456 (Figure F).

Experimental Validation and Analysis

Hemolytic Assay

A standard rabbit erythrocyte assay was used to initiate complement activation and to evaluate possible inhibitory effects of the virtual screening chemical compounds. In this assay, complement, as part of normal human serum, is activated by the foreign rabbit erythrocyte cells, resulting in cell lysis by the membrane attack complex. Addition of chemical compounds with potential inhibitory activities would result in decreased lysis, compared with experiments without inhibitors, as shown before in the case of compstatin analogues.[26−28]

Despite demonstrating good predicted binding properties, none of the 58 selected compounds from the primary approach and the 6 selected compounds from the alternative approach exhibited inhibitory activity in the hemolytic assays (Figure ). Compstatin analogues are fairly long, cyclic, and bulky peptides (13–15 amino acids of 1500–1900 Da molecular weight), so it may be unlikely for a single druglike compound (<500 Da molecular weight) to retain the needed intermolecular contacts necessary for inhibition. It is possible, however, that the identified compounds are capable of binding to C3c, without demonstrating inhibition of hemolytic activity. In such a case, possible combination of compounds using chemical synthesis methods could result in a larger molecule with potential binding affinity for C3c and complement inhibitory activity. Also, compound 29, and perhaps additional experimentally untested compounds from the library, may be used as a scaffold to design new and larger molecules with favorable binding and inhibitory properties.

Figure 6

Evaluation of hemolytic assay results. The results of the hemolytic assays are represented as a comparison of means using one-way analysis of variance between the controls and the drug compounds. The bottom row of the figure shows the comparison between the two controls for lysis (averaged from both types of positive controls for lysis) and no lysis (see Methods), demonstrating a significant difference in means. The other rows show that the difference between lysis and the inhibitory effects of the compounds is negligible. The means are for data acquired through triplicate hemolytic assays.

Microscale Thermophoresis Binding Assay

We measured direct binding of the top 10 compounds identified in the primary approach and the top 6 compounds identified in the secondary approach 2 to C3c using the microscale thermophoresis assay, as described in Methods. Of the selected compounds, only compound 29 demonstrated binding to C3c. Figure shows competitive binding of compound 29 to the C3c–competition peptide complex, where the competition peptide was labeled with the fluorophore Cy5 for detection of the thermophoresis signal. This competitive binding experiment shows that compound 29 binds to C3c, albeit with a low affinity. An accurate measurement KD was not possible because the dose–response binding curve requires additional data points to reach a plateau at higher concentrations, and this was limited by the stock concentration of the compound.

Figure 7

Concentration-dependent binding curve of compound 29 to C3c in competition with the competition peptide. Thermophoretic data is plotted as mean ± standard error of the mean (as error bars) from three replicate experiments, together with the fitted binding curve in red and the 95% confidence interval of the fitted binding curve represented as blue dots.

Refined Docking and MD Simulation Investigation of the Compound 29/C3c Complex Structure

As compound 29 demonstrated binding to C3c, a refined docking study was performed, in conjunction with MD simulations, to elucidate the structural features of the complex between compound 29 and C3c. Compound 29 was redocked to C3c using AutoDock Vina with a larger search space and higher exhaustiveness compared with our original docking studies, to decrease the probability of not finding the docked pose at the global minimum according to the Vina scoring function.[48] The generated docked poses were clustered, and from the clustering analysis, two clusters of poses were identified, of which one was oriented in a fashion similar to that shown in Figure D and the other was flipped along the short axis of compound 29. From each of the two clusters, the pose with the most favorable binding energy was extracted and selected for further investigation. An in-house docking protocol[61,62] was used to produce additional poses with variable orientations, using the two poses extracted from AutoDock Vina (see above) as initial structures for further MD-based refined docking investigation. As a result, an additional 48 poses (24 for each of the two initial structures) were generated. Values of 1.0, 2.0, 3.0, and 4.0 Å roffset were used for docking simulations using both the harmonic spherical potential and the quartic spherical potential.

Explicit-solvent 10 ns MD simulations were performed on the two selected poses generated by AutoDock Vina as well as the 48 docked poses generated from the docking protocol.[61,62] The molecular mechanics-generalized Born surface area (MM-GBSA) association free energy[63] was calculated in accordance to the docking protocol,[61,62] and the simulation with the docked pose generated by AutoDock Vina originating from the cluster of poses flipped along the short axis of compound 29 was identified as containing the most energetically favorable binding pose. Nevertheless, both poses may occur naturally. The “one-trajectory” MM-GBSA approximation was used, which assumes that the protein and ligand have identical structures in the complex and free (dissociated) states. On the other hand, it also eliminates contributions from intramolecular energies (bonded, intramolecular van der Waals, and Coulombic), which are taken into account in the “three-trajectory approximation”. However, the “three-trajectory approximation” may introduce large uncertainties in the relative affinities;[24,64,65] in the “one-trajectory approximation”, these contributions cancel out. Thus, given its limitations and advantages and analogous to our previous studies,[61,62,66−69] the “one-trajectory approximation” is useful in evaluating the difference of binding affinities between different poses of the same ligand-receptor and to identify the energetically most favorable binding pose.

The simulation containing the most energetically favorable binding pose was extended for a total of 50 ns to refine the interactions occurring between C3c and compound 29. The average RMSD of the compound, on the basis of the initial docked structure, across the entire simulation, is 6.13 Å, and the average RMSD with respect to the average structure of the compound excluding the initial 10 ns is 0.83 Å. The former value denotes that MD simulations enable the compound to adapt its binding and optimize its interactions with C3, whereas the latter value denotes that its binding is stable within the simulation trajectory. According to the 50 ns trajectory simulation, compound 29 maintains many polar and hydrophobic interactions to C3c residues that have been shown to be key to compstatin binding[24,38] (Figure C,D). The compound maintains at least one interaction to each of the four sectors of the compstatin binding site and polar interactions to three sectors (sectors I, III, and VI). Compound 29 maintains polar interactions to the side chain of R459, the side chain of R491, and the backbone oxygen of M457. In addition, the compound maintains a polar interaction to E462, which helps anchor the compound inside the binding pocket. Polar interactions between the compound and C3c are shown in panel A. Interactions between compound 29 and sectors III and IV of the compstatin binding site are overwhelmingly governed by polar interactions (Figure A). In contrast, interactions between compound 29 and sectors I and II are governed mainly by hydrophobic interactions (Figure B). A cation−π interaction between the compound and R456 expands a favorable hydrophobic region inside the compstatin pocket to include Y433, a hydrophobic residue on a sector outside of the compstatin binding pocket. This new sector is shown in green in Figure .

Figure 8

Molecular graphics images of the C3c–compstatin binding site in complex with compound 29 (A, B) and compstatin (C, D). The images of C3c in complex with compound 29 are taken from the explicit-solvent MD simulations. The images of C3c in complex with compstatin are taken from the crystal structure of C3c in complex with compstatin (PDB: 2QKI(22)). The compstatin binding sectors VI, III, I, and II are shown in gray, yellow, blue, and orange, respectively. The green sector represents a region that is not part of the compstatin binding site but which interacts with compound 29. Important hydrogen bond interactions that occur in both the C3c/compound 29 and C3c/compstatin complexes are indicated using black dotted lines (A, C). The nonpolar binding pockets formed by C3c in complex with compound 29 and compstatin are shown in surface representation.

Comparison to a Recent Virtual Screening Study

A recent study utilized virtual screening and a variety of binding, inhibition, competitive inhibition, and functional assay approaches to identify ligands binding at the C3b binding sites of factor B, Staphylococcus aureus complement inhibitor SCIN, and compstatin, with potential complement inhibitory activities.[51] This study identified seven binding ligands, but only one of them, called cmp-5, exhibited a rather weak complement inhibitory activity. Interestingly, although cmp-5 and its structural analogues bind on the C3-binding site of compstatin, its mechanism of inhibition differs from that of compstatin. Instead of inhibiting the cleavage of C3 to C3a and C3b by the C3 convertases, cmp-5 inhibits the cleavage of C5 to C5a and C5b by the C5 convertases. The authors of this study suggested that cmp-5 influences the interaction of C3b, when part of C5 convertases, with C5, thus affecting MAC, but not C3b or C4b, deposition. Altogether, this and our study demonstrate the difficulty, but potential promise, for use of virtual screening approaches, on the basis of the large amount of structural data, to identify low-molecular mass ligands with potential to inhibit the complement system activation.

Conclusions

In this manuscript, we provide a library of 274 molecules that were identified by pharmacophore-based virtual screening and were computationally predicted to have binding affinities for the compstatin binding site of C3c. These molecules comply with pharmacophore properties that are derived from the binding characteristics of compstatin with C3c and are expected to follow a compstatin-like mechanism of binding to C3. A subset of 64 chemical compounds from these 274 molecules was experimentally tested using hemolytic assays and although they were not found to have inhibitory effects, compound 29 demonstrated binding to C3c, thus establishing a proof-of-concept for this methodology.

Despite complement being implicated in various autoimmune and inflammatory diseases, very few complement-targeted therapeutics are currently in the clinic, such as Cinryze (Shire Plc)[70] and Soliris (Alexion Pharmaceuticals, Inc.).[71] Both of them are protein-based therapeutics, and they are two of the most expensive drugs in the market, an aspect that makes them cost-prohibitive for most patients. Protein-based therapeutics are typically administered intravenously and often suffer from low stability and bioavailability and high production costs. On the other hand, low-molecular-mass chemical compound therapeutics are typically orally available, have better stability and bioavailability properties, and are less costly for industrial production. Although the virtual screening chemical compounds selected for experimental evaluation in this study did not demonstrate complement inhibitory activity, it is possible that several of them retain binding capabilities. Future binding and inhibitory activity studies are needed to further explore the entirety of the library of 274 virtual screening molecules. If binding can be experimentally observed in future studies, chemical combination of two or more C3-binding compounds may be sufficient to reproduce the binding contacts of the long and bulky compstatin family peptides and produce a molecule with compstatin-like complement inhibitory activity.