Molecular Modeling Guided Drug Designing for the Therapeutic Treatment of Rheumatoid Arthritis

of all these receptors (IL-6, HLA-DR, and CD20) leads to the neurological, ocular, respiratory, cardiac, skin, and hematological manifestations that have been considered a potential therapeutic target for drug design. Various herbal, natural, and synthetic source inhibitors of interleukin-6 (IL-6), human leukocyte (HLA-DR), and CD20 were studied and reported previously. Reported inhibitors are compared to elucidate the best inhibitor for clinical trials, leading to the orally active drug. In this study, a computer-aided drug designing approach disclosed the potential inhibitors for all receptors based on their distinct binding a ﬃ nity. Moreover, drug suitability was carried out using Lipinski ’ s rule by considering the adsorption, distribution, metabolism, and excretion (ADME) of ligands. Results elucidated “ calycosin 7-O-glucoside ” and “ angeliferulate ” as putative ligands for IL-6 and HLA-DR, respectively. However, the pharmacokinetic properties (ADMET) revealed angeliferulate as an e ﬀ ete ligand for the biological system compared to calycosin 7-O-glucoside. Based on docking, drug toxicity pro ﬁ ling or pharmacokinetics, and MD simulation stability, this study highlights orally active therapeutic inhibitors to inhibit the activity of pivotal receptors (IL6, HLA-DR, and CD20) of RA in humans. After clinical trials, the resultant inhibitors could be potential therapeutic agents in the drug development against RA.


