Glu106 targeted inhibitors of ORAI1 as potential Ca2+ release-activated Ca2+ (CRAC) channel blockers - molecular modeling and docking studies.

Abstract Calcium release-activated calcium modulator 1(ORAI1) is an integral component of the calcium release-activated calcium channel (CRAC) channel complex and plays a central role in regulating Ca2 + concentrations in T-lymphocytes. It is critical for many physiological processes, including cell-proliferation, cytokine production and activation of the immune system. Loss of ORAI1 function is linked with rheumatoid arthritis (RA) and hence pharmacological blockers of ORAI1 could be potential therapeutic agents for the treatment of RA. In this study, we have used a high-throughput screening approach to inhibit the binding of Ca2+ toward ORAI1 and the interactions are verified through induced fit docking. The results hint that these compounds act by possibly binding with, and thereby blocking Ca2+-binding with ORAI1 (E106). The molecular dynamics (MD) simulations shows strong support toward the hit compounds by showing the ligand potency throughout the simulation timescale of 30 ns. We have thus identified a novel class of highly stable, potential lead compounds that directly bind with the selectivity filter region E106 and block Ca2+ binding on ORAI1. This resulting alteration in the pore geometry of ORAI1 due to the strong blocking mechanism of lead compounds will greatly diminish its function and the downstream activities that result from the same including decreased production of cytokines in autoimmune disorders. This study may lay the foundation for finding novel lead compounds for clinical trials that could positively modulate the course of autoimmune disorders with ORAI1 as its specific target.


Introduction
Rheumatoid arthritis (RA) is a chronic inflammatory autoimmune disease that mainly affects the joints of the human skeletal system (1)(2)(3). RA is cited to be the second most common autoimmune disorder (AID) with a prevalence rate of 1-1.5% in the world population (4). However, there is still a lacuna in understanding its cause and in effecting a therapy for cure of the same (5). Ongoing research has focused on several functional ion channel networks associated with RA and recent literature show that ion channels could prove to be potential drug targets for treating RA (6)(7)(8)(9)(10)(11). Of particular mention are the Calcium (Ca 2+ ) release-activated calcium (CRAC) channels that have proven to have a major role in the onset and progression of RA, which projects it as a potential target protein to develop suitable modulators against (12,13). CRAC channel is a major Ca 2+ selective store-operated channel activated by the binding of stromal interaction molecule 1 (STIM1) with calcium release-activated calcium modulator 1 (ORAI1) and allows the influx of Ca 2+ (14). This particular regulated control of Ca 2+ influx activates the nuclear factor of activated T-cells (NFAT) pathway and plays a vital role in normal T-cell activation, proliferation, gene expression and cytokine production (15)(16)(17)(18)(19). The overall mechanism of T-cell activation and subsequent Ca 2+ signaling pathway involved in CRAC channel activation are clearly elucidated in Figure 1.
As described in Figure 1, CRAC channels orchestrate the function of T-lymphocytes and mediate the long lasting influx of Calcium in abnormal conditions. Loss of function mutants of either STIM1 or ORAI1 has been linked to several autoimmune diseases including RA. Various ORAI1 mutations and genetic polymorphism of ORAI1 gene exhibit defective CRAC channel function in patients with severe combined immune deficiency and it is also reported that overexpression of ORAI1 is responsible for abnormal cytokine production of T-cells in RA patients. Additionally, Picard et al. (2010) has reported that mutation in STIM1 is associated with severe autoimmunity, whereas ORAI1 mutation causes only a major clinical syndrome of immunodeficiency (20)(21)(22)(23). Modulating ORAI1 function by hampering its binding with Ca 2+ could aid in regulating and controlling the progression of the disease. Several CRAC channel inhibitor compounds, including Synta-66, GSK-5503A, GSK-7975A, YM58483, SKF96365 and 2-APB have been found to be successful in controlling cytokine release in animal models. Unfortunately, none of the compounds with the exception of CM2489 had progressed to clinical trials due to lack of specificity and their cross reactivity with other ion channel signaling process (24)(25)(26)(27). Therefore, it becomes the need of the hour to identify potent, novel drug candidates that block Ca 2+ binding with ORAI1 with a high specificity. However, even now the 3D structure of ORAI1 in Homo sapiens has not been elucidated in detail. Hence, this study focused on the prediction of ORAI1 3D structure and screening of compounds that inhibit the binding of Ca 2+ with ORAI1 by specifically targeting the active site of ORAI1 through molecular modeling and molecular dynamics (MD) simulation approaches.

