Identification of new IL-7Rα small-molecule agonists: a multi-computational approach

ABSTRACT

Interleukin 7 (IL-7) is an essential cytokine that acts as a potent growth factor of T-cells and supports the growth of B-cell precursors. IL-7 binds to a heterodimeric receptor consisting of an IL-7 receptor alpha (IL-7Rα) and the common gamma chain receptor (γc) which is shared with IL-2, IL-4, IL-9, IL-15 and IL-21. The discovery of small-molecule agonists of cytokines would be of great pharmaceutical interest with the increasing scientific rationale. In this study, a series of molecular modelling methods, including field-based pharmacophore virtual screening, protein-protein docking and molecular dynamics simulations, led to the identification of two compounds (i.e. 1 and 2) of different classes that exhibit enhanced agonistic effects by activating the IL-7 signalling cascade. One of these compounds was selected as a hit and represents the first small-molecule agonist of IL-7Rα with single-digit micromolar activity. Moreover, the prediction model of the active compound to the IL-7Rα/γc interaction complex provides insight into the binding of a small-molecule agonist to its receptor.

KEYWORDS:

• IL-7
• IL-7Rα
• agonist
• virtual screening
• protein-protein docking
• molecular dynamics

Introduction

Interleukin 7 (IL-7) is a haematopoietic growth factor protein that is secreted by stromal cells in the thymus and bone marrow [1–3]. IL-7 is one of the essential cytokines for the homoeostatic maintenance and development of B and T lymphocytes [4,5]. Its signalling cascade starts with IL-7 binding to the IL-7 receptor, which is composed of common gamma chain receptor (γc) and IL-7 receptor alpha (IL-7Rα) [6,7]. Binding results in a cascade of signals, including the expression of multiple transcription factors (Jak3, STAT, PI3K and SRC) which are crucial for T-cell proliferation and activation [8,9]. IL-7 is considered one of the most important immunotherapy agents for treating various malignancies and viral infections including human immunodeficiency virus (HIV) [10–12]. Previous studies have shown that IL-7 signalling promotes osteoclastogenesis and fibroblast activation, which are involved in tissue damage in chronic inflammation [13]. IL-7 signalling also regulates antigen-reactive T cells, which suggests that abnormal IL-7 signalling can be related to the pathogenesis of autoimmune diseases, such as type 1 diabetes, rheumatoid arthritis, multiple sclerosis and ulcerative colitis [14]. IL-7 signalling also induces the production of TNF-α, IL-1β and IL-1α by monocytes, and these cytokines inhibit melanoma growth [5]. In addition, IL-7 signalling induces an increase in CD8 + T lymphocytes, which can attack cancerous cells.

Recent studies have suggested that the role of IL-7 in immunotherapy is not limited to treating viruses but can also be applied to treating cancers [15,16]. A human clinical study with cancer patients suggested that administration of IL-7 can transiently lead to an expansion of CD8+ and CD4 + T cell populations along with a decrease in CD4+ regulatory T cell populations [17]. This modulation of T cell homoeostasis can exert potential antitumor effects and has been intensely studied recently [18,19].

However, due to the high cost of manufacturing recombinant IL-7 proteins [20] and the limitations on delivery route (mainly intravenous) [21,22] searching for small-molecule agonists showing similar immunotherapeutic effects via activation of the IL-7 signalling cascade is worth undertaking. Here, we present our efforts in a small-molecule search for potential IL-7 agonists via intensive virtual screening and molecular assays.

Materials and methods

In silico molecular modelling studies

Shape- and electrostatic field- based virtual screening

IL-7Rα in complex with IL-7 was available, and the X-ray crystal structure of the IL-7Rα/IL-7 complex was obtained from the Protein Data Bank (PDB code 3DI3) [23]. The IL-7 protein was used to generate a pharmacophore model using ligand-based virtual screening within OpenEye Scientific software [24]. A pharmacophore query model was created from the binding site of IL-7 with IL-7Rα based on the binding motif and consists of a helix; its residues are amino acids (aa) 72–91. The 557,360-membered KCB library database was prepared using OMEGA v3.0.1, thereby generating a maximum number of 200 low-energy conformers for each molecule. The 3D shape similarity search of the helix as a query was performed using ROCS v3.2.2, and the top 10,000 hits were selected according to the ROCS_TanimotoCombo score. We then performed a 3D electrostatic potential field map similarity search for the screened 10,000 compounds using EON v2.2.0, and the top 1,000 hits were selected by the EON_ET_Combo score. A total of 182 compounds were ultimately selected by eye selection for the in vitro assay.

