Thermochemical, PASS, Molecular Docking, Drug-Likeness and In Silico ADMET Prediction of Cytidine Derivatives Against HIV-1 Reverse Transcriptase

In recent, millions of people are living with the human immunodeficiency virus type 1 (HIV1), which causes acquired immunodeficiency syndrome. HIV-1 reverse transcriptase (RT) is one of the main viral targets for HIV-1 inhibition. Pyrimidine nucleoside derivative, 3′-azido-3′-deoxythymidine (AZT) is a highly active nucleoside inhibitor of human immunodeficiency virus type 1 (HIV-1) reverse transcriptase (RT). In this work, hydroxyl (-OH) groups of cytidine structure were modified with different aliphatic and aromatic groups to get 5 ́-O-acyland 2 ́,3 ́-di-O-acyl derivatives and then employed for molecular modeling, molecular docking, biological prediction, and pharmacological studies. Herein, we relate the optimization of cytidine and its acylated analogues applying density functional theory (DFT) with B3LYP/3-21G level theory to explore their thermochemical and molecular electrostatic potential (MEP) properties. Prediction of activity spectra for substances (PASS) indicated promising antiviral, anti-carcinogenic, and antifungal functionality of these cytidine esters compared to the antibacterial activities. To support this observation, their cytotoxic prediction and molecular docking studies have been performed against HIV-1 reverse transcriptase (RT) (PDB: 3V4I). Most of the molecules studied out here could bind near the crucial catalytic binding site, Tyr181, Ile94, Ile382, Lys374, Val381, Val90, and Tyr34 of the HIV-1 reverse transcriptase (RT), and the molecules were surrounded by other active site residues like Gln332, Trp406, Asn265, Gly93, His96, Pro95, and Thr165. Finally, these novel molecules were analyzed for their pharmacokinetic properties which expressed that the combination of in silico ADMET prediction, toxicity prediction, and drug-likeness had shown a promising result. The study discusses the performance of molecular docking to suggest the novel molecules active against resistance mutants of RT and/or recombinant strains of HIV-1.


Introduction
Nucleoside antibiotics have been under investigation for many years [1]. Some of the most clinically effective antiviral agents currently in use are purine or pyrimidine nucleoside analogs [2]. The cytidine ( Figure 1) analog KP-1461 is an anti-HIV agent that acts as a viral mutagen [3]. Alteration of hydroxyl (-OH) group at 3՛ and 5՛ position increases the antimicrobial activity of pyrimidine nucleoside and brings about some potent antimicrobial agents [4,5]. As a result of screening synthetic compounds for potential antimicrobial activity, a study reported that azidothymidine (AZT) has potent bactericidal in vitro activity against various members of the family Enterobacteriaceae [6]. Azidothymidine (AZT) is also one of the most popular nucleoside derivatives (antiviral drug) in which 3՛ hydroxyl (-OH) of thymidine was replaced by an azide group and now it's used worldwide for the treatment of HIV infection [7].
AZT suppresses the mode of reverse transcription, a ticklish phase in the life cycle of the virus. Edoxudine is another thymidine-derived antiviral drug, strongly working against herpes simplex virus [3]. Various cytidine derivatives modified at the base or ribose exhibit antiviral or antitumor activities [8]. Furthermore, alteration of hydroxyl (-OH) group in nucleoside derivatives uridine and cytidine also have potential antimicrobial activity [9][10][11]. Recently, some quantum chemical studies reported modified thymidine esters are also consisting of potent thermodynamic stability and pharmacological properties [12]. Reverse transcriptase (RT) of human immunodeficiency virus type 1 (HIV-1) is an important target for the design of new anti-HIV drugs. Small molecules that target the viral proteins cause reverse transcription and viral DNA integration is of great importance for pre-exposure HIV-1 prophylaxis [13][14][15][16][17]. The presence of HIV-1 resistance, associated with amino acid substitutions in RT, requires the development of novel approaches that explain the relationships between a particular conformation of HIV-1 RT and drug resistance. In this view, the potent antiviral efficacy of several cytidine derivatives (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15) with various aliphatic and aromatic chains was investigated by molecular docking against HIV-1 reverse transcriptase (RT) (PDB: 3V4I) along with the prediction of activity spectra for substances (PASS). Besides, attempts have been taken to optimize the acylated cytidine derivatives to predict their physicochemical and thermochemical behavior on the basis of DFT (B3LYP/3-21G). Finally, these novel cytidine derivatives were analyzed for their pharmacokinetic properties which revealed that the combination of in silico ADMET prediction, toxicity prediction, and drug-likeness calculation had shown a significant result. Hence, this study scintillates on the development of antiviral lead molecules from cytidine derivatives against HIV-1 in-silico tools and studies their pharmacokinetic and toxicity properties.