Protein modeling and structure validation
The 3D structure of the H. sapiens ORAI1 is unavailable in protein PDB. Hence, the protein sequence of ORAI1 is retrieved from UniProt (ID: Q96D31) and subjected to Figure 1. Schematic representation of T-cell activation and Ca 2+ release-activated Ca 2+ (CRAC) channel gating mechanism in rheumatoid arthritis. [1] The antigen/major histocompatibility complex (MHC) proteins present to the T-cell receptor/CD3 complex. [2] The stimulation of T-cell receptor activates phospholipase (PLC W ) and results in the production of inositol-1,4,5-trisphosphate (InsP 3 ). [3,4] The generated IP 3 binds to its receptor (InsP 3 R) positioned in the membrane of the endoplasmic reticulum (ER) and this complex acts as a ligand gated ion channel to promote the release of Ca 2+ ions from ER store. [5] In resting lymphocytes, stromal interaction molecule1 (STIM1) is homogenously distributed in the endoplasmic reticulum with EF hand loaded with calcium ions. Following Ca 2+ depletion concentration in ER store sensed by STIM1, Ca 2+ unbinds from the EF hand. [6,7] STIM1 senses the depletion of Ca 2+ and forms oligomerization to form discrete puncta state in endoplasmic reticulumplasma membrane junction. [8] STIM1 captures diffusing ORAI1 channels. The N-and C-terminal of the ORAI1 interact with CRAC activating domain (CAD) on STIM1 leads to CRAC channel opening. [9] The Ca 2+ influx mediated by the pore forming subunit ORAI1 provides the sustained Ca 2+ flux and binds to calmodulin (CaM), which in turn activates the calcineurin and leads to dephosphorylation and translocation of the nuclear factor of activated T-cells (NFAT) into the nucleus. [10] NFAT and other transcription factors trigger the transcription of genes and involve in T-cell proliferation, differentiation and cytokine production in normal condition. [11] Under abnormal circumstances, up-regulation of Ca 2+ entry through CRAC channel in CD4 + T-cells leads to increase in cytokine production, which in turn produces inflammation in the synovial membrane in case of rheumatoid arthritis.
BLAST P search to find the suitable template structure for building the 3D structure of the protein. The percentage identity and query coverage against target sequence is less than the acceptable range (525%) for the homology modeling. Due to lack of suitable templates, the sequence is modeled via threading method by using model building server I-TASSER (http://zhanglab.ccmb.med.umich.edu/I-TASSER). I-TASSER is ranked based on confidence score, TM score, root mean square deviation (RMSD) and the standard deviation (28). The obtained structures were optimized and minimized through protein preparation wizard. The stable and average structure obtained from molecular dynamics simulations are further used for molecular modeling calculations and the final stable structure is provided in Figure 2. Structural Analysis and Verification Server (SAVS) is used to check stereo chemical qualities of a refined ORAI1 protein structure by analyzing the residue-by-residue geometry and overall structure geometry in the PROCHECK program (29). The Ramachandran plot displays the number of the residues in the most favorable, allowed, generously allowed and in disallowed regions and PROSESS server is used to check the global structural assessment by evaluating the bond quality (30). The overall quality of the model thus generated was validated using QMEAN, ANNOLEA and GROMOS. Further, the secondary structure and topology of the modeled ORAI1 were viewed using PDBsum.

Molecular dynamics simulations
The MD simulation studies are carried out for the modeled ORAI1 protein with the intention of understanding the stability and intra-molecular conformational changes in the protein that occurred during simulation. Simulations were carried out through academic version of Desmond molecular dynamics package with OPLS2005 force field for 30 ns (31,32). The modeled protein was solvated in orthorhombic box with 1-palmytoil-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC) lipids water molecules and neutralized by adding salt counter-ions. The distance between the box wall and complexes was set to410 Å to avoid direct interaction with its own periodic image. Prior to the simulation of the activation process, the system was energy minimized using the steepest descent method for a maximum force of more than 100 KJ/ mol À1 /nm À 1 , then gradually heated to 300 K over 100 ps under constant NVT, and equilibrated under constant NPgT (pressure 1 bar) with a time step of 2 fs and relaxation time was applied between 0.1 and 0.4 fs (33,34). The van der Waals (vdw) force was treated by using a cut-off of 12 Å and the final MD simulation was carried out for 30 ns. The Ca RMSD and root mean square fluctuation (RMSF) of ORAI1 protein were plotted for the entire simulation trajectory with reference to their respective first frames. Additionally, MD simulations were carried out for the top-ranked protein-ligand complex to understand the stability of the interacting residues during simulation using the above mentioned setup. Steepestdescent energy minimization method was applied to proteinligand complexes to remove close contacts, terminating when maximum force is found to be 5100 KJ/mol À1 /nm À1 . All systems were subject up to 100 steepest descent energy minimization steps before thermalization for a maximum force of more than 100 KJ/mol À1 /nm À1 . After thermalization MD simulations were run at the constant NPT (300 K) and NPT (pressure 1 bar) with a time step of 2 fs and the relaxation time was applied between 0.1 and 0.4 fs and NVT was carried out for 1 ns (35,36). This equilibrated system was used for molecular simulation with a time step of 30 ns. Ca-RMSD, RMSF and hydrogen bond interactions of the complex protein were plotted for the entire 30 ns simulations with reference to their respective first frames. Further, the residue interactions with the ligand in each trajectory frame were plotted to know the specific amino acid residue that is involved in the stability of the ligand inside the protein. In addition, the stability percentage of ligand protein contacts are plotted to identify the potency of the ligand within the protein during simulation.