Introduction
Rheumatoid arthritis (RA) is considered a systemic inflammatory disorder that can cause significant disability and joint disease and enhance the mortality rate [1]. The persis-tent reduction in RA mortality has been noticed, but the prevalence has increased the economic burden [2]. RA is an autoimmune chronic inflammatory disease with articular and extra-articular manifestations. Autoimmune disease is a type of connective tissue disease that mostly overlays the joint and epitenon synovium [3]. RA leads to joint hardness, inflammation, and probable loss of function in any joint lined by a synovial membrane. However, small joints of hands and feet are most frequently damaged [4]. Other extra-articular tissues, including the vasculature, lungs, muscles, skin, and heart, are also affected by RA inflammation [3]. RA patients may suffer an increase in coronary infarction (heart attack), the risk of arterial sclerosis (hardening of the arteries), and stroke [5]. Some other issues associated with RA could include endocarditis, inflammation in the cardiac valve, and fibrosis [6]. The prevalence of RA in Pakistan (1.6 per 1,000 people) is less compared to that in Bahrain (2.0), Djibouti (2.2), Lebanon (2.1), Qatar (1.8), Somalia (1.9), and Tunisia (1.8) [7]. However, the frequency of RA in North America and Europe [8,9] is comparatively higher than that in Japan [10], China [4], Indonesia, Malaysia, the Philippines, and even rural Africa [11]. The prevalence of RA is relatively constant at 0·5-1% worldwide [12].
A growing number of disease pathways have been discovered and reported [13][14][15][16], leading to potential target proteins for drug design [17][18][19][20]. Several studies illustrate RA as an inherent immune response [21], which could be genetic due to "shared epitope." The shared epitope is a motif sequence of five amino acids (70-74 residues) present on the HLA-DRβ chain [22][23][24], associated with severe rheumatoid arthritis in most cases [25]. In RA pathogenesis, the interleukin-6 receptor (IL6R) also plays a pivotal role in triggering RA [26]. Since RA mainly attacks synovial articulations and induces systemic inflammation over time, the insights into the pathogenesis of RA have led to the development of B cell-driven therapies as promising new therapeutic targets in autoimmune disorders. Besides autoantibody-secreting cells, B cells may have several important roles in the immuno-pathogenesis of RA by producing proinflammatory cytokines [27][28][29]. The anti-CD20 rituximab antibody has been used to treat patients with RA-refractory Tumor Necrosis Factor (TNF) [30]. In the pathogenesis of rheumatoid arthritis, many receptors and antigens are involved, and the type of inflammatory response may also be different depending upon the initiation of autoimmune response [31,32]. Several herbal, natural, and synthetic source inhibitors of interleukin-6 (IL-6), human leukocyte antigen (HLA-DR), and CD20 were studied and reported previously [3,[33][34][35]. However, a lack of comparison of the available inhibitors can lead to selection bias in a given biological system. To address this gap, previously reported RA inhibitors have been compared to identify the optimal choice molecule that can potentially lead to orally active drugs.
Computer-aided drug designing (CADD) [36,37] is the most common and widely used method in the modern age for drug discovery [38]. In this study, three receptor proteins (IL-6, HLA-DR, and CD20) are docked in silico to elucidate efficient ligands against these receptors. Molecular docking is aimed at predicting the experimental binding affinities of inhibitors (ligands) within the active site of the target receptor (protein). Thus, molecular docking discloses the interaction of ligand with active sites of pro-tein, whereas drug-likeness elucidates the potential of ligand for suitability in a biological system [39]. Hence, we predicted the drug-likeness of each ligand parallel to molecular docking analysis for explication of putative orally active ligand against rheumatoid arthritis. Physicochemical and molecular characteristics frequently have contradictory effects on pharmacokinetic and pharmacodynamic processes as well as medication safety. In tandem with target affinity, the current approach in drug discovery is to examine ADMET characteristics. In order to help the medicinal chemist in prioritizing drug candidates, the idea of "druglikeness" delineates acceptable limits of fundamental features expressed as simple rules of thumb. Absorption, distribution, metabolism, excretion, and toxicity (ADMET) of chemical entities play important roles in drug discovery and development. A good drug candidate should not only be effective against the therapeutic target but also have adequate ADMET characteristics at a therapeutic dosage. During clinical development, ADMET profiling may assist to reduce possible hazards [40]. That is why we finalized therapeutic potential agents based on their ADMET properties as well.

Characterization of Target
Proteins. The physicochemical properties of all receptor proteins (IL-6, HLA-DR, and CD20) were predicted using ProtParam online server [41]. ProtParam works on the Edelhoch method [42] to identify weight value, hydropathy values for extinction coefficients, instability index (II), and GRAVY value (grand average of hydropathy value) of proteins.

Tertiary
Structures of Target Proteins. The crystal structures of target proteins were retrieved from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB-PDB) [41]. The literature stated that IL-6, CD20, and HLA-DR are significant proteins involved in RA [12]. After analysis, IL-6 (PDB ID: 1il6) and CD20 (PDB ID: 2osI) were finalized based on high-resolution power and the lowest R value (R value determines how well the simulated diffraction pattern matches the experimentally observed diffraction pattern). Their tertiary structures are shown in Figure 1.
However, the tertiary structure of HLA-DR of Homo sapiens was not available in PDB, which was predicted through homology modeling. The primary sequence of HLA-DR Homo sapiens (UniProt ID: P01911) was retrieved from UniProtKB [43]. MODELLER [44] software was used to predict the 3D structure of HLA-DR through homology modeling. Predicted structure validation is necessary before doing docking to check the quality of the structure; therefore, the Ramachandran plot [45] was used for structural stereochemistry validation of our target proteins. The predicted tertiary structure of HLA-DR is shown in Figure 2. 2.3. Screening of Chemical Inhibitors. The initial step for pharmacophore designing is to identify the novel inhibitors against target receptors. The chemical compounds against active sites of IL-6 [46], CD20 [30,47], and HLA-DR [48] 2 Cellular Microbiology have known medicines used for RA treatment. Therefore, natural and synthetic derivatives of compounds having inhibitory potential against our target proteins were retrieved, and structures of these ligands were downloaded using PubChem [39,49,50]. Ligands selection was based on using keywords of our targets such as "CD20", "HLA-DR", and "IL-6" in PubChem, and we selected those chemical entities having biological inhibitory assays against required sites which are helpful in treating RA disease. PubChem retraces compounds based on their chemical formula and physicochemical properties [49]. In the current study, we mainly focused on three major receptors (HLA-DR, IL-6, and CD20) of RA [12]. Several herbal, synthetic, and animal source inhibitors (ligands) used in the study are listed in Table 1. In addition, X-ray crystalized tertiary structure of antibody "ofatumumab" (PDB ID: 6Y92) was retrieved from the Protein Data Bank (PDB) [51].

