Identification of Potential Small Molecule Allosteric Modulator Sites on IL-1R1 Ectodomain Using Accelerated Conformational Sampling Method

The interleukin-1 receptor (IL-1R) is the founding member of the interleukin 1 receptor family which activates innate immune response by its binding to cytokines. Reports showed dysregulation of cytokine production leads to aberrant immune cells activation which contributes to auto-inflammatory disorders and diseases. Current therapeutic strategies focus on utilizing antibodies or chimeric cytokine biologics. The large protein-protein interaction interface between cytokine receptor and cytokine poses a challenge in identifying binding sites for small molecule inhibitor development. Based on the significant conformational change of IL-1R type 1 (IL-1R1) ectodomain upon binding to different ligands observed in crystal structures, we hypothesized that transient small molecule binding sites may exist when IL-1R1 undergoes conformational transition and thus suitable for inhibitor development. Here, we employed accelerated molecular dynamics (MD) simulation to efficiently sample conformational space of IL-1R1 ectodomain. Representative IL-1R1 ectodomain conformations determined from the hierarchy cluster analysis were analyzed by the SiteMap program which leads to identify small molecule binding sites at the protein-protein interaction interface and allosteric modulator locations. The cosolvent mapping analysis using phenol as the probe molecule further confirms the allosteric modulator site as a binding hotspot. Eight highest ranked fragment molecules identified from in silico screening at the modulator site were evaluated by MD simulations. Four of them restricted the IL-1R1 dynamical motion to inactive conformational space. The strategy from this study, subject to in vitro experimental validation, can be useful to identify small molecule compounds targeting the allosteric modulator sites of IL-1R and prevent IL-1R from binding to cytokine by trapping IL-1R in inactive conformations.


Introduction
The interleukin-1 (IL-1) family of ligands and receptors are important regulators of innate inflammatory and immune responses [1]. Currently, eleven IL-1 ligands (or cytokines) and receptors were identified respectively [2]. These IL-1 ligands function as agonists and antagonists to IL-1 signaling by binding to IL-1 receptors including ten single pass transmembrane proteins and IL-18BP which lacks transmembrane domain [3]. Regulation of inflammatory signaling via IL-1 family of receptors can be mediated either by its binding to the activating cytokines or antagonistic ligands such as IL-1Ra which inhibit downstream signaling [4]. Some members of the IL-1 receptors are decoy receptors such as IL-1 receptor type 2 (IL-1R2) which lacks the cytoplasmic element to recruit downstream effector proteins. Additional level of regulation is controlled via the cleavage of the ectodomain of IL-1 receptors such as that found in IL-1 receptor type 1 (denoted sIL-1R1) [5] or the alternatively spliced soluble Suppressor of Tumorigenicity 2 (sST2) protein [6]. sIL-1R1 and sST2 have only the ectodomain of their respective receptors, IL-1R1 and ST2. They circulate in the extracellular milieu to modulate their respective binding cytokines concentrations and attenuate the inflammatory responses mediated by these ligands [5,7,8].
Roles of many members in the IL-1 family play in inflammatory diseases and auto-inflammatory diseases or disorders have been reported [9]. In one case, loss of function homozygous mutation in the endogenous IL-1 receptor antagonist (IL-1Ra) has been identified in infants with early death due to overwhelming systemic inflammation found in skins, joints and bones [10]. Contribution of the IL-1 family of proteins to cardiovascular diseases, Type-2 diabetes, rheumatoid arthritis and chronic inflammatory diseases have also been reported [9]. Current available therapy for treating IL-1 family related diseases is limited to Canakinumab, Gevokizumab and LY2189102 which are mainly anti-IL-1β monoclonal antibodies [9]. No small orally available drugs are available. In recent years, the expanding roles of the member proteins in IL-1 family play in multiple diseases have been discovered [2,11]. One example is that the elevated level of soluble ST2 in the plasma of patients receiving allogeneic hematopoietic stem cell transplantation predicted their graft-versus-host disease (GVHD) mortality [12,13]. The prognostic value of soluble ST2 was attributed to its role as a decoy receptor of IL-33. Binding of sST2 and IL-33 decreases the level of unbound IL-33 in plasma which drives type-2 T cells towards a type-1 T cell response during transplantation. Thus, there are needs to discover and develop new therapy targeting this important cytokine family of proteins for treating inflammatory and immune-mediated diseases.
Currently, FDA approved specific IL-1 targeting therapies are monoclonal antibodies. Orally available small molecule drugs have advantages in several aspects including ease of administration, cost effectiveness and flexible dosage schedule design. Given the chronic nature of the inflammatory disorders and diseases associated with IL-1 family of proteins [14], efficacious orally available drugs will benefit patients in their long-term treatment. To develop small molecule inhibitors for the IL-1 family of proteins, we focus on targeting the ectodomain of IL-1 receptor type 1 (IL-1R1) protein, the first cytokine receptor in the family discovered [15], in this work. Binding of IL-1R1 with IL-1α or IL-1β leads to subsequent recruitment of IL-1RAcP at the cell membrane [16] and downstream cytosolic NF-κB activation. A small molecule probe (NP1) was reported recently to inhibit the association between IL-1R and TLR4 and the downstream NF-κB activation [17]. However, it was suggested to bind to the cytosolic component of the cytokine receptor protein complex without clear indication of the protein target.
Available crystal structures of IL-1R1 indicated that IL-1R1 binds with IL-1β using three domains (D1,D2 and D3) and their contact interface (4180.5 Å 2 ) is much larger than those in typical small molecule binding sites (1500-3000 Å 2 ) [18]. While IL-1R1 binds with IL-1Ra involving small conformational changes of the D3 domain [19], significant changes were found when IL-1R1 binds with a peptide antagonist (AF10847) [20]. We hypothesized that the flexibility of IL-1R1 ectodomain observed from crystal structures may indicate the existences of multiple conformations that can potentially be exploited by small molecule inhibitors development. To investigate this inhibitor development strategy, we employed computational methods to first sample extensively IL-1R1 conformations based on accelerated molecular dynamics simulation algorithm [21]. Selected represented conformations determined by the hierarchy cluster method were analyzed by the Sitemap program [22] to identify "druggable" small molecule binding sites. Potential small molecule binding sites, denoted as allosteric modulator sites, not at the cytokine receptor/cytokine interaction interface were determined. The druggable binding sites were further assessed as binding hotspots via our recently developed cosolvent mapping method [23][24][25] which incorporated binding site flexibility. An in silico screening of a commercial fragment library was then used to target one modulator site in an unreported inactive IL-1R1 conformation. From the eight top-ranked compounds, we found four compounds bound to the allosteric modulator site for extended period of time and restricted the IL-1R1 dynamical motion to an inactive conformational space during the 8 ns MD simulations. The approach described here can be an attractive strategy to discover small molecule inhibitors targeting the challenging cytokine receptor proteins.