Binding site prediction
There is no accurate information on the binding site of the ORAI1 protein and hence a prediction of the same would serve to enhance the drug discovery process. The binding site of the ORAI1 is manually predicted through sitemap module using Sitemap version 3.0 (Schrodinger ,2015, LLC, New York, NY). The typical active and stable conformation of the ORAI1 protein structure chosen from the MD trajectories was subjected for the active site prediction. Druggability regions were identified based on the various physical descriptors like the size, the degree of enclosure, the degree of exposure, tightness, hydrophobic, hydrophilic, hydrogen bonding possibilities and linking site points that are more likely to contribute in protein-ligand binding (37,38). The top score active site region is found near to the Ca 2+ binding region and clear focus that inhibitors bound to this particular location will drift the Ca 2+ binding property. The charge distribution of the predicted active site region is viewed using electrostatic potential surface using PyMol (Schrodinger, 2015, LLC, New York, NY).

Structure-based virtual screening
Structure-based virtual screening is used to screen the potential lead compounds on the basis of active site features of the target protein. The top score active site region was chosen and the grid was successfully generated around the active site region by using receptor grid generation. To soften the potential for the non-polar part of the receptor, we scaled van der Waals radii of receptor atoms by 2.0 Å with a partial charge cut-off of 0.25. The small molecule library was obtained from the Maybridge database and around 14 million compounds were subjected toward the screening protocol. Initially, the compounds were screened with Lipinski filter, absorption, distribution, metabolism and excretion (ADME) profiling and heavy metal atoms are excluded (39)(40)(41)(42)(43). The screen accepts the best compounds with ADME, Lipinski filter and without heavy metal atoms are subjected to dock with funnel-based docking protocol (high-throughput virtual screening (HTVS), standard precision (SP) and extra precision (XP) docking). In each phase of screening, based on scoring parameters, the best compounds were chosen for the next phase. The compounds with least scoring value were eliminated and the top scoring hit compounds automatically forward to the next phase. Final Glide XP mode determines the compounds with lowest energy conformer with better scoring values in the predicted active site and 3D visualizations will focus the interactions between the protein and ligand (44).

Induced fit docking
The top score compounds screened using the Glide scoring function (G-score) were further refined through Induced Fit Docking (IFD) (45), here both the protein and ligand were in flexible condition. In this protocol, first the ligands were docked into the rigid protein using softened potential docking in the Glide program with the van der Waals radii scaling of 0.8 and energy minimization was carried out using OPLS-2005 force field with an implicit solvation model (46). As both protein and ligand in flexible condition with multiple receptor conformations will provide the maximum possibility of accurate interactions. Here, each docked conformers in previous step were subjected to side-chain trimming and backbone refinements through prime. Different possible conformations are generated through prime and multiple conformations are docked through Glide. The refined complexes were ranked by prime energy, and conformers within À20 kcal/mol of the minimum energy structure were passed for a final round of Glide docking and scoring (47). The side-chain orientations have been performed automatically with inclusion of prime in IFD. An IFD score that accounts for both the substrate bound protein-ligand interaction energy and the total energy of the system were calculated. Further optimization of these compounds may lead to synthesizable compounds for the development of a novel drug.

Free energy calculation
Binding energy calculations are more accurate than docking energy calculations. Hence, in this study, we used a molecular mechanics/generalized born surface area (MM/GBSA) approach to estimate relative binding affinity for the list of selected compounds against ORAI1. Prime uses a surfacegeneralized Born model employing a Gaussian surface instead of a van der Waals surface for better representation of the solvent accessible surface area (48,49). The binding free energy, DG bind, was calculated using the following equation.
where E complex , E protein and E ligand are the minimized energies of the protein inhibitor complex, protein and inhibitor, respectively.
where G solv (complex), G solv ( protein) and G solv (ligand) are the solvation-free energies of the complex, protein and inhibitor, respectively.
where G SA (complex), G SA (protein) and G SA (ligand) are the surface area energies for the complex, protein and ligand, respectively.