Materials and methods
To study the drug interactions with receptor proteins, molecular docking methods are the best suitable tools. The blind docking method employs a search throughout the whole surface of the protein molecule for binding sites. The following softwares systems were used in the present study: i) Gaussian 09, ii) AutoDock 4.2.6, iii) Swiss-Pdb 4.1.0 iv) Python 3.8.2, v) Discovery Studio 3.5, vi) PyMOL 2.3, vii) LigPlot þ v.2.2

Optimization of cytidine derivatives: DFT
In molecular drug design study, quantum mechanical (QM) methods have gained attention on the calculation of thermodynamic properties, molecular orbital features, dipole moment, and as well as interpretation of different types of interactions [18]. Molecular geometry optimization and further modification of all cytidine derivatives carried out using Gaussian 09 program [18]. Density functional theory (DFT) with Beck's (B) [19] three-parameter hybrid model, Lee, Yang, and Parr's (LYP) [20] correlation functional under 3-21g basis set has been employed to optimize and predict their binding properties.

PASS prediction
Web-based PASS (http://www.pharmaexpert.ru/PASSonline/index.php) was used for the calculation  of the biological spectrum of these cytidine derivatives [21]. Firstly, structures of the cytidine derivatives were drawn, and then, converted into their smiles formats by using SwissADME free web tools (http://www.swissadme.ch), which were used to predict biological spectrum using PASS online software. This program is designed to anticipate more than 4000 forms of biological activity including drug and non-drug actions and can be used to identify the most probable targets with 90% accuracy. PASS results are expressed as Pa (probability for active compound) and Pi (probability for inactive compound). Having probabilities, the Pa and Pi values vary from 0.000 to 1.000, and in general, Pa þ Pi 6¼ 1, since these probabilities are calculated freely. The activities with Pa > Pi are only considered as possible for a particular drug. The PASS prediction results were interpreted and used in a flexible manner, viz. (i) when Pa > 0.7, the chance to find the activity experimentally is high, (ii) if 0.5 < Pa < 0.7, the chance to find the activity experimentally is less, but the compound is probably not so similar to known pharmaceutical agents and (iii) if Pa < 0.5, the chance to find the activity experimentally is less, but the potentiality to find the activity structurally. As a result, the prediction of the activity of the spectrum of a drug is known as its intrinsic property.

Protein preparation and visualization
The 3D crystal structure of HIV-1 reverse transcriptase (PDB: 3V4I) (Figure3) was retrieved in pdb format from the protein data bank [22]. All hetero atoms and water molecules were removed using PyMol (version 1.3) software packages [23]. Swiss-Pdb viewer software (version 4.1.0) was employed for energy minimization of the protein [24] . Then optimized cytidine ligands were subjected for molecular docking study against HIV-1 reverse transcriptase (3V4I (Figure 3). In fine, molecular docking simulation was rendered by PyRx software (version 0.8) [25] considering the protein as a macromolecule and the drug as a ligand. AutodockVina was employed for docking analysis, and AutoDock Tools (ADT) of the MGL software package was used to convert pdb into a pdbqt format to input protein and ligands. The size of the grid box in AutoDockVina was kept at 87.8639, 71.8873, and 88.2770 Å for X, Y, Z directions respectively. After the completion docking, both the macromolecule and ligand structures were saved in. pdbqt format needed by Accelrys Discovery Studio (version 4.1) to explore and visualize the docking result and search the nonbonding interactions between ligands and amino acid residues of receptor protein [26]. The protein validation was checked by PROCHECK online server and it gives 93. 40  pubmed/1853201?dopt=Abstract). PDBsum online server was also used to check the validation of the main protease receptor with Ramachandran plot (Figure 2) which reveals that 93.40% residues in the allowed region.

Docking validation protocol
The docking validation was checked by extracting the co-crystallized ligand (3'-Azido-3'-Deoxythymidine-5'-Triphosphate) of the HIV-1 reverse transcriptase (PDB: 3V4I) and re-docking it into the same position. The lowest energy pose obtained on re-docking and the co-crystallized ligands were superimposed using PyMOL 2.3, and its root means square deviation (RMSD) was calculated between these two superimposed ligands. To validate the docking process, The RMSD must be within the reliable range of 2 Å [27,28]. It had been done to enhance ligand enrichment, which is necessary to test the docking procedure. Figure 4 represents the interfaces of HIV-1 reverse transcriptase that predicted from PDBsum. https://doi.org/10.37358/RC.21.3.8446

In silico pharmacokinetics ADMET and drug-like parameters prediction
To point out potential drug candidates, the ADMET properties were developed for the preliminary prediction of the pharmacokinetics and toxicity, and drug-like properties in the discovery drug process. In silico study gives a way to the accession of pharmacokinetic parameters (Adsorption, Distribution, Metabolism, Excretion, and Toxicity; ADMET) [29], its absorption in the human intestine, penetration of the blood-brain barrier and the central nervous system, the metabolism indicates the chemical biotransformation of a drug by the body, total clearance of drugs and the toxicity levels of the molecules. The Lipinski's rule of five properties was obtained from the SwissADME server (www.swissadme.ch /index.php) [30]. Prediction of the drug-likeness of the designed cytidine esters was also assessed by rule-based filters from Lipinski, Ghose, Veber, Muegge and Egan, and the synthetic accessibility difficulty scale was 1-10. The toxicity properties of ligands were assessed through pkCSM and SwissADME online server.

Results and discussions
In the present study, we have modified a total of fourteen cytidine derivatives with different aliphatic and aromatic chains (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15) and attempt to perform a quantum chemical study to realize the mode of their binding properties. Initially, we predicted biological activities using the PASS program. The observed activities were then rationalized by calculating their physicochemical (DFT method), cytotoxicity (predicted), molecular docking, and with the combination in silico ADMET, toxicity, and drug-likeness properties.

Computational evaluation of antimicrobial activities: PASS
We have predicted the biological spectrum using online PASS software (http://www.pharmaexpert. ru/PASSonline/index.php) of all the cytidine derivatives 2-15. The PASS results are designated as Pa and Pi which are presented in Table 1. It was manifest from prediction Table 1, for cytidine derivatives 2-15 showed 0.18 < Pa < 0.39 for antibacterial, 0.26 < Pa < 0.48 for antifungal, 0.53 < Pa < 0.70 for antiviral and 0.53 < Pa < 0.97 for anti-carcinogenic. This revealed that these molecules were more active against viruses and fungi as compared to bacterial pathogens. Attachment of additional aliphatic acyl chains (compound 2-5 and compound 8-13) increased antibacterial activity (Pa ¼ 0.392) of cytidine (Pa ¼ 0.395), whereas attachment of aromatic and bromo-substituted benzoyl groups did not improve reasonably (compound 6, 7, 14 and 15). In the case of fungal pathogens, the activity is improved to (Pa ¼ 0.484) of cytidine (Pa ¼ 0.345). But remarkable results were observed for antiviral activities where acyl chain derivatives are revealed to improve values than the halo-benzoyl derivatives. But compound 8 which has a tri-phenyl group exhibited the highest antiviral activity (Pa ¼ 0.708). These novel cytidine esters were found highly active against tumor protein producing gene TP53 (0.70 < Pa < 0.95). We have also predicted anti-carcinogenic property to support the activity against TP53. Thus, PASS prediction exhibited 0.53 < Pa < 0.97 for anti-carcinogenic which indicated that the cytidine derivatives were more potent as anti-carcinogenic agents than previous antimicrobial properties. Significantly, antibacterial, antiviral, and anti-carcinogenic properties of cytidine derivatives with saturated acyl chains (compound 2-5 and compound 10-13) were found more promising than the halo-benzoyl derivatives (6, 7, 14, and 15) except tri-phenyl derivative (8) [31].

Cell-line cytotoxicity (CLC) prediction
Web-based PASS (http://www.pharmaexpert.ru/PASSonline/index.php) was used to predict the cellline cytotoxicity of modified cytidine esters. Activity against lung adenocarcinoma, non-small cell lung cancer (stage-3), and pancreatic carcinoma have been predicted to suggest the maximum nontoxic bioactive drug molecule. It was evident from prediction Table 2, for cytidine derivatives (2-15) showed activities 0.65 < Pa < 0.73 for lung adenocarcinoma, 0.60 < Pa < 0.65 for lung carcinoma and 0.56 < Pa < 0.62 for pancreas carcinoma. There is a clear concept from the predicted data that these molecules have almost equal potentiality to work against these three cancer cells. All acyl chain substituted compounds (2-5 and 10-13) exhibited promising results than the triphenyl substituted (6), 4-t-butylbenzoyl substituted (7), Cl-benzoyl substituted (14) and cinnamoyl substituted (15). But derivative (8) which has a tri-phenyl ring exceptionally showed better activity than other aromatic substituted derivatives against tested cancer cells. We hope to conduct such studies for the more drug-related validation of these promising cytidine derivatives.

Molecular electrostatic potential map
Molecular electrostatic potential (MEP) is widely preferred as a map of reactivity that reveals the most suitable region for organic molecules to perform electrophilic and nucleophilic reactions of charged point-like reagents [32]. It aids to explore the biological recognition process and hydrogen bonding interaction [33]. MEP counter map provides a simple way to predict how different geometry could interact. The MEP of the title compound is obtained based on the B3LYP with basis set 3-21G optimized result and shown in Figure 5.  The prime significant of MEP is that it simultaneously displays positive, negative, and neutral electrostatic potential regions as well as molecular size, shape by color grading, and very useful in the study of molecular structure with physicochemical features relationship [34]. MEP was predicted to forecast the reactive sites for electrophilic and nucleophilic attacks of the optimized structure of cytidine 1 and its derivatives (6, 8, and 11). The different colors of electrostatic potential indicate different values. The potentiality of the attacking zone decreases in the sequence of blue ˃ green ˃ yellow ˃ orange ˃ red. The maximum negative area is displayed by red color where electrophiles can easily attack and the maximum positive area indicated by blue color which is suitable for nucleophilic attack. Moreover, the green color showed zero potential zones.

Docking validation study
In order to study the ability of docking algorithms to determine the conformation of the proteinbound ligand, re-docking of the co-crystallized Ligand was employed to validate the accuracy of the docking procedure. Figure 6(a) clarifies the superimposed view between the docked ligand conformation and the co-crystallized ligand conformation and the RMSD value is < 2Ǻ. The complex was then found to interact with the same amino acid residues compared to the ones reported in the present study. The bulky symmetric molecules can be exchanged in the binding site during docking, as the case in this investigation; the RMSD would be at a very high level. On the contrary, the small compounds can gain low RMSD easily even when placed randomly. Some reported studies [35 -37] have suggested a new benchmark for the quality of docking poses based on visual inspection. For visual inspection, Figure 6 (b) shows the 2D visualization of the interactions between a generated docking pose and the experimental ligand conformation. The results of this visual inspection show the same interactions as in the experimental binding mode, as observed in Figure 8 and Figure 10. This result indicates that the https://doi.org/10.37358/RC.21.3.8446 RMSD score alone is not a reliable parameter for the quality of docking poses for docking validation and the use of visual inspection as a new reference is essential. This partly verified the proficiency and accuracy of the docking procedure.

Molecular docking studies
Molecular docking is one of the most common methods used in structure based drug design to analyze the interaction between a small molecule and a protein at the atomic level [38]. By using the Autodock Vina 4.2.6, molecular docking was performed to suggest the best candidates among the selected cytidine derivatives 2-15 based on their binding affinities. All selected molecules were subjected to docking into the same binding pocket of HIV-1 reverse transcriptase (RT) (PDB: 3V4I) using similar optimized docking conditions to identify their binding mode. The results of the docking exploration showed that all cytidine derivatives, along with the parent molecule, gain binding affinities ranging from -4.5 to -8.3 kcal/mol. As shown in Table 3, compounds 2-7 and 9-12 showed low binding affinities compared to the parent compound, cytidine, while compounds 8, 13, 14, and 15 exhibited high binding affinity. These results revealed that modification of -OH group at position 2´,3´,5´of cytidine structure, along with an aromatic ring or an aliphatic chain molecule improved the binding affinity, while the insertion of halo-benzoyl groups like Cl and Br, made some fickleness in binding affinities. However, modification with either an aromatic ring or halogenated aromatic rings increased the binding affinity. The docked pose evidently showed that the drug molecules bind within the active site of the HIV-1 reverse transcriptase (RT) (3V4I) macromolecular structure (Figure 7). Figure 8 and Figure 10 shows that cytidine esters 7, 8, 14, and 15 (binding affinities -6.3, -8.1 -8.3, and -7.2 kcal/mol) bind firmly through hydrophobic bonds with the catalytic binding site Tyr181, Val35, Val90, Val179, Val381, Lys32, and Ile31, where, these residues exhibited alkyl, pi-alkyl pi-sigma pi-pi T-shaped and amide-pi stacked interaction. This pi-pi interaction revealed the tight binding with the active site. Again, Gln91, Gln161, Tyr144, Ile142 (shorter distance 1.88046 Å), Val90, and Thr165 which are the highly specific binding pocket of HIV-1 reverse transcriptase (RT) found to form a hydrogen bond. An electrostatic bond (Pi-cation) was also found with Lys172 and Lys374 at compound 7 and 15. It is evident from the structural contrast compounds 7-8, and 14-15 have an additional aromatic (4-tert-butylbenzoyl, triphenyl ring, halogenated ring, and cinnamoyl ring) substituent in the parent structure, indicating a high density of electron in the molecule leading to a comparatively higher binding affinity -6.   On the other hand, cytidine derivatives 2-6 substituted by various aliphatic chains (C10-C18) exhibited fluctuations in binding affinities (-6.0, -5.6, -4.5, -6.2, and -6.3 kcal/mol). Here, also observed that compound 6 which has an additional tri-phenyl ring showed a greater binding value. These derivatives showed some similar binding pockets with residues Tyr34, Tyr405, Asn265, Ala355, Arg356, Val261, Trp406, and Ile94 through both hydrogen and hydrophobic bonds. Most of the residue interacts at a distance of 2.23-2.90 Å. This result revealed that these potential inhibitors interact strongly with the mentioned residues HIV-1 reverse transcriptase. Compound 9-13, which are modified with a long aliphatic chain (hexanoyl, heptanoyl, lauroyl, and myristoyl), enhanced the binding affinities by specific recognition events at post-glycosidic linkage atomic positions two to three with Ile94 and Tyr181 and Tyr319 of the tyrosine gate. Interestingly, these derivatives specifically bind through the hydrogen bond with the specific pocket Gln91, Gln161, Gln182, Pro95, His96, Val90, Val381, Leu92, Ile382, Ile382, Trp383, and Lys172. Among all of these residues, Gln182 was found in close interaction 1.885Å.  Carbon atoms are most likely to make hydrophobic contacts with Ile94 which was observed in heptanoyl substituted cytidine esters 10. Two different interaction pi-donor hydrogen bonds and amide pi-stacked with His96 and Val381 were found in compound 13. Compound 7-9, 11 and 13-15 displayed the maximum π-alkyl, π-cation, pi-donor hydrogen, amide pi-stacked, pi-pi T-shaped and π-π interactions with the Tyr181, Ile382, Ile135, His96, and Val381 indicating the tight binding with the active site. https://doi.org /10.37358/Rev.Chim.1949 Rev. Chim., 72 (3)  It may suggest that Ty181 is considered as the principal component of the PPT, responsible for the accessibility of small molecules to the active site. Binding affinity and binding specialty are improved in the case of (2, 5, 6, 7, 8, 10, 13, and 14) due to significant hydrogen bonding. It was observed that modifications of the -OH group of cytidine (1) enhanced the π-π interactions with the residues of the active site while increasing their polarity resulted in the formation of hydrogen bonding interactions. The most prominent H-bonds were obtained for the compounds (5 and 12), forming with Ala355, Gln332, Thr338, Arg356, and Asn265. In contrast, compounds 2, 4, and 7 formed similar numbers of H-bonds with the same residues. Hydrogen-bonds execute a vital function in shaping the specificity of ligand binding with the receptor, drug design in chemical and biological processes, molecular recognition, and biological activity [39]. The hydrogen bond surface and hydrophobic surface of derivative 14 represent in Figure11. It was realized that the analyzed cytidine derivatives bound within some of the catalytic binding sites such as Valine (Val35, Val90, Val381, and val179), Tyr181, Tyr144, His96, Pro95, Lys101, Lys172, Ile31, Ile94, Asn135, Gln91, and Gln161 of the HIV-1 reverse transcriptase (RT), which is responsible for which causes acquired immunodeficiency syndrome. Figure 9 represents a two-dimensional LigPlot image of HIV-1 reverse transcriptase (RT) complex from PDBsum. The calculated binding affinities varied in the range of (-4.5 to -8.3 kcal/mol) suggesting the molecules can spontaneously interact within the binding site of HIV-1 reverse transcriptase. Among all the molecules, the inhibition activity of the compound (8, 14, and 15) was found to be the highest (-8.1, -8.3, and -7.2 kcal/mol). The results were summarized into a set of structural changes to be used in HIV-1 reverse transcriptase inhibitor design: cytidine esters gave an improved affinity and inhibitory potential, because of their relative flexibility combined with a favorable interaction with the active site.

Pharmacokinetic prediction
All newly designed cytidine derivatives have promising drug activities. Therefore, to confirm that the modified compounds are durable drugs, we employed the in silico pharmacokinetic parameters ADMET. The pkCSM online server [40] was employed to calculate in silico ADMET properties (Table  4). Absorption score less than 30% indicates poor absorption, here compounds 4-6 and 9-13 are displaying a value of 100% and 14-15 displaying >90% which reveals an excellent absorption in the human intestinal system. The volume of distribution (VDss) is thought high if the score is greater than 0.45. Moreover, bloodbrain barrier (BBB) and central nervous system (CNS) permeability standard scores (>0.3 to < -1 Log BB and > -2 to < -3 Log PS), respectively. For a given drug candidate a Log BB < -1 is poorly distributed to the brain, while Log BB >0.3 is the viable to cross BBB and LogPS > -2 thought to permeate the CNS, while Log PS < -3 are tough to move in the CNS [41]. It was noticed that most of the derivatives have the best significant probability to pass the barriers except derivatives 2, 3, and 8. The enzymatic metabolism ensures the chemical biotransformation of a designed drug in the body, which plays a key role in the transformation of drug compounds. In the body, drugs produce several enzymatic consequences, which plays an important role in conversion of drug compounds [42]. It is required to consider their metabolism of drugs, which may have different physicochemical and pharmacological properties. The cytochrome P450 (CYP450) plays a major role in drug metabolism because the major liver enzyme system is engaged in phase 1 metabolism. Some selective CYP genes CYP1, CYP2, CYP3, and CYP4 families are seen to be involved in drug metabolism, with CYP (1A2, 2C19, 2D6, and 3A4) causing the biotransformation of greater than 90% of drugs conducting phase I metabolism. residues Val35, Val90, Val381, Val179, Tyr181, Tyr144, His96, Pro95, Lys101, Lys172, Ile31, Ile94, Asn135, Gln91 and Gln161 of the HIV-1 reverse transcriptase (RT), which is responsible for which causes acquired immunodeficiency syndrome. These compounds showed several non-covalent interactions, such as hydrogen bonding, hydrophobic, and electrostatic interactions. These blind molecular docking analyses may provide a potential approach for the application of antivirals drugs as expected inhibitors HIV-1 reverse transcriptase. The docking validation process revealed that RMSD is in the standard range. Visual inspection exhibited very convincing results in the molecular docking validation process. In fine, these derivatives were analyzed for their pharmacokinetic properties which expressed that the combination of toxicity prediction, in silico ADMET prediction, and drug-likeness had shown promising results because the newly designed molecules have improved kinetic parameters and it maintains all drug-likeness rules as well as an interesting result in terms of biological activity. So, it could be concluded that most of the selected antivirals show promise and can be used to design effective antiviral agents against HIV-1. In this promising investigation, more drug-likeness in vitro and in vivo studies such as nontoxic concentration towards healthy cells may be conducted in the future.