MD simulation setup
In preparation of the simulations, we have used the MOE program [30] to determine the protonation state of ionizable groups on IL-1R1 under the standard physiological condition. PMEMD from Amber (version 12) [31] was used for molecular dynamics simulations. The Amber 99SB force field parameters [32] were used for the amino acids. To prepare the topology and coordinate files, counter ions were added to neutralize the charges in IL-1R1 before it was placed in a 13Å octahedral box of water. The TIP3P [33] water model was used. A 3000-step minimization (steps 1-1000 using conjugated gradient followed by 2000 steps steepest decent) was first carried out. After minimization, a 500 ps constant volume and constant temperature (NVT) simulation was performed to raise the temperature of the system to 298K while constraining backbone atoms with a 5 kcal/mol/Å 2 force constant with reference to the crystal structure. A second 200 ps constant pressure and constant temperature (NPT) simulation at 298 K was performed while constraining backbone atoms with a 2 kcal/mol/Å 2 force constant with reference to the crystal structure. The system was then equilibrated for 1 ns at 298K without any constraints followed by the production run. All the MD simulations were in the isobaric isothermal (NPT, T = 298K and P = 1 atm) ensemble. The SHAKE [34] algorithm was used to fix bonds involving hydrogen. The PME method [35] was used and the non-bonded cutoff distance was set at 10Å. The time step was 2 fs, and neighboring pairs list was updated every 20 steps.
For performing accelerated MD simulations, a 2 ns conventional MD (no modification of the potential energy function) was conducted to determine the average values of the potential energy for the total system (V total ) and the dihedral angles energy of the proteins (V dih ). The threshold potential energy was defined as V total + N atoms /5 (N atoms : total number of atoms in the system). Three threshold potential energies of the dihedral angle motion were used called AMD1, AMD2 and AMD3 with V dih + (number of residues) times 3.5, 4.2 and 4.9 respectively. These parameters were taken as suggested by previous works [21]. The V total and V dih parameters used in each machine were based on the averaged potential energy values calculated respectively. A total simulations time of 118 ns were collected from 2 ns equilibrium production run, 40 ns run of AMD1 parameters, 10 ns run of AMD2 parameters, 10 ns run of AMD3 parameters performed at the GORDON cluster supported by XSEDE, 30 ns run of AMD1 parameters in a local GPU machine and 26 ns run of AMD1 in an eight cores local machine. The lengths of the simulations performed at the GORDON cluster were limited by the computation times allocated to this work by XSEDE. Additional simulation runs performed at the local machines were included to ensure the coverage of the conformational space guided by the projections of conformations to the principal component spaces discussed later in the Methods section.
The procedures of performing cosolvent mapping simulation and analysis have been reported previously [23][24][25]. In this work, we used phenol as the probe molecule because the phenol group is frequently found in the fragment screening against proteins involved in proteinprotein interaction [36]. Cosolvent mapping analyses were done using the trajectories of 4 ns cosolvent MD simulations.
Force field parameters of the small molecules were derived using the Antechamber module in the Amber program suite. The protocol for generating the point charge parameters is as follows: The docked pose of each molecule was minimized at the RHF level using a 6-31GÃ basis set with Gaussian09 [37]. The electrostatic field potential calculated from Gaussian09 was used to derive the point charges at each atom site based on the RESP fitting procedure [38].

