In silico investigation of the interactions of certain drugs proposed for the treatment of Covid-19 with the paraoxonase-1

Abstract Coronavirus disease 2019 (Covid-19) has caused one of the biggest pandemics of modern times, infected over 240 million people and killed over 4.9 million people, and continues to do so. Although many drugs are widely recommended in the treatment of this disease, the interactions of these drugs with an anti-atherosclerotic enzyme, paraoxonase-1 (PON1), are not well known. In our study, we investigated the interactions of 18 different drugs, which are claimed to be effective against covid-19, with the PON1 enzyme and its genetics variants L55M and Q192R with molecular docking, molecular dynamics simulation and free energy calculation method MM/PBSA. We found that ruxolitinib, dexamethasone, colchicine; dexamethasone, sitagliptin, baricitinib and galidesivir, ruxolitinib, hydroxychloroquine were the most effective compounds in binding PON1-w, PON1L55M and PON1Q192R respectively. Mainly, sitagliptin, galidesivir and hydroxychloroquine have attracted attention by showing very high affinity (<-300 kJ/mol) according to the MM/PBSA method. We concluded that the drug interactions should be considered and more attention should be paid in the use of these drugs. Communicated by Ramaswamy H. Sarma


Introduction
Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2), a positive sense and single-stranded enveloped RNA virus, causes coronavirus disease 2019 . This virus enters the oral and nasal cavities, binds to the receptors on the epithelial cell membrane and can also travel further down the airway into conductive airways. Thus, the virus entering the cell through receptor activates humoral and cellular immunities by releasing its RNA to the cytoplasm (Chams et al., 2020;Ko car et al., 2021). The immune response is strongly modulated by oxidative stress and inflammatory processes. High production of free radicals by the immune cells triggers oxidative stress and oxidizes biomolecules including lipids, proteins, nucleic acids (Iddir et al., 2020). The oxidation products are potent activators of scavenger receptors and induce further inflammation. Cytokine storm causes the dyslipidemic conditions (Sorokin et al., 2020). Although blood levels of total cholesterol, low density and high density lipoprotein cholesterol (LDL-C and HDL-C, respectively) decrease mainly in severe Covid-19, high levels of membrane cholesterol contributes SARS-CoV-2to entry into host cells (Ko car et al., 2021). Dependent on inflammation-induced oxidative stress, reverse cholesterol transport deteriorates (Sorokin et al., 2020). A reverse cholesterol transporter HDL, exerts also some other functions by its antioxidant, anti-inflammatory and antiviral properties. Therefore, HDL functionality is important against Covid-19 (Cho et al., 2021).
Paraoxonase 1 (PON1, EC 3.1.8.1, an aryldialkylphosphatase) is mainly carried by HDL in circulation. It is a calcium ion-dependent enzyme that has capability to hydrolyze lactose, thiolactone, aryl esters, paraoxon and different organophosphates. PON1 activity differs according to its genetic polymorphism. Its most important polymorphisms are PON1Q192R (glutamine or arginine at position 192) and PON1L55M (methionine or leucine at position 55). The Q isoform gives much better protective effect against LDL oxidation than the R isoform (Kotur-Stevuljevi c et al., 2020). PON1 in the cells increases cholesterol efflux capacity and affects the redox status (Raz et al., 2020). It shows the protective roles against cardiovascular diseases by its atheroprotective properties such as protecting the cell membrane, LDL and HDL from oxidation due to oxidative stress and breaking down lipid peroxides in LDL and HDL (Kotur-Stevuljevi c et al., 2020;Mar ın et al., 2019). In addition, it activates endothelial nitric oxide synthesis and thereby in vasodilatation and prevents platelet aggregation (Mar ın et al., 2019). This function is important mainly in severe Covid-19 disease in which blood hemostasis may be impaired due to endothelial cell dysfunction, and hypercoagluation and thromboembolic disorders may occur (Chams et al., 2020). Moreover, it shows antiviral properties by contributing HDL to lower cholesterol levels in lipid rafts (Ko car et al., 2021). Furthermore, HDL with high PON1 exerts antiviral activity by suppressing the replication of SARS-CoV-2. Cho et al. (2021) showed that native HDL level was associated with high paraoxonase activity and with the decreased cytopathic effect. Because of the important functions described above, we think it is important to identify molecules that can affect PON1 activity.
Molecular docking is an inexpensive computer-based computational process used to determine the interaction of selected ligands with amino acids in protein structure in low energy conformation. In this method, which is widely used in the pharmaceutical industry, a large number of candidate molecules are screened and eliminated through various docking protocols, including Autodock (Shah et al., 2020). In addition to the molecular docking approach, molecular dynamic (MD) simulation and free energy calculations are performed to increase the calculation accuracy. The molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) approach to calculate binding free energy is often preferred due to its relatively high accuracy and low cost (Wang et al., 2019).