Similarity-based virtual screening using two-dimensional fingerprints

To identify more analogues of active compound 1, we performed a second round of virtual screening using the fingerprint similarity method. The similarity between single pairs of compounds was calculated using Pipeline Pilot 2020 [25]. The KCB library database was screened using a Tanimoto coefficient ≥ 0.8 compared with active compound 1, and 9 analogues were ultimately selected by eye selection for the in vitro assay.

Protein-protein docking

To understand the IL-7Rα/γc heterodimeric complex in the presence of compound 1, we performed a rigid-body protein-protein docking simulation using the ZDOCK protocol in Discovery Studio 2020 [26]. The structure of γc extracted from the structure of the IL-2Rβ/γc complex (PDB code 5M5E) [27] was used. The receptor-ligand complex model for IL-7Rα bound to compound 1 was retrieved by superimposition of compound 1 onto IL-7, which was a pharmacophore query model for ligand-based virtual screening. The IL-7Rα/compound 1 and γc protein structures were assigned the CHARMM force field for protein-protein docking procedures. The top poses of docking models between IL-7Rα/compound 1 and γc were chosen by visual inspection with the ZDOCK score from the top 500 poses. Five protein-protein docking models were selected by high scoring and eye selection.

Molecular dynamics simulation

For further refinement of IL-7Rα/γc docking models and the molecular interaction analysis of compound 1, molecular dynamics (MD) simulations were performed using Desmond v6.4 with the OPLS3e force field in Schrodinger Suite 2020–4 [28]. System Builder was used for solvation, employing predefined TIP3P water molecules in an orthorhombic box with dimensions of 15 Å x 15 Å x 15 Å, and the overall complex was neutralized by adding Cl- counterions at least 25 Å from ligand compound 1. The NaCl concentration was 0.15 mol/L. MD simulations 500 ns in length were produced under periodic boundary conditions in the NPT ensemble at normal temperature (300 K) and pressure (1.01325 bar). Recording intervals of 1.2 and 500 ps were used for energy calculations and trajectory analyses. The final proposed IL-7Rα (bound to compound 1)/γc docking model for the equilibrium state of the MD simulations was represented using Discovery Studio 2020.

Biological assays

Materials

The PathHunter IL-7 dimerization assay kit was purchased from Eurofins DiscoveRX (Luxembourg). Human interleukin 7 protein (IL-7) was obtained from Sino Biological Inc. (Beijing, China).

IL-7 agonism assay

PathHunter-expressing dimerized cells (DLD-1-based) were prepared as described in the manufacturer’s protocol. In brief, cells were thawed and seeded on 96-well plates at a density of 10,000 cells/well and maintained at 37°C in a humidified incubator with 5% CO₂ for 48 hours in 100 µL of cell culture medium. Cells were then treated with 10 µL of compounds at the designated concentrations (1, 10, 100 and 1000 μM) and IL-7 in an incubator at 37°C. After 24 hours of incubation, 110 µL of detection reagents was added, and the plates were further incubated at room temperature. After 1 hour of incubation, chemiluminescence signals were monitored with a spectrophotometer [29].

Statistical analysis

Statistical analysis was performed using OriginPro 8.5 [30]. The half maximal effective concentration (EC₅₀) of compound 1 and IL-7 was estimated using a sigmoidal logistic function, type 1, in Origin 8.5.

$y=\frac{a}{1+e^{-k(x-x_c)}}$

Each experiment was performed in triplicate, and values are expressed as the mean ± standard deviation (SD).

Results and discussion

Ligand-based virtual screening by shape and electrostatic pharmacophore