Active Site Prediction.
InterPro [52,53], an online webbased tool, was used to investigate the active sites of target proteins for molecular docking against reported ligands.
2.5. Molecular Docking Analysis of Inhibitors. Virtual screening (VS) is one of the standard steps for drug discovery before wet-lab experiments. This process involves the estimation of the binding affinity of the drug candidate towards a target protein. Virtual screening also exhibits possible binding modes against active sites of target proteins. Most prominent drug candidates show promising binding affinity towards target protein and can be screened out using High-Performance Computing (HPC) infrastructure [54]. A molecular operating environment (MOE© 2009.10 version) was selected for molecular docking to identify a putative ligand for receptor inhibition [50,55]. MOE© prioritizes efficient binding geometries based on S-score's numerical value for docked poses among various available resources listed elsewhere.
Molecular docking [36,56] analysis was implemented to finalize the hits for pharmacophore designing. For this, the 2D structures were obtained from the PubChem [57] using their respective PubChem ID, and energy minimization was done using the Marvin Desktop Suite, ChemAxon, Hungary [58]. The energy minimized structures were then converted to Microsoft Access Database (MDB) formats and further docked with the target structure of MOE©. Binding energy, molecular interactions (at the active site), and structural confirmation were the prime parameters to finalize the inhibitors. The target proteins (IL-6, CD20, and HLA-DR) were prepared by removing hetero-atoms and adding polar hydrogen and Kollmann charges using MOE©. The average value was found for binding energy as an S-score, and all the molecules above the average value (-6 kcal/mol) were chosen for pharmacophore designing [59]. This 3D structure was converted to MDB format, and a total of 30 confirmation poses of each ligand were created against respective receptor proteins. The first ten results (based on pose confirmation and binding energy) were analyzed using MOE© ligand-receptor interaction and pocket map tools.
2.6. Enrichment Analysis of the Docking Library. The chemical properties of selected active compounds were determined by enrichment analysis. Useful decoys (DUD) were narrowed down for enrichment analysis based on the chemical properties of compounds, as explained elsewhere [39]. A  subgroup of commercialized chemical compounds, 31,200 in number, was randomly selected from a broader dataset. The same molecular weighted DUD distribution was chosen to narrow down the active compounds further and alleviate the notable tendency of the scoring component to support large compounds. The set of decoys was refined by filtering functional groups and cutting out both rotational bonds and molecular weight to guarantee that separation was mainly based on the scoring function. When we combined the randomly chosen subset with 450 receptor-specific ligands and excluded redundant structures, a total of 3100 compounds were obtained and prepared for additional docking investigations.

Molecular Property Distribution of Arbitrarily Chosen
Compounds. Before conducting virtual screening, the molecular properties viz polar surface area [60], rotary bonds, molecular weight (MW), and logD of selected compounds against particular targets were already computed. In addition, hydrogen bond acceptors [61] and donors (HBD) were also calculated. This was needed to decide an evident foundational inconsistency between the known active compounds assortment and the arbitrarily picked compounds. However, drug-likeness was applied for the underlying irregular arrangement of the compounds. As a result, no apparent discrepancies were found between known active substances and randomly chosen compounds in the chemical properties around the predefined ranges.
2.8. Pharmacokinetic Properties and Toxicity of Ligands. The SwissADME© tool was utilized to predict the drug-likeness of our inhibitors [62]. Different drug-likeness parameters, including excretion, metabolism, distribution, and absorption of the compounds for the biological system, were evaluated [63]. The potential ligands were further assessed based on absorption, distribution, metabolism, excretion, and toxicity (ADMET) to reduce failure in drug discovery [64]. AdmetSAR© is a web-based tool containing information on toxicity, carcinogenicity, and compliance of medicine with Lipinski's regulations. The finalized ligands' SMILES were then posed for toxicity tests in the AdmetSAR© programmer [65].
2.9. Molecular Dynamic Simulation. Based on molecular interaction and visual evaluation of docking data, the topranked complexes (calycosin 7-O-glucoside-IL6, angeliferulate-HLA-DR, and ofatumumab-CD20) were chosen for the Molecular Dynamic (MD) simulation investigation. The MD simulation was run on a supercomputer cluster with HPE-DL385 [66], Gen10, AMD EPYC 64 cores, RAM 128 GB, 500 GB SSDs, and NVIDIA Tesla T4 16 GB Accelerator on Schrodinger's Desmond module [61]. A water-soaked solvent solution was used to make predictions. Desmond's system builder tool was used to design the water-soaked solvated system. The TIP3P water model was looked at to see if it may help address the problem. A box with periodic boundary conditions and a buffer distance of at least 10 Å from the protein's outer surface was used to produce the orthorhombic simulation. To neutralise the system, a suitable amount of counterions was introduced. The isosmotic condition of the simulation box was maintained by adding 0.15 M NaCl. A predefined equilibration method was followed before the simulation's production run. At a temperature of 300°K and a pressure of 1.013 bar, the MD simulation was run. The Simulation Interaction Diagram was used to examine the MD simulation trajectory.  Table 2. IL-6 has a molecular weight of 20980.97 Daltons and isoelectric pH 6.22. The GRAVY value -0.498 indicates that IL-6 has hydrophilic property. The instability index of 59.45 predicted that IL-6 is a stable protein containing an aliphatic index of 84.43, as shown in Table 3.