Principal component analysis and hierarchy cluster analysis of the conformations
Five crystal structures of IL-1R1 bound to different ligands were used in the principal component analysis (PCA) implemented in the Bio3D version 2.0 [39]. In the PCA, amino acid sequences of the crystal structures were aligned and gaps detected based on MUSCLE [40]. Cα atoms coordinates of non-gap amino acids in the sequence alignments were first superimposed. Those with largest displacements were iteratively removed to reach a core set of amino acids to represent the structurally "static" amino acids among the variations of conformations analyzed. Based on the alignment of the core set of amino acids, the averaged displacements of all nongap amino acids were calculated to construct the covariance matrix. Diagonalization of the matrix gives eigenvectors and eigenvalues. The eigenvectors are referred as principal components (PCs) and the coefficients of the eigenvectors characterize the collective atomic displacements for each principle component. The eigenvalues of the matrix gives the contribution of each eigenvector to the overall displacements in the covariance matrix. Such an analysis can be considered as the linear transformation of the coordinate system from Cartesian bases to collective movement bases. The PC gave the physical representation of collection motion in the protein dynamics because the Cα atoms were connected (or correlated) via backbone peptide bonds. Detailed procedures can be found in the Bio3d website. Snapshots of conformations obtained from each MD simulation were then first aligned to the initial conformation and the conformational changes in the simulated trajectory were mapped to the associated principal components shown in figures. The PCA using the five crystal structures showed the first three principal components contributed 94, 4.6 and 0.7% to the overall structural displacements. Mapping of conformations from the trajectory to two principal components informed a two dimensional subspace characterization of the conformational changes and/or transition.
For hierarchy cluster analysis, a second PCA of 1180 aMD generated IL-1R1 conformations taken at every 100 ps from the 118 ns trajectory were first performed without including the crystal structures. For this hierarchy cluster analysis, we used the first 38 PCs, which accounted for 90% of backbone structural variation of the 1180 IL-1R1 conformations. In the cluster analysis, the Euclidean distances between the 38 PCs were used to construct the dendrogram as implemented in the hclust function of the R program [41]. By using a cutoff value of 40 Å at the dendrogram, we obtained 25 cluster groups. The conformation closest to the center of each cluster group was used to represent each cluster group, also called representative conformation.

Sitemap Analysis
Three crystal structures (chain A from 4GAF [19], chain R from 1IRA [28], chain B from 1G0Y [20]) and twenty five representative conformations of IL-1R1 determined from the hierarchy cluster analysis were used in the Sitemap program [22] to identify potential small molecule binding sites in each conformation. Each conformation was first processed in the Glide program [42] from the Maestro 9.7 program suite (Release 2014-1) [42] followed by the Sitemap analysis using default parameters. Up to five binding sites were saved. Among the twenty five MD generated conformations, only conformation 399 gave four binding sites. Druggability index of each binding site was based on the Dscore value in the Sitemap program.

In silico screening and the fragment library
Conformation 289 was used in the in silico fragment library screening. The virtual screening workflow procedure in the Glide program was adopted in the screening calculation. Val117 of IL-1R1 (P2 site) was selected as the center of the docking site and the inner and outer rectangular box sizes were set to 10 and 30 Å respectively in the grid generation. The Maybridge fragment library satisfied the rule of three [43] were obtained from the Maybridge website with a total of 2783 compounds. Structures of the Maybridge fragment library were prepared using the Ligprep program in the Maestro program suite using the default setting. Two compounds (both adamantine derivatives) failed in the Ligprep program and were fixed manually. A total of 4409 compound conformations were obtained and used in the in silico screening calculations. The poses of the top eight ranked compounds based on the XP docking function from the in silico screening were inspected and used in the subsequent MD simulations and analyses.