Protein modeling and validations
The importance of molecular modeling of the ORAI1 3D structure rests with the fact that there is an absolute lacuna in the knowledge of the exact structure of the protein given the absence of any crystal structure data. As for the focus of template structure, 525% of similarity tends to follow threading methods by using I-Tasser server. (2010) thoroughly studied the important role of N-and Cterminal of ORAI1 protein and notably found that nearly 107 residues in the CRAC activation domain region (CAD) of STIM1 bound with both terminals of ORAI1 and significantly participate in the gating mechanism of CRAC channels. However, just deleting the N-terminal region of 1-64 amino acid of ORAI1 or mutation in those regions does not affect the CRAC gating function (50)(51)(52)(53). Therefore, in this study 1-64 amino acids were trimmed and the remaining residues were subjected for modeling studies. I-Tasser provides the model by choosing the multiple templates as per the regions allocates in similarity. Multiple template-based threading are choosy and trims the required identity and provides five successive models; among that, model 2 with a high confidence score (À2.52) and with lesser number of bad contacts was chosen for further validation process. As per literature studies, the pore region of the CRAC channel consists of both high and low binding affinity sites. At first, the incoming Ca 2+ ions bind at the high-affinity calcium binding site (E106) and the following Ca 2+ ions bind at the low-affinity site. Due to the electrostatic repulsion between the ions, the calcium ion bound at the high-affinity site forces the one bound at the low-affinity site down the pore, thereby leading to a flow of current through the CRAC channel (54,55). Our aim is to find out the potential blocker which blocks the influx of calcium during the open state of CRAC channel. Therefore, in our model to acquire active state of the protein, Ca 2+ was charged and embedded in the E106 regiona high-affinity calcium binding site for ion selectivity and the opening of the channel. Further, reliability and overall quality of the predicted structure is evaluated by three validation servers. The stereo chemical property of the modeled ORAI1 is evaluated by Ramachandran plot (PROCHECK, Los angeles, CA) and it shows 88.6% amino acids in favor region, 7.6% in allowed, 1.4% in general and 2.4% of amino acids in the disallowed region (Supplementary Figure S1). Global structure assessment is evaluated by PROSESS (Alberta, Canada) server and it shows that the covalent bond quality, noncovalent bond quality and torsion angle quality are within the acceptable range (Supplementary Figure S2). After that, SWISS-MODEL (Basel, Switzerland) server is used to validate the overall quality using QMEAN (Basel, Switzerland), ANOLEA (Basel, Switzerland), GROMOS (Basel, Switzerland) (Supplementary Figure S3) and the overall results obtained from the above analyses reveal that the 3D model of ORAI1 with Ca 2+ interaction, as shown in Figure 2, is reliable and satisfactory for further analysis. Further modeling studies are carried with the stable and least energy minimized conformation obtained from MD simulations.

Architecture of ORAI1: protein structure and topology
The architecture ORAI1 protein is composed of four transmembrane a-helices (TM1-TM4) with two extracellular and one intracellular loop. These transmembrane helices of each ORAI1 subunits are arranged in three concentric rings and form the CRAC channel. TM1 subunits form the inner ring and this is liable for the ion pore region. The TM2 and TM3 a-helices form the middle ring, which provides the structural integrity to pore region. Followed by that TM4 helices form the outer peripheral ring, which is far away from the pore region (10). The secondary structure and the functional environment of the predicted ORAI1 are elucidated using PDBsum (Cambridgeshire, United Kingdom). The topology view of the modeled ORAI1 protein has 8 helix, 36 beta turns and 7 gamma turns and it is composed of 62% of helix, 25% of turns and 38% of the coil region (Supplementary Figure  S4). Starting with the N-terminus Ser 65-His 69 forms the beta turns in the cytosolic region, which is followed by Ser70-Ala103 that forms the a-helix and this region is significantly involved in gating mechanism of the CRAC channel by interacting with STIM1. Met104-Gly118 forms the beta turns, the residues that lie in these regions actively participate in Ca 2+ selectivity and also in the pore forming unit. Followed by those, Leu119-Asn147 forms the helical structure and Ile148-Asn156 contribute in the beta turns of the ORAI1 in the cytoplasmic region. Continued with that, Leu157-Ser159 forms the short a-helix and Val160-Ser163 forms the beta turns. Further, the residues from Pro164 to Leu200 form the transmembrane helical region and Pro201-Thr230 form the beta turns in the extracellular region. Pro231-Leu276 forms the helical region, followed by which Ala277-Glu278 forms the beta turns. Again in the C-terminal, amino acids from 272 to 279 form the a-helix which interact with STIM1 and are involved in the CRAC channel opening mechanism. Further, Ala277 and Glu278 form the short beta turns and Phe279-Asp291 forms the a-helix. In the end, residues His292-Ala301form the beta turns in the cytosolic region.