Ligands
The potential substrates of PON1 used in this study listed in Table S1, with their molecular structures and IUPAC names .
Potential drugs in medications of Covid-19 chosen for our study were given in Table S2 (Chams et al., 2020;Cho et al., 2021;Elekhnawy et al., 2021;Keni et al., 2020;Nissen et al., 2021). All drugs were obtained from the PubChem ligand database . ChimeraX software was used to demonstrate protein-ligand interactions in three dimensions. (Goddard et al., 2018). Discovery studio visualizer 2020 was used to draw two-dimensional protein-ligand interaction diagrams (BIOVIA, Dassault Syst emes). Graphpad prism 6 software was used for plotting (San Diego, California USA, www. graphpad.com).
Modeling of PON1-w, PON1Q192R and PON1L55M polymorphic structures As the PON1 protein, 1.99 Angstrom (Å) resolution nonmutant structure, 3SRE encoded X-RAY structure was preferred (Ben-David et al., 2012). The X-ray crystal structure form of the PON enzyme (3SRE) obtained from the protein databank is prepared by separating it from ligands and solvents with the Chimera program (Pettersen et al., 2004). MODELLER software was used for homology modeling (Webb & Sali, 2014). Flexible missing residues in the range 72-81 in the 3SRE model were added by loops/refinement tools in MODELLER. Non-terminal missing structure option was selected, and 20 models were generated for each wild and mutant structures. The amino acid sequence of the human serum PON enzyme was used as a template. The sequence was obtained from the Uniprot protein sequence database (code: P27169) (The UniProt Consortium, 2017). Since the PON1Q192R and PON1L55M polymorphic structures of the PON enzyme will also be analyzed, the amino acids in the 55 th and 192 th regions have been changed in accordance with the PON1Q192R and PON1L55Mpolymorphisms. With the MODELLER program, 20 models were obtained for each template (PON1-w, PON1Q192R and PON1L55M), and tested using Molprobity server for determination of protein structure quality and the highest scored model was preferred (Chen et al., 2010).