To discover novel small-molecule agonists that promote IL-7Rα/γc heteodimerization, we explored the IL-7 helix domain region where IL-7Rα interacts to the three-dimensional X-ray crystal structure of the IL-7Rα/IL-7 complex (PDB code 3DI3) was utilized to generate a shape- and field-based pharmacophore model based on aa 72–91 that are involved in the binding interface of IL-7 to IL-7Rα (Figure 1) for the similarity search. Ligand-based virtual screening was performed using OpenEye Scientific software. The Korea Chemical Bank (KCB) library, which consists of 557,360 compounds, was prepared for the 3D similarity search using OMEGA v3.0.1, which was used to generate up to 200 low-energy conformers for each molecule. For virtual screening using a large chemical database, two computational strategies were employed. First, a 3D shape similarity search of peptide molecules consisting of aa 72–91 in IL-7 as a query was performed using ROCS v3.2.2, and the top 10,000 hits for query molecules were selected based on the ROCS_TanimotoCombo score. Second, a 3D electrostatic property similarity search was executed for the 10,000 selected compounds using EON v2.2.0. After the positive and negative field contour maps with EON_ET_Combo scores were analysed, the top 500 hits were selected. An additional visual inspection was carried out to remove the PAINS [31]. As a result, 182 compounds were ultimately selected to examine their IL-7R agonist effects in vitro.

Figure 1. Ligand-based virtual screening using 3D shape and electrostatic fields. The IL-7Rα binding motif, a helix consisting of aa 72–91 from the IL-7 crystal structure (PDB code 3DI3), was utilized as a template for 3D similarity screening. Representation of electrostatic potential maps for the helix consisting of aa 72–91. The positive contours are coloured blue, and the negative contours are coloured red. A small-molecule compound obtained from EON is coloured yellow

In vitro assay for IL-7 agonistic activity

The IL-7 agonistic activities of the 182 selected compounds were measured and compound 1 was most potent and promoted for further biological activity study. We used a commercially available IL-7Rα/γc dimerization assay kit and detected compound/ligand-induced dimerization of two subunits of IL-7Rα and γc. Engineered DLD-1 cells coexpress one receptor subunit fused to an enzyme donor and a dimerizing partner receptor fused to an enzyme acceptor, and binding of an agonist to the receptor subunit induces dimerization with the partner receptor, forcing complementation of the two enzyme fragments. This complementation results in the formation of a functional enzyme, which can generate a chemiluminescent signal that can be monitored. The EC₅₀ of compound 1 (Figure 2a) was 9.86 μM (Figure 2b), while the EC₅₀ of IL-7 was 0.71 nM (Figure 2c). The other hit (compound 2 in Figure S1) was derived from virtual screening; this compound from a structurally different class than that of compound 1 also induced dimerization of IL-7Rα and γc, although it showed less potency (EC₅₀ = 48.8 μM) than that of compound 1. The detailed structure and biochemical activity of compound 2 are presented in Figure S1. To this end, compound 1 from the primary virtual screening was selected for further studies.

Figure 2. (a) Compound 1 structure, (b) the biochemical activity of compound 1 and (c) the biochemical activity of IL-7

Similarity-based virtual screening using 2D fingerprints

To validate compound 1 as a primary hit, we performed a second round of virtual screening based on two-dimensional fingerprint similarity. The KCB library was screened based on a Tanimoto coefficient ≥ 0.8 compared with active compound 1 using Pipeline Pilot. Nine analogues were chosen by eye selection for further IL-7 agonistic activity measurements. However, none of the newly selected compounds 183–191 were as effective as compound 1, and their agonist activities are summarized in Table 1.

Table 1. Structure-activity relationships of IL-7 agonist activities

Compound	R^1	R^2	R^3	R^4	Conc. (µM)	Activity (%)^a
1	F	CH3	H	CH3	10	56.2
183	F	H	H	CH3	10	18.1
184	F	F	H	CH3	10	−19.5
185	F	Cl	H	CH3	10	−15.1
186	H	H	CH3	CH3	10	17.8
187	H	CH3	H	CH3	10	21.2
188	H	OCF3	H	CH3	10	−8.8
189	H	H	H	CH3	10	−16.2
190	F	H	H	H	10	−14.2
191	F	H	H	CH2CH2OH	10	42.2

^a10 µM single concentration for each compounds was used for screening.

3D shape and electrostatic similarity analysis of compound 1