MD simulations
MD simulation for the apo structure is performed to understand the structure stability and conformational flexibility of the modeled ORAI1 protein. As CRAC channel protein comes under the family of a-helix bundle of the integral membrane protein, most of the TM residue of ORAI1 interacts with the lipid bilayer environment. In resting lymphocytes, STIM1 is homogeneously distributed in the endoplasmic reticulum with EF hand loaded with calcium ions and ORAI1 is scattered across the plasma membrane (PM). Following Ca 2+ depletion concentration in endoplasmic reticulum (ER) store sensed by STIM1 forms oligomers and recruits diffusing ORAI1 in the PM to form a discrete puncta state in endoplasmic reticulum-plasma membrane junction (12). Given the important role of protein-lipid interactions in structural stability and in the gating mechanism of the channel, here we want to mimic the same environment to explore the specific interaction of ORAI1 protein with phosphate groups of the plasma membrane which facilitate the ORAI1 oligomerization for the CRAC gating mechanism. Therefore, the ORAI1 protein was embedded in the POPC (56)(57)(58)(59). The dynamic behavior of ORAI1 protein is calculated by the C-a atom RMSD of the final conformation with that of the initial structure. The average RMSD value shows that, there is more fluctuation in the initial 5 ns, which rise up to 4.6 Å and for the next 5 ns it further rise up to 5.5 Å . But in the time scale of 10-15 ns and 15-20 ns, it reaches up to 6 Å and remains stable at 6 Å after it reaches an equilibrium state in 20-30 ns (Figure 3). Average RMSD plot and timeline representation of secondary structure elements of modeled ORAI1 throughout the simulation trajectory are showed in Supplementary Figures S5 and S6. Further, the fluctuations in each residue are shown in RMSF plot with reference to protein secondary structure (Supplementary Figure S7). It shows the fluctuations were high in the loop region compared with other secondary structure elements. Apart from the loop region, the architecture of the protein remained moderately similar even after simulation. From the trajectory analysis, it is hypothesized that the specific binding interaction between the TM segments of charged amino acids such as negatively charged residues like glutamic acid and aspartic acid in the TM1 near the selectivity region and the positively charged Lys and His in the hydrophobic region and further Ser, Tyr and Thr AA of the side-chains charged amino acids collectively entail in interacting with phosphate group of the lipid layer, which facilitates in the ORAI1 oligomerization and this mechanism assists in opening of the CRAC channel. Together, these results confirm that the protein is stable under simulation and based on the RMSD value the average stable and least energy minimized conformation of the protein obtained at 6 Å was chosen for further screening studies.

Binding site prediction
Given the critical role of the CRAC channel mechanism in many autoimmune disorders, ORAI1 becomes the promising target for CRAC channel blockers. Several potent CRAC channel blockers were synthesized and their inhibitory effects were well studied through electrophysiology and calciumbased fluorescence measurements. Yet, the binding site information and the detailed mechanism of existing CRAC inhibitors that involve in blocking I CRAC have not yet been reported. As the identification of active site region of the protein is the key step in drug designing, stable conformation of ORAI1 protein resulted from MD studies is subject to Sitemap V3.0 to predict the active binding site location. Basically, the ORAI1 pore region distinguishes into four divisions namely extracellular selectivity filter region, hydrophobic region, basic section and cytosolic region and more details are already discussed in architecture of the ORAI1 protein. In detail, the pore is $55 Å long and the amino acids participating along the channel pore are E 106, V 102, F 99, L 95, R 91, K87, R 83, Y 80 and W76. Out of which six glutamate (E) 106 residues from each ORAI1 subunit form the selectivity filter region called Glu ring with high-affinity Ca 2+ binding site and it majorly involves in opening of the channel (15). Based on the interaction specificity on these regions, Site Map come up with top five ligand binding sites in 3D structural format and the best score binding site is interestingly, as expected, found in the N-terminal side of the TM1 a-helix near to the more negatively charged pore region which occupies 23.230 surface area for hydrogen bond acceptor, 209.360 surface area for hydrogen bond donor, 244.924 area for hydrophilic and 35.420 surface area for hydrophobic region. This indicates that the binding site region is more favorable for the hydrophilic compounds. The predicted druggability binding sites within 4 Å were E106, Tyr 115, Leu 120, Val 107, Val 105, Val 102 and Met 104 and the druggability pocket is clearly illustrated in Figure 4.

Electrostatic potential analysis
The binding of Ca 2+ ion in the extracellular region is driven to the intracellular cytoplasmic region, this is due to the electrostatic repulsion between the ions in the pore region. To view in detail, the electrostatic potential (ESP) surface of the protein were created using PyMol. It showed that, more negatively charged (red) residues were observed in the predicted active site region and in the entrance of the pore region that is in the extracellular hydrophilic region. Followed by that, the majority of the ORAI1 proteins is occupied with neutral surface (gray) present inside the hydrophobic region of the cell membrane. The hydrophilic region of the intracellular space is composed majorly with positively charged amino acids (blue). The electrostatic gradient is contoured from +59.748 (red) to À59.748 (blue) kT/e. The structure along with its surface potential is represented in Figure 5 as per color scheme of ESP.