Results and Discussions Crystal structures between IL-1R1 and different ligands
Three different conformations of IL-1R1 ectodomain have been determined via X-ray crystallography including the complexes of IL-1R1/IL-1β, IL-1R1/IL-1Ra and IL-1R1/AF10847 (Fig. 1). While IL-1β is a cytokine that activates IL-1R, IL-1Ra is an IL-1 receptor antagonist and AF10847 is a peptide antagonist [20,44]. The crystal structure of IL-1R1/IL-1β ( Fig. 1A) showed that three domains of IL-1R1 (D1-D3 domains) interact with IL-1β using extensive surface contacts. Their contact surface area (contributed from both proteins) is 4180.5 Å 2 based on the structure of PDB entry: 4DEP calculated by NACCESS [45]. The interaction interface between IL-1R1/IL-1β with IL-1RAcP is less with a surface area of 2971.1 Å 2 . As reported, the antagonistic cytokine, IL-1Ra, interacts primarily with D1 and D2 domains of IL-1R1 but much less with the D3 domain of IL-1R1 (Fig. 1B) [19,28,46]. From our analysis, the D3 domain of IL-1R1 rotates 30 degrees away from that in the IL-1β bound IL-1R1 conformation; thus, making little contacts with IL-1Ra. The rotation of the D3 domain led to less contact interface area between IL-1R1 and IL-1Ra (3367.4 Å 2 ). In both cytokine bound structures, the IL-1R1 ectodomain adopts a clamp form to bind with the spherical shape cytokine ligands that fold into a conserved 12 stranded beta sheet structure. Although the antagonistic helical peptide, AF10847, is much smaller in size (21 amino acids), a dramatic conformational change in the D3 domain of IL-1R1 facilitated via the loop between the D2 and D3 domains occurred and resulted in favorable interaction between IL-1R1 and AF10847 (Fig. 1C). Based on the structure of IL-1R1/AF10847, the contact interface area is 2505.4 Å 2 . The large interaction interfaces between IL-1R1 and different ligands also correlate with their high binding affinities. For examples, the IC50 values between IL-1Ra, AF10847 and IL-1R1 were reported to be 1.6 and 2.6 nM respectively [47]. Using a surface plasma resonance (SPR) based kinetic binding assay, the K D values for IL-1R1 and IL-1β, IL-1Ra were reported to be 2.0 and 0.33 nM respectively [19]. Comparison of these three ligand-bound IL-1R1 structures (Fig. 1D) also suggested the loop between the D2 and D3 domains is much flexible than that between the D1 and D2 domains. This flexible loop permits the variation of the relative orientation between the D1-D2 domains and D3 domain and allows recognition of IL-1R1 to different ligands.

Conformational sampling of IL-1R1 based on conventional and accelerated MD simulations
The crystal structures of IL-1R1 with IL-1β, IL-1Ra and AF10847 gave static representations of IL-1R1 conformations stabilized by each ligand molecule. To explore and sample the conformational spaces of IL-1R1 extensively, we performed conventional and accelerated MD simulations of IL-1R1 using the IL-1R1 conformation in the IL-1R1/EBI-005 crystal structure. EBI-005 is an IL-1β chimeric protein that binds to IL-1R1 ectodomain at a higher affinity than IL-1β (K D = < 0.014 versus 2.0 nM) [19]. To analyze the conformations of IL-1R1 obtained from the MD simulations, we first performed the principal component analysis (PCA) of the five available crystal structures of IL-1R1. Based on the PCA analysis, the first two principal components (PC1 and PC2) account for 98.6% of backbone conformational variations associated with the five crystal structures. As shown in Fig. 2, projections of the five IL-1R1 structures onto PC1 and PC2 subspace gave a clear distinction of the IL-1β-and IL-1Ra-bound IL-1R1 conformations from the AF10847-bound IL-1R1 conformation. Specifically, IL-1β-, IL-1Rabound IL-1R1 conformations were mapped to a similar region (PC1 = -70, PC2 = -20) whereas the AF10847-bound IL-1R1 conformation was projected to a further out region (PC1 = 330, PC2 = 4).
Although only five crystal structures of IL-1R1 in different ligand-bound conformations are available currently, the first two principal components of the PCA analysis can be informative to characterize the conformational changes of IL-1R1 in the MD simulations. Based on the mapping, conformations of IL-1R1 started with the EBI-005-bound conformation were found trapped in a region (PC1 = 140, PC2 = -80) in the 37 ns cMD run (S1 Fig.). IL-1R1 started with the AF10847-bound conformation remained confined at its initial conformational spaces during 20 ns cMD simulations (S1 Fig.). For the same 30 ns of simulation times, the IL-1Ra-bound IL-1R1 conformations overlap with part of the region visited by the EBI-005-bound IL-1R1 simulation. Results of cMD simulations indicated their limitation to overcome free energy barriers between different conformational microstates of IL-1R1 ectodomain because no overlapping conformations were found between EBI-005-bound and AF10847-bound IL-1R1 simulations.
In the aMD simulation, the potential energy wells were modified or raised to facilitate the barrier crossing of protein conformations between different microstates. Here, we used three different parameters denoted as AMD1, AMD2 and AMD3 where the modified potential wells become progressively shallower from AMD1 to AMD3. In Fig. 2B-D, we found the parameters used for AMD1 and AMD3 allowed sampling of IL-1R1 conformations to the proximity of the AF10847-bound IL-1R1 structure but not those using the AMD2 parameters. IL-1R1 conformations obtained from the GPU machine were found to map to different conformational spaces not accessed by other trajectories. Results from the aMD simulations indicated that the modifications to the potential energy wells permits extensive conformational sampling of IL-1R1 in a relatively short simulations times (118 ns in total). Projection of the IL-1R1 conformations from aMD simulations to the three dimensional subspace (PC1, PC2 and PC3) was provided in S2 Fig. which showed the extensive coverage of the conformational spaces including all ligand-bound conformations from the crystal structures. We further analyzed the conformational changes of IL-1R1 from the aMD simulations by studying backbone fluctuation of the D1, D2 and D3 domains. In all IL-1R1 conformations, the D1 and D2 domains of IL-1R1 showed very limited backbone variations where the rootmean-square deviations (RMSD) of the backbone atoms vary between 2-3 Å from the crystal structure (Fig. 3A). Larger backbone motions of IL-1R1 were observed in the D3 domain where the distribution of RMSDs peaks at 35 Å from the EBI-005-bound IL-1R1 crystal structure. The large RMSD values in the D3 domain do not result from the structural unfolding. It is attributed to the large rotation of the D3 domain anchored by the flexible loop between the D2 and D3 domains of IL-1R1. To characterize the rotational motion of the D3 relative to the D1-D2 domain, we defined a coordinate system as the following (Fig. 3B). The Cα atom of T217 at the C-terminus of the D2 domain is selected as the origin. The z-axis is defined as the vector connecting the origin and the center-of-mass of the D2 domain. The y-z plane is defined as the plane both the z-axis and the vector between the origin and the center-of-mass of the D1 domain lie on. Thus, the origin, the center-of-mass of the D1 and D2 domain are on the y-z plane. The x-axis is the axis perpendicular to the y-z plane. Relative orientation of the D3 domain to the D1-D2 domain is then represented by the vector connecting the origin and the center-of-mass of the D3 domain (Fig. 3B). Based on this coordinate system, the longitudinal motion of the D3 domain encompassed 110 degree (20 ─ 130) while the latitudinal motion covered a similar range of 100 degree (-30 ─ 70). The vectors from the origin to the center-of-mass of the D3 domain were mapped to a sphere shown in Fig. 3B to elucidate the large conformational variations of the D3 domain in the aMD simulation. Conformations of the IL-1Ra (1IRA)-, EBI-005(4GAF)-and AF10847(1G0Y)-bound IL-1R1 mapped on the sphere confirmed that they were included among those in the aMD simulations.