Results and Discussion
CD20 has a molecular weight of 33077.32 Daltons and isoelectric pH 5.04. The GRAVY value -0.048 exhibits the hydrophilic nature of IL-6. The instability index of 64.69 predicted IL-6 as a stable protein containing an aliphatic index of 91.62, as shown in Table 4.

3.2.
In Silico Molecular Docking. Molecular docking of IL-6, HLA-DR, and CD20 target proteins was accomplished by using MOE©. The results verified that adopted inhibitors were inside the active pocket sites of the target proteins with possible interactions. However, some have shown an important binding affinity with active site residues such as calycosin 7-O-glucoside, ononin, angeliferulate, and ofatumumab. An inhibitor with a lower value (high in negative) of free binding energy is supposed to establish strong interaction  Table 5. We docked 15 synthetic compounds against IL-6 and HLA-DR and one antibody against CD20 antigen. Thus, MOE© generated 200 possible poses with numerous S -scores (binding energy score). Out of 200, the top pose of each ligand is enlisted in Table 5. Out of these top poses, calycosin 7-O-glucoside, angeliferulate (S-score = −10:9376), and ofatumumab antibody (S-score = −10:3737) have shown maximum binding energy affinity (S-score = −10:751) against IL-6, HLA-DR, and CD20, respectively. A threshold for binding energy score was set as -6.0 kcal/mol [59], and the ligand has to follow at least three conditions out of Lipinski's rule of five. Docked poses of high binding affinity ligand-protein complex are shown in Figure 3. 50 and Lipinski's rule of five) were implemented to differentiate between active and   inactive ligands. Those compounds were considered active in which IC50 values ranged between 1 and 10 μM, while the rest were declared inactive. According to Lipinski's rule of five (LRo5), the chemical compounds violating more than one condition of Lipinski's rule were supposed to be inactive for biological implications and vice versa. The graphical representation of enrichment analysis of active and inactive ligands is shown in Figure 4.