Molecular dynamic simulation
MD calculations took place in two stages. In the first stage, 10 ns MD simulations of PON1-w and its variants were performed, and in the last ns, snapshots were taken at intervals of 100 ps and molecular docking calculations were carried out. The second stage of the MD simulation was started with the snapshot that had the best correlation with the experimental data obtained from the study of Draganov et al. (2005). In the second stage, 100 ns MD simulations of protein-compound complexes were carried out with top-scored three Covid-19 drugs for each protein.
First stage MD simulations before molecular docking MD calculations were carried out with GROMACS 4.5.5 software in 4 nodes, each with 24 cores (AMD Opteron), in  computer clusters on TRUBA (Pronk et al., 2013). In the study, the 'Leap Frog' algorithm was preferred as the integrator and all steps were performed with 2 fs intervals. MD simulations were run at periodic boundary conditions (PBC) in the form of 'Rhombic Dodecahedron' with 1 nm between the edge of the protein and solvent environment boundary. A pre-minimization process was carried out in a vacuum environment to eliminate clash contacts in the structure. 0.15 mM Na-Cl were added to neutralize the system. In the minimization phase, a total of 5000 steps steepest descent minimization algorithm was prepared. After minimization, the system was equilibrated with MD simulation at 100 ps in the canonical (NVT) ensemble and 1 ns in the isothermal-isobaric (NPT) ensemble. With the LINCS constraint algorithm, all bonds and atomic positions are fixed at 1000 kJ/mol force. In the NVT phase, the temperature was set to 300 K and the v-rescale algorithm was used. In the NPT phase, the temperature and pressure were adjusted as 300 K and 1 atm, with v-rescale and Parrinello-Rahman coupling algorithms, respectively. In the production phase, the parameters in the NPT phase were used and position restrictions were removed and the 10 ns molecular dynamics simulation for each structure was carried out, and the data were obtained (Marto n ak et al., 2003). All trajectories were saved in the 10 ps intervals for analysis at later stages. In order to understand that the simulation was carried out successfully, root mean square deviation (RMSD) calculations were made using the 'g_rms' tool. The g_rmsf tool was used to report the fluctuations of each residue in the protein. The radius of gyration, a measure of the compactness of the protein, was calculated with the g_gyrate tool. Amber ff99SB-ILDN was chosen as the 'forcefield' as it is one of the most accurate force fields (Lindorff-Larsen et al., 2010). As the water model, TIP3P was preferred as it is one of the most accurate models in the literature. The contents of each polymorphic structures in MD simulations are shown in Table 1. Our study systems contained about 18-19 thousand water molecules, 70 Na þ , 56 Cland 2 Ca 2þ (cofactor of the PON1 protein) ions. The mean molecular weight of the proteins is 39.75 kDa.

Second stage MD simulations after molecular docking
MD simulations were run at periodic boundary conditions (PBC) in the form of 'Rhombic Dodecahedron' with the 'Leap Frog' algorithm with the GROMACS 2020.3 software package with 2 fs time intervals (Abraham et al., 2015). Thanks to its spherical structure, the Rhombic Dodecahedron form increases the computational performance up to 40% compared to the rectangular box (Abraham et al., 2015). The corner and cube distance of the protein-ligand complex was set at 1.2 nm. Neutralization was performed by the addition of 0.15 mM Na-Cl. 50,000 steps of energy minimization was performed according to the Steepest Descent algorithm. MD simulation with a length of 300 ps was performed in the NVT stage. With the LINCS constraint algorithm, all bonds and atomic positions are fixed at 1000 kJ/mol force. 1 ns long MD simulation was performed at the NPT stage. With the LINCS constraint algorithm, only bonds are restricted (Hess et al., 1997). Atomic positions are released. The Nos e-Hoover algorithm was used as a temperature coupling algorithm and the Verlet algorithm was used as a cut-off scheme and the temperature was set to 300 K (Evans & Holian, 1985;Grubm€ uLler et al., 1991). As the pressure coupling algorithm, the Parrinello-Rahman algorithm with isothermal compression was used and the pressure was set to 1 atm (Marto n ak et al., 2003). Unlike the NPT phase in the production phase, all restrictions were removed and 100 ns MD simulation was performed. Van der Waals and short-term electrostatic interactions were cut at 10 Å and 'Particle-mesh-Ewald' was used as long-range electrostatics (Darden et al., 1993).

Molecular docking
Studies in the literature have shown that Y71 and R292 residues have flexible sides which act as gates in protein-ligand interaction (Ben-David et al., 2012;Fairchild et al., 2011;Hu et al., 2009). Therefore, in this study, these residues were adjusted flexibly. The residues other than Y71 and R292 were considered rigid. Autodock 4.2.6 software was used for molecular docking and the docking procedure (Ben-David et al., 2012;Morris et al., 2009). In AutoDock 4.2.6, an average of 39 A o cube was applied for the grid volume by targeting the active center of the PON1 enzyme for docking process. Electrostatic area was calculated according to 1 A o grid area. Ca atoms in the catalytic center were set to þ1 charge. A total of 100,000 step minimization procedure was performed for each ligand with the force field (Force Field) of mmff94 (Merck Molecular Force Field) and the 'steepest descent' optimization algorithm for all ligands. The ligands were charged with gasteiger using the ADT program, the 'detect root' command was determined as root, and all bindings except amide bonds were set freely. The grid was calculated to cover all the atoms responsible for catalytic activity in the protein by taking the catalytic region of the target as the center and not to restrict the rotation angles of the ligand at this center. In the MD simulation, after 9 ns, protein snapshots at every 100 ps were docked with their own substrates for each polymorphic structure to test the experimental validation. In this validation study, Pearson correlation coefficient values were obtained by comparing docking scores with experimental enzyme activities in the literature (Draganov et al., 2005). Each polymorphic structure showing the highest correlation the experimental data of Draganov et al. (2005) was selected to calculate the affinity of Covid-19 drugs. The snapshot with the highest correlation for each polymorphic structure was selected for docking calculations. A total of 31 compounds were processed, 14 of which were PON1 substrates and 17 were Covid-19 drugs. The flowchart of the study was demonstrated in Figure 1.