The Electrostatic Tanimoto combo (ET_combo) scores for the 182 virtual hits ranged from 0.158 to 0.278, while the value for compound 1 was 0.240 (the electrostatic field and shape similarity scores were 0.092 and 0.148, respectively), and the compound ranked in the top 10 among the 182 hits. The electrostatic potential maps of aa 72–91 in IL-7 and compound 1 are depicted in Figure 3b. The potential maps revealed distinct differences in the electrostatic surfaces between the two molecules. Regarding the similarity to the α-helix consisting of aa 72–91 in IL-7 as a query, compound 1 is radically different with respect to its 2D and 3D structures but similar in the 3D shape to a Shape Tanimoto (ST) similarity value of 0.148. Although no part in compound 1 was aligned to the contour map with aa 72–91 of IL-7, the sequence adjacent to the piperazine group was comparable (Figure 3b). The other three parts were not electrostatically similar to aa 72–91 of IL-7 but could overlap well by shape (Figure 3a). Compound 1 had an ST value of 0.148, which was higher than 0.092 with the Poisson-Boltzmann Electrostatic Tanimoto (ET_pb) value. Overall, our field-based pharmacophore virtual screening results suggest that compound 1 and the helix consisting of aa 72–91 of IL-7 may be structurally dissimilar and therefore bind to different binding sites. Full details for 182 ligand-based virtually screened compounds, including their structures, EON ranks, ET_combo scores, ST scores, ET_pb scores and primary screening activity data, are provided in Supplementary Table 1.

Figure 3. Three-dimensional shape and electrostatic contour map comparison between the α-helix consisting of aa 72–91 of IL-7 crystal structure and compound 1. (a) overlay of the aa 72–91 helix (red stick model) and compound 1 (yellow ball and stick model) and (b) structures of the aa 72–91 helix and compound 1 with electrostatic contour maps

Protein-protein docking and molecular dynamics simulation for the prediction of the IL-7Rα/compound 1/γc heterodimer

To explore the structural basis of IL-7Rα/γc heterodimerization, we performed protein-protein docking and MD simulations. The X-ray crystal structure of the IL-7Rα/γc complex was not yet available at the time of this work. As a surrogate model, we used the γc protein structure from the IL-2Rβ/γc complex structure (PDB code 5M5E). According to our previous ligand-based pharmacophore model and the release of the IL-7Rα/IL-7 complex (PDB code 3DI3), we hypothesized that compound 1 binds to the IL-7Rα/IL-7 binding interface. The initial binding mode of compound 1 to IL-7Rα was retrieved by superimposition of compound 1 onto IL-7. Then, a protein-protein docking simulation was run to predict the binding model of the γc protein with IL-7Rα/compound 1 by employing the ZDOCK protocol in Discovery Studio. As shown in Figure 4a, γc binds to IL-7Rα complexed with compound 1 in our protein-protein docking study based on high scoring selection. This was supported by relaxation of the docking poses conducted with MD simulation using Desmond in Schrodinger 2020–4. The binding mode of IL-7Rα with γc remained similar in 500 ns MD simulations relative to the initial structure. Interestingly, compound 1 binds deeper into the IL-7Rα/γc binding interface and closer to γc (Figure 4b). Further analysis of the structural changes of IL-7Rα/γc in the presence of bound IL-7 revealed no significant changes in the dimeric conformations between the X-ray crystal structure and the structure determined by MD simulation. The helix consisting of aa 72–91 (Asp74, Leu77, Leu80, Lys81 and Ile88) of IL-7 forms hydrophobic interactions with Val58, Leu80 and Ile82 and hydrogen bonds with Ser31 and Lys77 in the binding site of IL-7Rα (Figure S2b).

Figure 4. Molecular dynamics simulations of IL-7Rα/compound 1/γc complex models. (a) the protein-protein docking model of IL-7Rα bound compound 1, the pose of which is the superimposition of the 3D conformation by ligand-based virtual screening and γc, and (b) the equilibrium state complex model after 500 ns MD simulations. IL-7Rα subunit: blue ribbon model; compound 1: yellow ball and stick model; γc subunit: green ribbon model

Molecular interaction analysis of compound 1