Virtual screening and molecular docking studies
A combined approach of virtual screening and IFD was used to identify the potential inhibitor drug against ORAI1 activity. Around 14 million compounds from the Maybridge database were subjected to Virtual Screening Workflow (VSW). Initially, Qikprop (Schrodinger, 2015, LLC, New York, NY) program was used to evaluate the ADME properties and this is the most important step in the drug development process to eliminate the false-positive molecules with poor   pharmacokinetics property. Further, the ADME filter subjects to the Lipinski filter, and exclusion of heavy metal compounds reduce the compounds, which may show high toxicity and screened compounds are forwarded to docking analysis. Initial docking was performed through HTVS and here the compounds having less binding efficiency toward the ORAI1 active site are eliminated and around 83 compounds were sorted out and these resulting compounds further underwent SP method. Based on the Glide score and Glide energy values obtained from the SP method, the best compounds were forwarded to extra precision molecular docking process, up to which the receptor is fixed and the ligand is flexible. Finally, the resultant hit compounds were again re-docked using IFD approach where both the ligand and the receptor are flexible. IFD resulted with top four compounds namely SB01990, SPB06836, KM06293 and RHO1882 with the IFD score value of À457.109, À454.993, À454.463 and À453.771 kcal/ mol, respectively. Two dimensional conformations of the hit compounds are illustrated in Supplementary Figure S8 and from the structure, it is observed that these compounds comprise six-member heterocyclic pyridine groups. It is well defined that pyridine-containing drugs are considered to be the one of the best selling drugs in pharmaceutical industries due to their wide range of biological activity and for the simpler metabolism. Even though these metabolites are toxic to humans, these compounds are easily and rapidly excreted from the body. Therefore, understanding the pharmacokinetics properties and structural optimization of these drugs will be more effective for the oral administration.
Protein-ligand complex through hydrogen bond interaction plays a crucial role in drug design by determining its specificity of ligand binding. So, we accessed the length and orientation of the hydrogen bond interaction between ligand and receptor and the distance in the acceptance range of 1.50-2.74 Å indicates the presence of favorable interactions. Compound SB01990 shows the good IFD score of À457.109 kcal/mol with three hydrogen bond interactions with Val105, Val107 and Met104. SPB06836 shows the IFD score of À454.993 kcal/mol with two hydrogen bond interactions with Val105 and Val107. KM06293 and RH01882 show the IFD scores of À454.463 and À453.771 kcal/mol, respectively, with direct hydrogen bond interactions with Glu106 which indicates our ligands strongly bind in the Ca 2+ ion interacting sites. Molecular interactions of these top ranked compounds with ORAI1 are illustrated in Figure 6 (a-d). On the whole, E106, Gln 108, Lys 198, Ser 124, Tyr 115, Leu 120, Val 107, Val 105, Val 102 and Met 104 are the important amino acids involved in the atomic interaction between the protein-ligand complexes. Table 1 explains the detailed results of ligand-receptor binding affinity through IFD score, G-score, number of hydrogen bonds and it also provides the information about the interacting atoms that are involved in accepting and donating the atom, along with the length of each bond. Taken together, these four hit compounds showed better docking scores with satisfactory energy levels; further, these compounds are subjected to binding energy calculations and molecular dynamics simulation to study the stability of the ligand interaction in the active site region of the ORAI1.

MM/GBSA calculation of ligand-receptor complex
The efficiency of the ORAI1 protein complex with screened compounds were again confirmed by estimating the binding energy using MM/GBSA method and the results are shown in Table 2. The binding energy with more negative values indicates good binding affinity. From Table 2, it is observed that the screened compounds show satisfactorily strong binding energies (DG bind ) with À31.679, À23.047, À40.5 and À26.806 kcal/mol, respectively for SB01990, SPB06836, KM06293 and RH01882. While comparing the binding energy components it is observed that major favorable contributors to ligand are Coulomb energy (DG Coulomb ), van der Waals (DG vdw ) energy and nonpolar solvation terms (D GsolvSA ) energy, whereas, polar solvation energy (DG solvGB ) and covalent energy (DG covalent ) are not in favor of binding affinity with more positive values. From the results it is clearly evident that DG Coulomb energy majorly contributes for ligand binding and overall outcome indicates that ligands are strongly adopted inside the active site of the ORAI1 protein.