Small molecule binding site analysis of IL-1R1 using Sitemap
We first performed the Sitemap analysis to identify the small molecule binding sites in the crystal structures of three different ligand-bound IL-1R1 conformations. Two common binding sites were identified in the EBI-005-and IL-1Ra-bound IL-1R1 structures (cf. Fig. 4A1, 4A2, 4B1, 4B2). Both sites are located at the interface between D1 and D2 domains. A similar binding site at the same location was detected in the AF10847-bound IL-1R1 structure (Fig. 4C2). Two smaller binding sites close to the loop between D2 and D3 domains were detected in the EBI-005-bound IL-1R1 structure (Figs. 4A3 and 4A4). The top two ranked binding sites based on the Dscore values in the AF10847-bound IL-1R1 are much larger in sizes (volume: 738 and 439 respectively) than most other binding sites with the typical volume of less than 200 (Table 1). For inhibitor development consideration, small molecules bound to S1, S2 in EBI-005-and IL-1Ra-bound IL-1R1 conformations may directly inhibit IL-1R1 binding with the endogenous proteins that required the D1-D2 domains for high binding affinity such as IL-1Ra [28]. In contrast, small molecules bound to S3 and S4 of IL-1Ra can potentially impact on the relative orientation of D3 domain to the D1-D2 domains of IL-1R1. Because AF10847 antagonizes IL-1R1 by stabilizing IL-1R1 in an alternative conformation, small molecules bound to S1 of the AF10847-bound IL-1R1 may act as gluing molecules to trap Il-1R1 at this alternative conformation. To assess the potential of developing small molecule inhibitors targeting these sites, we used a druggability index discussed by Halgren's Sitemap paper [22] which suggested undruggable binding sites have Dscore values < 0.83. Among the 15 binding sites detected as shown in Fig. 4, S1 Fig. and S2 Fig. are consensus druggable sites in IL-1R1 for small molecule inhibitors development according to the three crystal structures analyzed here (see Table 1). The S1 and S2 sites both have larger hydrophilic than hydrophobic surface area.    Subsequent Sitemap analysis was performed on the 25 representative conformations to identify small molecule binding sites. As shown in Table 2, 32 of the 121 detected binding sites (27%) gave the Dscore ! 0.83 which is the same as those detected using the three crystal structures (27%). The distribution of the Dscore values for all binding sites was found to peak at 0.7 in Fig. 6. When analyzing the volume of each binding site, conformation 490 has a very large binding site with a volume of 1330 Å 3 which is an outlier when compared with the sizes of the binding sites in all other conformations. The average volume of the druggable binding sites (Dscore ! 0.83) detected from the 25 representative conformations is 285.54 Å 3 which is greater than those in the EBI-005-and IL-1Ra-bound IL-1R1 conformations (137.54 Å 3 ). The differences of the binding site volumes are equivalent to a ligand of 34 versus 11 atoms. Although the peak of the volume distribution is at 150 Å 3 , an extended tail population between 300-500 Å 3 can be found (Fig. 6). The volume sizes at the binding sites may be helpful for selecting appropriate sizes of compounds when building a library of compounds for screening. The Dscore value is also linearly proportional to the volume of the binding site although a steeper increase was seen after Dscore = 0.83 in Fig. 6. For the druggable binding sites with Dscore ! 0.83, the binding site volumes are around 200 Å 3 and higher. All probe points in the druggable sites of the 25 conformations from Sitemap were shown in grey color in three figures in Fig. 7. Most points were found clustered at three common regions (labeled P1, P2 and P3 in Fig. 7) which may be useful for developing inhibitors to prevent its binding with known ligands or small molecule modulators to affect the IL-1R1 ectodomain conformation (cf. Fig. 1). For example, P1 pocket can be used to develop inhibitors because it is located at the interface between the D1 and D2 domains which interact with IL-1β, IL-1Ra and AF10847. Small molecules binding to the P3 pocket located at the interface between the D1 and D3 domains may stabilize the inactive conformation of IL-1R1 similarly to that achieved by AF10847. The P2 pocket is located at the loop region between the D2 and D3 domains and was identified as binding sites in at least two different relative orientations of D3 to D2 domains as those shown in Figs. 4A3, 4A4 and 7A. Among the 25 conformations, those with the D3 domain orienting opposite to the D1-D2 domain (see Fig. 7A) have not been observed by experiments.