Free energy calculations by MM/PBSA method
Free energy calculations were made according to MM/PBSA approach for the 3 compounds with the highest docking score for each polymorphic structure and wild type protein structure (Wang et al., 2019). MM/PBSA calculations were carried out with the g_mmpbsa tool that integrates with the GROMACS simulation package (Kumari et al., 2014). The binding free energies were calculated by making 100 measurements in 100 ps intervals within the last 10 ns of the 100 ns MD simulation. Partial charges were taken from force field parameter in Poisson-Boltzmann (PB) calculations. The dissolved dielectric constant was set as 2 in vacuum electrostatic calculation, and the solvent dielectric constant was set as 80 in vacuum electrostatic calculation. The Solvent-Accessible Surface Area (SASA) was used to estimate the non-polar contribution. The entropy change was estimated using the Schlitter formula (Schlitter, 1993).
Binding free energy is calculated based on the following equations (Duzgun & Eroglu, 2020):

Results and discussion
Homology modeling 20 models for human PON1 polymorphic structures (PON1-w, PON1L55M, PON1Q192R) were obtained, by using 3SRE protein structures as the templates and by homology modeling using the MODELLER Software (Ben-David et al., 2012;Webb & Sali, 2014). The values of the build quality statistics of these proteins were given in Table 2. The best structures were chosen according to the 'MolProbity Score' which is a value obtained by the combination of atomic collisions, rotation angles and Ramachandran evaluations normalized with the resolution of the X-RAY structure (Chen et al., 2010). There is a negative relationship between the score and the quality of the building that a score value of 1 represents the best but a value of 100 the worst (Chen et al., 2010). Although Ramachandran favored ! 98% is considered ideal, the values > 92% are acceptable. Below 1% is considered ideal for poor rotamers (Hu et al., 2009). 3SRE coded rabbit-human hybrid rePON1-G2E6 structure (approximately 84% amino acid compatibility), and 100% human PON1, PON1L55M and PON1Q192R polymorphic structures obtained by homology modeling from PDB.org site were aligned by using Chimera software (Ben-David et al., 2012;Pettersen et al., 2004). Since the L1 node region is extremely flexible, the greatest deviation in alignment occurred in this region. In other regions, pretty much perfect harmony was observed (Figure 2). Hu et al. (2009) found that the L1 and L2 regions in the PON1 enzyme are more active than the other regions and play a critical role in enzymatic activity by blocking the catalytic region.
RMSD, RMSF and Rg in pre-docking molecular dynamics simulation RMSD is a measure of the stability of the system and the fitness of the protein structure. The RMSD value is usually found by calculating the shift relative to the position in the alpha carbon structures of the protein (Maiorov & Crippen, 1994). The system has become stable if the RMSD value reaches its equilibrium and does not show much drift. The fluctuations below 0.1 nm showed that the system was in equilibrium and the proteins were stable. The 10 ns RMSD graph is shown in Figure 3.
The most characteristic feature of a protein is that it has a certain flexibility. While this flexibility is observed with MD simulations, the RMSF curve is used to graphically show the plasticity per residue. The RMSF curve expresses the mobility of each residue on the protein within the structure during the MD simulations (Maiorov & Crippen, 1994). Figure 4 shows the RMSF value of each residue. Regions denoted by L1, L2 and L3 represent knot regions on protein structures, while H1 and H2 represent helical structures. The most flexible region is H1 at the N-terminal end of protein structures, while the L1 region represents the second most flexible region. The L1 region is also close to the active site (Hu et al., 2009). In a hybrid in vivo and in silico study in which Maiorov and Crippen (1994) investigated the effect of titanium dioxide nanoparticle on serum PON1 activity in rats, the PON1 structure showed flexible L1 (residues 72-79) and L2 (residues 290-300) loop regions with MD simulation. This result is consistent with our findings. In addition, in our predocking MD simulation study, we calculated the L1 loop region as the most flexible region in wild type. While the decreased flexibility in the PON1Q192R structure was observed, it was almost rigid in the PON1L55M polymorphic structure (Figure 4). The fact that the L1 region is close to the catalytic center and that it acts as a gate has shown us that the polymorphisms affect the catalytic mechanism of the enzyme (Hu et al., 2009;Mortazavi et al., 2021). No significant change was observed in other regions. The reason why the fluctuation especially in the L1 region is higher in PON1w compared to other mutant structures may be the steric effect of mutations in other structures that limit flexibility in the protein structure. Mota et al. (2019) showed that the PON1-w type had higher enzyme activity than the PON1L55M mutant.
B-factor refers to the displacement of atoms within a protein structure. If there are flexible regions in the protein structure, the displacement range of the atoms will be greater in this region (Kuzmanic & Zagrovic, 2010). Protein structures can be colored according to the B-factor by the graphics programs such as Chimera (Pettersen et al., 2004). Figure 5 shows the flexible regions of PON1 protein and its polymorphic structures in terms of B-factor. The blue and red   The radius of gyration (Rg) shows that whether the structure has reached a stable conformation in MD simulations. If the protein is well folded, the Rg value remains at a certain constant value. The greater the deviation indicates that the protein changes conformation or begins to deteriorate (Lobanov et al., 2008). Rg values for each polymorphic structure in our study are shown in Figure 6. As can be seen in this figure, all three protein structures are conformationally stable.