ADME
In this analysis, preceding the virtual screening, all compounds were calculated for the prediction of octanol/water partition co-efficient (QPlogPo/w), skin permeability (Qplogkp), percentage of human oral absorption, prediction of IC50 value for the blockage of HERG K + channels (Q PlogHERG) and prediction of binding to human serum albumin (Q PlogKhsa).The ADME properties calculated for the final four lead compounds which passed the acceptable range are shown in Table 3 with QPlogPo/w, 5.272-11.4041; Qplog kp, À2.956 to À4.456; percentage of human oral absorption, 62.205-76.939; Q PlogKhsa, À0.403 to À0.776; Q PlogHERG, À1.316 to À2.852 range and it satisfies the Lipinski rule of five. Therefore, these lead compounds are theoretically considered to have the drug-like properties. Recommended ranges for the predicted physicochemical descriptors are clearly shown in Supplementary Table S1.

MD simulations of protein-ligand complexes
MD simulation was carried out for the final four ORAI1ligand complexes to understand its structural stability and dynamic properties of the ORAI1 and newly screened complex for 30 ns. The MD simulations for bound structure of ORAI1 are analyzed by using the Ca-RMSD, RMSF, hydrogen bonding interaction and ligand atom interactions with the protein residues. Ca-RMSD is calculated to measure the average change in the structure, conformation with reference to their respective first frames. The RMSD of the SB01990 shows that within 2 ns of simulation it rises up to 6 Å and further reached to 9 Å in 5 ns. However, shortly after reaching 10 ns it attains the equilibrium state at 7 Å and it is consistent throughout the end of the simulation. Compared with SB01990, SPB06836 has more fluctuations up to 10 ns but after that no significant change was observed and it is stable till the end of the simulation. The more compact and stable conformation was observed for KM06293 and without more fluctuation it attains the most stable conformation at 6 Å and this was maintained throughout the 30 ns simulation. This observation is consistent for these three compounds, which indicates that compounds that bind in the binding pocket of ORAI1 are much more stable. In contrast, RHO1882 shows increasing fluctuations till the end of the simulations and it indicates that the ligand has diffused away from its initial binding site during simulation. Figure 7(a) shows the Ca- RMSD of the protein complexes with the hit compounds SB01990, SPB06836, KM06293, RHO1882 and its average RMSD values were 7.54 Å , 9.23 Å , 6.85 Å and 9.28 Å , respectively (Figure 7b). However, the drastic change in the RMSD values indicates that large conformational changes occurred during the simulation to acquire the quite stable conformation. Information on the hydrogen bonding interaction throughout the simulation is monitored, because it is more important in drug designing for its significant role in drug specificity, metabolization and adsorption. The result of the hydrogen bonding interaction of each protein ligand complex is shown in Figure 8.