Small molecule binding hotspot analysis of IL-1R1 based on the cosolvent mapping method
An emerging inhibitor development is to target potential binding sites not directly involved with protein-ligand binding (or allosteric binding sites) [48][49][50]. We investigated these types of binding sites in 32 druggable sites of 25 representative IL-1R1 conformations. Druggable sites were detected at the P2 and P3 sites in four conformations, i.e. 289, 496, 524 and 903 (Fig. 8A  and 8B). While the P2 site was assessed druggable in conformations 289, 496 and 903, the P3 site were determined druggable in conformations 496 and 903 (see Fig. 8A, 8B and Table 2). To account for the binding site flexibility and its influence on the binding site druggability, we employed the cosolvent mapping method [23][24][25] for the binding hotspots assessment. Here, we used the phenol as the cosolvent molecule to probe these four IL-1R1 conformations because the phenyl group is frequently found in fragment screening library used to evaluate targets involved in protein-protein interaction [36].
Based on the 4 ns cosolvent MD simulations, we found that conformations 496 and 903 were stable during the simulations when their conformations were mapped to the PC1 and PC2 subspace (Fig. 8C). In contrast, conformations 289 and 524 underwent continuous conformational changes and deviated from their initial conformations. While conformation 289 quickly became trapped at a nearby region, conformation 524 wandered to a larger conformational subspace projected by PC1 and PC2 (Fig. 8C). The cosolvent mapping analysis further characterized the location of the hydrophobic hotspots preferentially bound by the phenyl molecules at the P1, P2 and P3 sites in these four conformations ( Fig. 8D and S3 Fig.). Binding modes of phenol molecules at the P1 site in conformation 524 indicated greater flexibility of the P1 pocket as they interacted with IL-1R1 at multiple locations in the same region (S3 Fig.). For conformations 496 and 903, less mobility of the phenol molecules at the P2 and P3 sites were found (S3 Fig.). For conformation 289, the phenol molecules are found to be confined at smaller pockets in the P1 and P2 sites ( Fig. 8D and S3 Fig.).
In the cosolvent MD simulations, the specific site found to interact with the probe molecules with less mobility reflects ideal complementary and can be considered as a hotspot region. Assessment of the results from cosolvent MD simulations indicated that the P2 site of conformation 289 and 903 and the P3 site of conformations 496 and 903 are binding hotspots using phenol as probes. The druggable property of the P2 site in conformation 524 was transient as the metastable conformation 524 deviated substantially from its initial conformation during the cosolvent MD simulation. The stability of each conformation was further evaluated by the 8 ns conventional MD simulations in which similar trends of conformational changes were observed (cf. Fig. 8C and S4 Fig.). Because conformation 289 represented a novel conformation of IL-1R1 not reported previously, we then focus on identifying candidate molecules that docked favorably into the P2 site of conformation 289 and assess their impacts on the IL-1R1 conformations.
Targeting the allosteric modulator site of conformation 289 using in silico fragment library screening followed by MD simulations of the top eight ranked ligands with conformation 289 Based on the Sitemap evaluation and the cosolvent mapping analysis, the P2 site of conformation 289 warrants further investigation to determine how small molecules bound to the P2 site will affect IL-1R1 adopting the 289 conformation. Because the comparable size of P2 (352 Å 3 ) to fragment compounds, we performed in silico screening to target the P2 site of conformation 289 using the Maybridge fragment library (a total of 2783 compounds). The Glide docking program was used to rank the predicted potencies of the compounds to the P2 site of conformation 289 based on their docked poses. The highest eight ranked compounds (denoted L154, L537, L951, L1192, L1206, L1882, L3097 and L3241) were selected and subject to additional 8 ns MD simulations to evaluate the effects of the ligand binding on IL-1R1 conformations.
In Fig. 9, we calculated the root-mean-square deviations (RMSDs) of the backbone atoms in the D1-D2, D3 domains and the heavy atoms of the compounds in reference to the initial ligand poses and conformation 289 in the 8 ns MD simulations. All IL-1R1 conformations in the MD simulations were aligned to the D-D2 domain of the initial conformation 289. Similar to the observation in the conformational sampling calculations, the D1-D2 domain in all eight fragment-bound IL-1R1 conformations exhibited less than 3 Å backbone deviations (Fig. 9A) except those with L154. The movement of the D3 domain is much greater even with the binding of these fragment ligands. The RMSD distributions of the D3 domain are less than 20 Å for L1192, L3097, L1206, 30 Å for L951, L1882, L537, and 40 Å for L154, L3241 (Fig. 9B). When monitoring the deviations of the ligands from their initial binding site locations, we found four ligands (L951, L1882, L1192 and L537) have major RMSD peak distributions at less than 5 Å whereas the other four (L1206, L3097, L154 and L3241) have broader distributions of RMSDs at greater than 5 Å (Fig. 9C). The RMSD values of L3097 along the trajectory (Fig. 9D) informed that L3097 escaped from the binding site. In contrast, L951 and L537 remained at the initial docked positions throughout 8 ns MD simulations (the peak distribution of RMSDs at around 2.5 Å). L1882 was close to its initial docked position only up to 4.8 ns whereas L1192 was up to 3 ns.
To investigate the conformational changes of these ligand-bound IL-1R1 conformations, we mapped them to the PC subspaces derived from the crystal structures. Two dimensional mapping to PC1/PC2 and PC1/PC3 were provided in Fig. 10. The mapping to PC1/PC3 subspace was provided because we found dynamical motion of IL-1R1 along the PC3 in the simulations. The results showed that IL-1R1 conformations obtained from four ligands (L951, L1882, L1192, L537) exhibiting smaller average RMSD values were projected to much localized regions in the PC1/PC2 and PC1/PC3 subspaces. In contrast, larger conformational spaces were visited by IL-1R1 bound with L1206, L154 and L3241. Although L3097 escaped from the docking site at around 0.9 ns, yet it caused IL-1R1 to trap at a local microstate. Based on the mapping of conformations to PC1 and PC2, the ligand-bound IL-1R1 conformations do not overlap with ligandfree IL-1R1 conformations but resembled those obtained from the cosolvent MD simulation.
Two dimensional mapping in Fig. 10 implicated certain overlap of conformations between L951-and L537-bound IL-1R1 and IL-1β-bound IL-1R1. When we examined the IL-1R1 conformations, we found that they do not resemble the conformations of IL-1R1 conformations revealed by the crystal structures. This agrees with the large RMSD values of the D3 domain in IL-1R1 shown in Fig. 9 and the limitation of the two dimensional map. To elucidate the modulations of IL-1R1 conformations by these eight ligands, mapping of conformations to the three dimensional subspaces by PC1, PC2 and PC3 was used and shown in Fig. 11. IL-1R1 conformations were either locally trapped when bound with L1192 or transited along PC2 degree of freedom for L951 and L537. For L1882, IL-1R1 conformations were locally trapped for the first 4.8 ns and moved to the subspace closer to the IL-1β-bound conformation after L1882 escaped away from the binding site (Fig. 9D). Among the four other ligands, the initial interaction between L3097 and IL-1R1 caused the IL-1R1 to be trapped at a localized subspace even after L3097 dissociated from the P2 site. Broader conformational changes of IL-1R1 were found for the other three ligand-bound IL-1R1 conformations in which the ligands dissociated from the P2 site relatively early during the simulations. Ligands bound to the P2 site in IL-1R1 can Potential IL-1R Allosteric Modulators and Auto-Inflammatory Diseases potentially modulate the IL-1R1 conformations and impact on IL-1R1 binding with other protein ligands. Thus, L951, L1882, L1192 and L537 are attractive candidates for future experimental validation and evaluation.
Figs. 10 and 11 indicated that L951, L1882, L1192 and L537 bound IL-1R1 conformations underwent less conformational exploration in the PC subspaces (Figs 10 and 11). We further analyzed the protein-ligand poses of these four ligands at 2 ns of MD simulation where the ligands did not deviate much from their initial poses (RMSDs < 5) and interaction between IL-1R1 and ligands were established. We found a consensus binding site motif at the P2 site of IL-1R1 that interacts with the ligands via a similar type of protein-ligand interaction (Fig. 12). All four ligands interact with the binding pocket formed by R272, R271, V117, E202, E171, R174 and K172 between the D2 and D3 domains of IL-1R1. Salt-bridge interactions between Arg, Lys and Glu at the interface between the D2 and D3 domains establish a mall binding pocket for ligand binding. At the binding site, hydrogen bonds are formed between the hydroxyl groups of L951, L1882, L1192, L537 and E171, R272, E202 respectively. The carbonyl group of L951 and the thiophene group of L1882 form additional hydrogen bonds with R271 and R272. Furthermore, the aromatic rings of L951, L1192 and L537 interact with R271 and R272 via cation-aromatic interactions [51]. In this binding site, V117 is the key hydrophobic amino acid to interact with all ligands.