Molecular docking
In Figures 1-3, it has been shown in which conformation that all substrates interact with PON1-w, PON1L55M and PON1Q192R polymorphic structures, respectively. In PON1-w and PON1L55M polymorphic structures, Y71 residue was observed to form open conformation for each substrate, while Q192R showed closed conformation in polymorphic structure. Looking at Figure 2, it was seen that the closed form of Y71 had a positive effect on the affinity of the substrates. F292 residue was shown on the graphs as it was left  flexible in the molecular docking process. It was observed that it remained outside the active center and did not interact with the substrates (Hu et al., 2009). The detailed analysis of the active region in Apo form, showed that Y71 displayed open and closed conformations by acting as a gate in front of the active site. Conformational transitions at residue Y71 below 10 ns show that the free energy barrier between the two transitions is very low (Hu et al., 2009). Our molecular docking scores were compared with the experimental PON1substrate inhibition values obtained by Draganov et al. (2005) according to the Pearson correlation coefficient. For each protein, snapshots were taken at 100 ps intervals in the last 1 ns and when their interactions with PON1 substrates were examined, PON1-w at 9.8 ns, PON1L55M at 9.4 ns and PON1Q192R variants at 9.2 ns, with Pearson correlation coefficient values À0.94, À0.76 and À0.84, respectively ( Figure  7). With this method, which can be classified as ensemble docking, we investigated the snapshot of the conformational changes in the protein active center during the MD simulations, which has the best correlation with the experimental data. Ensemble docking has emerged as a frequently used method in early drug identifications in the pharmaceutical industries and the literatures (Amaro et al., 2018;Evangelista Falcon et al., 2019). In our study, a high Pearson correlation coefficient of over 0.9 was observed between PON1-w structure and ligand interactions, while this value for polymorphic      variants was in the range of 0.75-0.85. The reason for this was thought to be that the highly flexible L1 region in MD simulation of the polymorphic forms was not as mobile as the wild type (Figures 4 and 5).
After MD simulations, the interactions of the mentioned Covid-19 drugs on the snapshots of the PON1 proteins, which showed the best correlations with the experimental data, were calculated by molecular docking. As seen in Figure 8, the compound with the best score in PON1-w structure was ruxolitinib, while dexamethasone in the M/L55 variant and galidesivir in the Q192R variant. It can be said that the drugs such as ruxolitinib, dexamethasone, sitagliptin, galidesivir, which have the highest scores in interaction with proteins, may be potent candidates for inhibiting PON1 enzyme activity.
The positive scores of the drugs with relatively high molecular weight compared to PON1 ligands in the docking calculation showed us that these compounds are thermodynamically unable to interact with PON1. This means that these drugs can be safe drugs as they cannot inhibit PON1 enzyme. These drugs were ivermectin, ritonavir and azithromycin ( Figure 8). There are some studies on that shows the effects of these drugs on PON1 activity. Supporting to our results, Maselli et al. (2014) found that lopinavir/ritonavir did not significantly affect PON1 activity in HIV-1-infected patients. On the other hand, Pastryk, Rusek & Bełtowski (2016) investigated the effects of ritonavir on PON1 activity in rats and showed that ritonavir reduced PON1 activity. Another study reported that azithromycin inhibits PON1 enzyme competitively by interacting with the amino acids in enzyme active site (Akin & Dilek, 2018). According to the study of Can et al. (2007), the treatment of children with bronchial asthma with montelukast, which had the low affinity in our study, increased PON1 activity. The other study with montelukast found mixed inhibition of PON1 activity (G€ okc¸e et al., 2019).
When the interactions of drugs within PON1-w structure were examined, ruxolitinib, dexamethasone and colchicine showed the best molecular docking scores. Ruxolitinib formed hydrogen bonds with Asn224 and Pro72 while exhibited pi-alkyl interaction at Leu240, Leu267, Ile291 and Val346. There is no study, in literature, related to PON1-ruxolitinib interaction. Dexamethasone interacted toAsp269 residue in the active site of PON1-w by forming hydrogen bond (Figures 9 and S4). Various in vitro studies showed that dexamethasone significantly decreased PON1 activity but increased PON1 mRNA expression explained that this drug noncompetitively inhibitedPON1 (Arslan et al., 2012;Beydemir et al., n.d.;Sinan et al., 2006). On the other hand, Hu et al. (2009) reported that Asn224 and Asn168 were important residues in the catalytic site whereas Pro72 residue had a weaker role. Besides, it was shown that Asp269 residue in the paraoxon played an important role for increasing PON1 inhibition activity. When the interactions between the colchicine and PON1 were examined, it was observed that hydrogen bonds were formed with Asn224 and Asn168, pi-alkylationswithVal346, Ile74 and Ile291 and pi-sigma bond with His285. Demir et al. (2014) explained that the colchicine inhibited PON1 activity semi-competively in patients with Behc¸et's Disease while it activated in healthy people.
According to molecular docking studies on the PON1L55M enzyme, it was observed that dexamethasone, sitagliptin and baricitinib drugs had the best scores. Sitagliptin formed hydrogen bond with Glu53 and showed pi-alkyl bond with Leu28, Leu267 and Ile291 residues. In addition, halogen interaction was observed between Leu69 and this drug. Trocha et al. (2019) reported that sitagliptin increased PON1 arylesterase activity in rat livers subjected to ischemia and reperfusion. On the other hand, while dexamethasone formed hydrogen bonds with Pro72 and Asn224 residues, baricitinib formed it withHis115, Asn168, Asn224 and Phe292.Additionally,it was observed that baricitinib also formed a pi-sigma bond with Ile291 (Figures 10 and S5). No study has been found regarding the activity of baricitinib on PON1.
Galidesivir, ruxolitinib and hydroxychloroquine showed the best molecular docking scores for PON1Q192R ( Figures  11 and S6). Galidesivir interacted with the hydrogen bonds at the residues of Glu53, Pro72, Asn168, Asn224, Asp269 and Phe292, and exhibited pi-alkyl bonds at Ile291. Ruxolitinib made hydrogen bonds with Glu53, Asp269 and His285. Apart from that, it was observed that the compound formed a large number of pi-alkyl bonds. Hydroxychloroquine formed hydrogen bonds with Leu69 and Pro72 residues and made pi-alkyl bonds with Tyr197, Phe222, Leu241 and Ile291.Uzar et al.found that hydroxychloroquine decreased PON1 activity muscle tissue of rats and they concluded that this drug induces oxidative stress (Uzar et al., 2012). An increase in PON1 activity was observed in a study investigating the combined therapy of methotrexate, sulfasalazine and hydroxychloroquine on patients (Charles-Schoeman et al., 2017).
When the RMSD, RMSF and GYRATE graphs were examined as a result of the MD simulation of the interactions of PON1-w and other variants of PON1L55M and PON1Q192R with selected compounds, it was observed that the proteincompound complexes were in equilibrium and stable. It was observed that fluctuations in the RMSD graph were at a maximum level of 0.1 nm and the Rg interval in the GYRATE graph did not exceed 0.5 nm. On the RMSF graph, residues between 72-80 and 294-298 belonging to the loop region were observed to be mobile (Figure 12).
When the binding free energies of protein-compound interactions as a result of molecular dynamics simulation were calculated according to the MM/PBSA method, it was observed that colchicine was the most effective compound for PON1-w (-123.657 kJ/mol). Interestingly, Sitagliptin (-509.409 kJ/mol) and followed dexamethasone (-110.622 kJ/ mol) showed very high affinities for PON1L55M. On the other hand, galidesivir and hydroxychloroquine exhibited strong affinity (-389.753 and À360.998 kJ/mol, respectively) for PON1Q192R. 2-Coumaranone, the natural substrate for PON1, showed affinity less than 100 kJ/mol, for all three forms of proteins ( Figure 13 and Table 3). These results showed us that selected compounds can bind more effectively in all polymorphic structures than 2-coumaranone. The fact that Sitagliptin has very low binding free energy in PON1L55M form made us think that the risk of atherosclerosis may arise with PON1 inhibition in a treatment that can be done with this drug. The fact that dexamethasone and ruxolitinib have similar binding free energies in different variants, evokes that sitagliptin may have very low values in other variants.