Drug Liability and Toxicity
Evaluation. The ligand's antagonistic interactions with receptor protein do not ensure the inhibitor is an active oral drug. Consequently, ADME and drug-likeness play a critical role in inspecting the credentials of ligands for biological systems. Egan's BOILED-Egg method and LRo5 were utilized to evaluate ligands' pharmacokinetics. LRo5 relies upon edge upsides of four physicochemical properties: molecular weight ðMWÞ ≤ 500 g /mol, lipophilicity ≤ 5, hydrogen bond donors ðHBDÞ ≤ 5, and hydrogen bond acceptors ≤ 10. The radar diagram demonstrates that both selected ligands satisfied the LRo5, which indicates the suitability of these chemical compounds (ligands) for bioavailability Figure 5. However, immunoassay is necessary to cross-check our findings related to the ofatumumab monoclonal antibody.
Calycosin 7-O-glucoside showed human intestinal absorption, obeyed Lipinski's rule, behaved as noncarcinogenic, and led the liability of human ether-a-go-go-related gene (hERG) channel inhibition, estrogen receptor binding (ERB), androgen receptor binding (ARB), and thyroid receptor binding (TRB). These liabilities validate our identity hit (calycosin 7-O-glucoside) as an orally active therapeutic agent against RA. In contrast, angeliferulate did not show human intestinal absorption and exhibited liver toxicity. Calycosin 7-O-glucoside is considered a better therapeutic compound than angeliferulate due to its high binding affinity, good physicochemical properties, and liability with body channels, as shown in Table 6.
There are three general classes of medications typically utilized to treat rheumatoid joint pain: nonsteroidal calming specialists (NSAIDs), corticosteroids, and disease-modifying antirheumatic drugs (DMARDs) [20]. NSAIDs and corticosteroids have a short-span onset of action, while DMARDs have delayed effects, and it can take several weeks to show a clinical impact. DMARDs include methotrexate, sulfasalazine, leflunomide, etanercept, infliximab, adalimumab, certolizumab pegol, golimumab, abatacept, rituximab, tocilizumab, anakinra, and antimalarials. Other immunomodulators are occasionally used, including azathioprine (Imuran) and cyclosporine [67]. Since cartilage harm and hard disintegrations typically happen within the initial two years of the disease, DMARDs are preferred from the beginning and throughout sickness, usually when an analysis is affirmed. Pain-relieving drugs are also used at the start as they are concurrently in diminishing torment until DMARDs produce results [68,69]. Though DMARDS significantly ameliorate the disease, they have multiple common side effects like abdominal pain, chills or high fever, dizziness, hair loss, headache, light sensitivity, itching, hepatic disorders, and low blood counts. Besides these, dry cough, fever, or difficulty in breathing may result from lung inflammation requiring special care in COPD patients [70].
With all the advancements in DMARD research, there is still a need for a greater understanding of its underlying mechanisms, as by modifying the mechanism of drugs' mode of action towards studying the exact mechanism of disease at molecular levels, we can facilitate the development of new drugs and produce a paradigm shift in treatment [71].
3.5. Molecular Dynamic Simulation. Docking alone cannot provide full insight into the binding mode, stability, and dynamics of proposed ligands. As a result, we used the Desmond module of Schrödinger to run MD simulations for different nanosecond frames based on the stability point of  The protein-ligand complex combined RSMD during MD simulation remains stable throughout the simulation time frame of 200 ns. From this stability, it is elucidated that with respect to time, this docked complex is structurally stable due to hydrogen bonds.
While the MD simulation RMSD trajectory of HLA-DR and angeliferulate complex has revealed that docked complex is dynamically stable with RMSD value of 1.5 Å to 3.5 Å (RMSD change < 3 Å) throughout a 150-nanosecond time period (Figure 7). The combined RMSD of the protein-ligand complex during this MD simulation also remained stable throughout the simulation time frame of 150 ns as combined trajectories did not show any notable fluctuations (only shown at a short time frame from 82 ns to 86 ns) in trajectory or RMSD change. This stability elucidated that with respect to time, the docked complex is structurally stable after making hydrogen bonds and other interactions. The combined RMSD of the protein-ligand complex during MD simulation remained stable throughout the simulation time frame of 150 ns. CD20 and ofatumumab MD simulation RMSD trajectory plot elucidated this docked complex interaction stability throughout 100 ns simulation time with 1000 frames of snapshots. The CD20-ofatumumab complex RMSD plot has shown that the trajectory of RMSD deviates just 1.1 Å from the complex initial position at the docked time (Figure 8), and this deviation is acceptable and reflects the stability of the docked complex with respect to time. The shift in the protein's RMSD of the order of 1.8 Å-2.7 Å during the simulation suggests that the protein underwent moderate conformational change throughout the simulation.
Based on molecular docking, MD modeling, and toxicity profiling, calycosin 7-O-glucoside, angeliferulate, and   9 Cellular Microbiology ofatumumab were discovered to be effective and selective potential agents to treat RA disease. Concerns have been raised about the risk of the usage of these agents; for example, it has been shown to increase the risk of cardiovascular disease and the severity of infections. Recent research suggests, however, that administration of these agents can reduce the incidence of cardiovascular disease and has a neutral infection risk. During clinical studies, it is advised that the other stated medical condition be considered to better examine its potential as a medication against HLA-DR, IL-6, and CD20 targets.
In silico comparative analysis is recognized as an initial step for insight into ligand-receptor complexes and drug targets that lead to a better selection process of potential disease inhibitors. Our findings show "calycosin 7-O-glucoside" as a potential orally active drug due to its better pharmacokinetic profile and significantly higher binding affinity against HLA-DR (-10.0867) and IL-6 (-10.751).

Conclusion
Briefly, combined results of molecular docking, Lipinski's rule of five, and drug toxicity assessment exhibited calycosin 7-O-glucoside (S-score = −10:751) to be of most significance among selected ligands against RA. Angeliferulate (S-score = −10:9376) showed an effective binding score against target proteins but less human intestinal absorption. Moreover, the hepatoxicity of both compounds predicted adverse effects of these compounds on the liver after long-term usage. In contrast, the monoclonal antibody ofatumumab exhibited potential therapeutic efficacy against receptor (CD20) of RA. In the light of current findings, calycosin 7-O-glucoside has all those characteristics which should be considered during drug discovery and clinical trials.

Data Availability
The corresponding author will provide all data produced during the current study on a reasonable request.

Conflicts of Interest
The authors have declared no conflict of interest.

Authors' Contributions
HAR and AAB conceptualized, designed, interpreted, and supervised the study and contributed to the final draft of the manuscript. SAA and MI performed molecular docking analysis. MH, MU, and QUA carried out drug-likeness methods and interpreted the ligands' pharmacokinetics. MM, MI, SAA, and MH contributed to writing the manuscript in discussion with all other authors.