Conclusion
Existing crystal structures of IL-1R1 with three different ligands including IL-1β, IL-1Ra and AF10847 indicated the large conformational flexibility of the loop between D1-D2 and D3 domains. We hypothesized the existence of transient small molecule binding sites among multiple conformational states accessed by IL-1R1 which may be suitable for inhibitor development. To identify these binding sites, we performed efficient conformational sampling of IL-1R1 using the accelerated MD algorithm. In a total of 118 ns aMD simulations started with a super cytokine (EBI-005)-bound IL-1R1 structure, we observed that the D3 domain of IL-1R1 exhibited a wide angular motion relative to the D1-D2 domain facilitated by the flexible loop between the D2 and D3 domains. These include the IL-1R1 conformations with close contacts between the D1-D2 and D3 domains mimicking closely to the AF10847-bound IL-1R1 conformation. In addition, we identified unreported IL-1R1 conformations with the D3 domain oriented to a position different from either IL-1Ra-or AF10847-bound IL-1R1 conformations.
Using the Sitemap analyses, we identified small molecule binding sites in the IL-1R1 conformations. The druggability index based on the Dscore value of these binding sites was then used to select "druggable" sites in these conformations for further evaluation. Four binding sites in three IL-1R1 crystal structures and 32 binding sites in 25 representative IL-1R1 conformations obtained from simulations gave the Dscore values ! 0.83 suggesting that they are potentially druggable. These druggable binding sites were found at three common locations in IL-1R1 including 1) the region between the D1 and D2 domain (the P1 Site), 2) the location between the D2 and D3 domains (the P2 site) and 3) the interface region between the D1 and D3 domain (the P3 site) when IL-1R1 adopts an inactive conformation.
Because the P1 and P3 sites can also be identified from existing crystal structures, we focused on investigating compounds potentially bound to the P2 site, an allosteric modulator site. Among the 25 representative conformations classified by the hierarchy cluster analysis, the P2 site in four conformations (289, 496, 524 and 903) gave Dscore values ! 0.83 indicating they are suitable for small molecule ligands binding. To account for the flexibility of the binding sites not included in the Sitemap analysis, we performed the cosolvent mapping analysis. Only the P2 site of conformation 289 and 903 were confirmed to be small molecule binding hotspots.
Because conformation 289 is a novel conformation not reported previously, an in silico screening by targeting the P2 site of conformation 289 was conducted using the Maybridge fragment library. Eight highest ranked ligands were selected and subjected to additional 8 ns conformation 289/ligands MD simulations. Four ligands (L951, L1882, L1192 and L537) bound to the P2 site of conformation 289 for an extended period of time in the simulations and restricted the IL-1R1 conformations to a relatively localized region in the principal component subspaces. We further identified key residues in IL-1R1, including R271, R272, E202, E171 and V117, forming the binding pocket at the P2 site to interact with the ligands. The results suggested small molecule ligands can potentially bind to the P2 site of IL-1R1 and modulate the IL-1R1 conformations. Because conformation 289 differs from the IL-1R1 conformations observed in the crystal structures, potent chemical probes developed to target the P2 site of conformation 289 can potentially restrict IL-1R1 to an inactive conformation which abrogates the downstream IL-1 signaling. Additional in vitro experimental evaluation can validate our approach of discovering small molecule compounds targeting the allosteric modulator site in the ectodomain of IL-1R1. This strategy is attractive for the development of small molecule therapeutics targeting this difficult target that involved in protein-protein interaction and can potentially be applicable to other interleukin receptor proteins.