Ligand atom interactions with the protein residues
Protein-ligand interaction is the initiation step for every biological reaction. Therefore, it is important to encounter the changes in the ligand shape, surface potential, hydrogen bond interaction, hydrophobic bond interaction and ionic interaction between the protein ligand complexes throughout 30 ns simulation. Figure 9 (a-d) represents the number of interactions and it also points out the specific residues that interact with the ligand during the course of a trajectory, in which the residue shows light orange shade to dark orange shade representing the increasing number of contacts (H-bonds, hydrophobic, ionic, water bridges). On comparing ORAI1-ligand interaction, it is significantly observed that all three lead compounds (SB01990, SPB06836 and KM06293) remarkably interact with Glu106 with more number of contacts and this interaction is consistent throughout the 30 ns simulation. This provides the strongest evidence that the residue E106, a high calcium selective binding site region is strongly interacting with the ligand. To highlight the importance of this interaction in inhibition activity it is worth to mention the study carried out by Isabella Derler et al. (25), they characterized the CRAC channel inhibitory effect of the compound GSK-7975A through electrophysiological approach and it exhibited strong I CRAC inhibition. Following that, they studied the effect of these compounds in ORAI1 E106D pore mutant and interestingly found that action of GSK-7975A and 2-APB compounds on ORAI1 E106D pore mutant does not show I CRAC inhibition even at higher concentrations. To further clarify the effect of these  CRAC channel blockers in different altered pore geometry and to find out the precise binding region of these compounds, they mutated the outer vestibular region of ORAI1 D110/112/114A, but these mutations do not show much impact and substantially inhibit the CRAC channel as similar to wild-type ORAI1. From this they hypothesized that the compound GSK-7975A binds closely to the selectivity filter region and this mechanism might involve in CRAC channel inhibition. Intensely, to confirm their hypothesis, our study shows that the lead compounds KM06293, SB01990, SPB06836 have strong and direct interaction with the calcium selectivity filter region E106 and thus, it will have the greatest impact in blocking CRAC channel mechanism by altering the ion permeation and pore size. Therefore, E106 is majorly represented as a molecular determinant for the screened compounds to block the influx of Ca 2+ ions. In contrast, RHO1882 had interaction with E106 for the first few nanoseconds of simulation, but after that we did not observe any interaction with E106 till the end of the simulation. Apart from strong Glu106 interaction, compound SB01990 shows good interaction with Val 105 and Val107 throughout the end of the simulation and three side-chain hydrogen bond interactions with Tyr 115, Ser124, Thr 127 and hydrophobic interaction with Leu120 and Phe123. Compound SPB06836 shows very strong and stable hydrogen bond interaction with Val 105 and Val107. KM06293 shows the strong and stable interaction with Ala100, Met101 and Met104 only at the end of the simulation. RHO1882 forms strong hydrogen bond interaction with Lys 198, Lys 204 but it did not last after 20 ns and we observed that no residues were in contact with the ligand at the end of the simulation. As our aim of the study is to identify the novel CRAC channel inhibitor drugs that have the ability to block the Ca 2+ influx, the resultant hit compounds were particularly monitored for the Ca 2+ interaction changes that occur during the protein-ligand complex process. From the protein-metalligand interactions, it is significantly observed that, instantaneously after the ligand interaction in the binding pocket of ORAI1 resulted with the breakdown of the hydrogen bond  interaction tied between the Ca 2+ and the E106 residue. At the same time, this Ca 2+ ion forms the strong bidentate interaction with the ligand and as it can be seen in Supplementary Figure  S9(a-d), SB01990 forms strong Ca 2+ ion interaction with 98% and 99% stability and SPB06836 with 100% and 92% stability and KM06293 forms the Ca 2+ interaction with 74% and 82% stability and on the contrary, the RH01882 does not show any affinity toward calcium. Evidence from overall in silico simulation data suggest that these virtually screened lead compounds (SB01990, SPB06836 and KM06293) have made three major conformational changes on the binding site region of ORAI1 protein that will apparently alter the ion permeation and inhibit the influx of Ca 2+ : (i) due to the conformational changes in the acidic residue of the calcium selectivity region by the strong ligands interaction in the E106 region, it will greatly diminish the Ca 2+ selectivity of CRAC channels; (ii) based on protein-metal-ligand interactions, it is observed that these ligands show extremely strong affinity toward the Ca 2+ ions and thus it is presumed to behave like a chelating drug property that will specifically trap Ca 2+ ions which is prepared to enter the gate of the pore region that results in non-conducting CRAC channel; (iii) results from the induced fit docking energy, MD studies and free energy binding values lends support that these ligands have stable and stronger relative binding affinities toward the active site region of the ORAI1, which results in predominantly masking the highaffinity Ca 2+ binding site (E106). In addition, we have compared the existing reported known compounds with the compounds screened in this study. The results suggest that the scoring values of the above screened compounds are significantly better than the earlier reported ones (Supplementary Table S2). As a result, these compounds will primarily act as a barrier for the Ca 2+ to bind in the selectivity filter region, thus it would result in dramatic changes in the ion conduction channel property and the domino effect in blockade of the CRAC channel activity.

Conclusion
Aberration of CRAC channel function leads to several autoimmune diseases which makes modulators of CRAC channels one of the sought-after potential classes of lead molecules in the drug discovery pipeline for autoimmune disorders. In this study, we have screened and evaluated novel, potent ORAI1 inhibitors that specifically block the interactions between the Ca 2+ and ORAI1 using an in silico molecular modeling and docking approach. In addition, we also provide detailed structural information on the conformation and interaction of specific amino acids that are involved in inhibition of calcium influx. Our study starts with search of perfect 3D model of ORAI1 and further, structure-based virtual screening results with top screened compounds (SB01990, SPB06836, KM06293 and RH01882) with acceptable molecular modeling descriptors like docking energy, binding energy and ADME properties. Remarkably, from IFD, MD simulation and Prime MM/GBSA studies we did observe our screened ligands strongly bind with the E106 region and this alteration in the Ca 2+ selectivity region leads to lose its original activity that will greatly diminish the influx of Ca 2+ through CRAC channel. Consequently, this course of action blocks the influx of calcium entry, which in turn drastically affects the NFAT pathway results in subsequent reduction in differentiation of inflammatory T-cells and cytokine production. Taken together, we have shown that these three pyridinebased lead compounds (SB01990, SPB06836 and KM06293) with best pharmacokinetics ''drugs-like'' properties act as a novel and potent inhibitor chemical template that will effectively inhibit the influx of Ca 2+ by specifically targeting the ORAI1 active site. Further studies on synthesis and in vivo evaluations of these compounds could be a promising novel approach for the design and development of potent CRAC channel inhibitors and this drug could possibly be the new ''hero'' in the treatment of AID, including RA.