Conclusion
In this study, we identified the interactions of the potential Covid-19 drugs (azithromycin, baricitinib, chloroquine, colchicine, dexamethasone, favipiravir, galidesivir, hydroxychloroquine, ivermetcin, lopinavir, montelukast, nitazoxanide, remdesivir, ribavirin, ritonavir, ruxolitinib, sitagliptin) with wild PON1-w and polymorphic forms PON1Q192R and PON1L55M by using molecular docking method. It was found that ruxolitinib, dexamethasone, colchicine; dexamethasone, sitagliptin, baricitinib; galidesivir, ruxolitinib, hydroxychloroquine were shown higher affinity for PON1-w, PON1L55M, PON1Q192R, respectively. The lowest docking scores for three PON1 structures were obtained with ivermectin, ritonavir and azithromycin (Figure 8). We suggested that considering the protective effect of PON1 enzyme on cardiovascular diseases and its antiviral effect on SARS-CoV-2 which causes to Covid-19, attention should be paid in the use of the drugs that bind PON1 enzyme with high affinity. Although the study is validated with an experimental study during the docking phase, further in vitro and in vivo studies are required to advance it.

Acknowledgments
The numerical calculations reported in this paper were partially performed at TUBITAKULAKBIM, High Performance, and Grid Computing Center (TRUBA resources).

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
The author(s) reported there is no funding associated with the work featured in this article. Authors' contributions ZD is responsible for the preparation and implementation of the methodology and writing the manuscript. BVK and AO are responsible for designing the concept and reviewing the manuscript. IK is responsible for the analysis of data and critical evaluation of the manuscript.

Availability of data
All data are available with reasonable request.

Code availability
N/A