Through the application of MD simulations, the insight gained through protein-protein docking was further extended [32]. We performed an in-depth analysis of the binding interactions formed between compound 1 and the IL-7Rα/γc complex for molecular analysis of ligand recognition, in which the binding mode of compound 1 was also obtained from the corresponding MD model. As shown in Figure 5, the 3,4‐dihydro‐1,3‐pyridooxazine ring of compound 1 forms significant π-π stacking interactions with IL-7Rα around protein residue Tyr192 and alkyl hydrophobic interactions with IL-7Rα around protein residue Phe193. Notably, the N4-nitrogen atom of oxazine forms hydrogen bonds with the hydroxy group of the Tyr182 residue from γc. The triazolopyrazine ring of compound 1, located between the two proteins, shows very strong π-π stacking interactions with Tyr139 of IL-7Rα, in addition to other alkyl hydrophobic interactions with Pro207 of γc (see the protein-ligand interaction diagram in Figure 5b). Moreover, the methyl group of the 1-fluoro-3-methylbenzene ring in compound 1 also forms alkyl hydrophobic interactions with Ile82 of IL-7Rα, and the 1-fluoro-3-methylbenzene ring makes ca-π interactions with Leu208 of γc. The terminal methylpiperazine substituent is oriented towards γc by forming ca-π interactions with Tyr103. Taken together, our molecular modelling studies suggest that the important residues for small-molecule agonist binding promote IL-7Rα/γc dimerization throughout MD simulation.

Figure 5. Proposed binding model and interaction profile analysis of compound 1 for IL-7Rα/γc proteins. (a) predicted binding mode of compound 1 (yellow ball and stick model) after 500 ns MD simulations. The IL-7Rα and γc proteins are displayed as blue and green ribbon models, respectively. The hydrogen bonds are shown as green dashed lines, and hydrophobic interactions are represented by pink dashed lines. (b) contact histograms of IL-7Rα (B chain)/γc (C chain) proteins with compound 1 are represented as a coloured stacked bar plot between 200 and 500 ns MD simulations. Hydrogen bonding is displayed by yellow-green coloured bars, and hydrophobic interactions are displayed by lavender coloured bars for clarity

Conclusion

To date, to the best of our knowledge, no IL-7Rα agonist with low molecular weight has been described in the literature. In this study, we report the results of virtual screening with the aim of identifying agonists capable of promoting the IL-7Rα/γc complex. This is the first attempt to target IL-7Rα by an agonist approach, and this study proposes an in silico approach for the discovery of new chemical entities for IL-7. A combined computational approach, including virtual screening and protein-protein docking and MD simulations, has been applied to identify small-molecule agonists other than IL-7 from KCB chemical library databases. A total set of 191 compounds (182 first-round virtual hits by ligand-based virtual screening and 9 second-round virtual hits by a 2D fingerprint similarity search with compound 1) were tested to validate our results in vitro, of which two compounds were shown to have agonistic effects, promoting IL-7Rα/γc dimerization. Compound 1, as a first-round virtual screening hit, exerts single-digit micromolar efficiency on heterodimerization by binding to IL-7Rα/γc. Following an in silico characterization of the IL-7Rα/γc complex interface, a protein-protein docking simulation was performed. In addition, the use of an MD simulation method enabled the discovery of a small-molecule agonist and its binding-determining structural features against the IL-7Rα/γc complex, which will aid in the further development of agonists. In the present research work, a series of compounds with similar functional groups were tested in vitro, and the effects of substituents on compound 1 were determined. The strategy described herein offers a feasible and effective approach to discover other bioactive compounds using in silico virtual screening. In this study, we identified the first class of small-molecule agonists of the IL-7Rα/γc complex, and we are currently searching for more analogues capable of higher-level IL-7Rα/γc dimerization.

Acknowledgement

This work was supported by Korea Research Institute of Chemical Technology intramural funding (grant numbers KK2062–132 and SI2131–50). The chemical library used in this study was provided by Korea Chemical Bank and by Hee Jung Jung Ph.D. of KRICT for donating the compound.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed at: https://doi.org/10.1080/1062936X.2021.1969684.

Additional information

Funding

This work was supported by Korea Research Institute of Chemical Technology intramural funding [grant numbers KK2062–132 and SI2131–50].