A Computational Study of Ivermectin and Doxycycline Combination Drug Against SARS-CoV-2 Infection

In the present study, we have described how by using molecular docking and molecular dynamics (MD) simulation studies the combination drug of ivermectin and doxycycline can be used as a potential inhibitor for Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV) virus. In lieu of unavailability of specific cure of coronavirus disease of 2019 (COVID-19) till now various possibilities for individual and combination drugs have been explored by the medical practitioners/scientists for the remedial purpose of CoV-2 infections. 3C-like protease (3CL pro ) is the main protease of SARS-CoV-2 virus which plays an essential role in mediating viral replication in the human body. 3CL pro protein can serve as an attractive drug target. In this work, we have studied drug: 3CL pro interactions by in-silico molecular docking and MD simulation approaches. Common and easily available antiviral drugs ivermectin, doxycycline and their combination can regulate 3CL pro protein's function due to its easy inhibition.


Introduction:
In the year 2020, the COVID-19 disease has spread globally and it has become an ongoing pandemic. Reported by the World Health Organization (WHO), due to this pandemic disease, more than 35,659,007 numbers of active patients with 1,044,269 people have already died till 10 October 2020 (https://covid19.who.int/). WHO declared the COVID-19 as a global health emergency. This disease is caused by a member of the coronavirus family [1]. Coronavirus was first found in 1930 in domestic poultry [2]. After that they were identified as causing several diseases in humans such as; respiratory illness, neurological, liver diseases, etc. [3]. Till now seven categories of this virus were identified.
Among the seven categories of coronavirus, four causes only common cold with mild symptoms and in very rare cases pneumonia, respiratory infections in infants and older people [4]. The other three categories are Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV) [5], Middle East Respiratory Syndrome Coronavirus (MERS-CoV) [6] and lastly the new one known as SARS-CoV-2 [7] identified in 2003, 2012 and 2019 respectively. The international committee on taxonomy of viruses declared this new novel coronavirus as SARS-CoV-2 [8]. The SARS-CoV-2 is a single-stranded RNA virus and belongs to the Coronaviridae family having genome sequences of 79.5% sequence matching [9,10]. This shows that bats may be the carrier of this virus. The uniqueness of this virus is the presence of spike glycoproteins on its surface which gives a crown-like appearance of the virus structure. The crown-like spike protein surface of this virus can be easily visible with the help of electron microscopes.
These spike proteins are a very significant part of SARS-CoV-2 [11] virus as they can easily interact with the human proteins which coats the inside of the nose and the cells of lungs. The interaction of spike protein and human protein causes change in spike protein of CoV-2 shape and causes the human receptor cell to swallow up the virus. Through the receptor binding domain (RBD), glycoproteins of the viruses start binding and entering to the host cells. The key receptor for SARS-CoV-2 in humans is angiotensin converting enzyme 2 (ACE2) [12]. After entering the host cell, different human protease like airway trypsin-like protease (HAT), cathepsins and trans membrane protease serine 2 (TMPRSS2) divide the glycoproteins of the virus and so the conformational alteration of the virus structure occurs.
From this phase the transformed virus replicates itself very fastly through some cyclic processes [12] and starts infecting the neighboring cells like lung, heart, brain cells and many others. From studies, scientists showed that the spike glycoproteins of coronavirus attach on the cell surface of the ACE2 receptor in the human body and allows the virus's genetic material to enter the human cell [13]. Virus's genetic material proceeds to hijack the metabolism of the cell and help the virus to divide.
To overcome this disease the whole world is in a race to find vaccines/drugs to attack this virus.
Through clinical trials around 200 drugs and vaccines (approved by Food and Drug Administration).
Covaxin, INO-4800, mRNA-1273, NVX-CoV2373, BBV152 etc. are some candidate vaccines that are currently under trials for COVID -19 [14]. In the eleven months since the SARS-CoV-2 virus werediscovered, the scientific community has put forward an extraordinary effort that has resulted in the creation vaccines. Similarly examples of some FDA approved drugs for COVID-19 are atazanavir, remdesivir, ritonavir, lopinavir, chloroquine, hydroxychloroquine (HCQ), cyclosporin, favipiravir etc. [15 -19 ]. Screening of potential drug from different medicinal plants extracts for SARS-CoV-2 is also going on [20][21][22]. Now according to most common treatment protocols since there is no detected and approved drug for COVID-19, patients with severe COVID-19 symptoms are usually treated by different purposed antiviral drugs as trial basis. Most of the above-mentioned drugs are usually antiviral in nature and are used for various viral diseases like: HIV medication, influenza, MERS and SARS diseases or for enhancing the immune system of human life [23][24][25][26]. Nowadays to identify potential drugs for various diseases, the concept of drug repurposing is widely used. Drug repurposing is an approach to find out the new uses for already available drugs that are originally developed for specific diseases [27]. Drug repurposing process has already proved to be very effective since many drugs have multiple protein targets and genetic factors; molecular pathways which can be shared by diverse diseases. For many years repurposing of drugs have been used such as favipiravir drug used for influenza virus, sofosbuvir drug used for hepatitis C virus have a strong repurposing prospective against Zika and Ebola [28], drugs oseltamivir, lopinavir, nelfinavir, atazanavir and ritonavir have been used for the treatment SARS and MERS [29,30]. But these drugs have their own toxicity related issues. On the other hand, some immunomodulatory plasma-based therapies are in use. Some food nutrients, herbal medicines having antiviral and immunity building properties are considered as an alternative of COVID-19 therapies [31,32]. In the same way, a repurposing of combination drugs with ribavirin, lopinavir, and ritonavir have already been anticipated for the COVID-19 patients [40]. Lopinavir and ritonavir combination are already in use for HIV treatment. However, the efficacy of the vaccines are almost 70-80%. So there is an urgent and strong requirement for a newly invented drug/repurposed drug/combination drug to fight the disease.
A combination drug includes two or more than two active ingredients mixed in a single dose form. For many years combination of drugs has been used for treating diseases such as aspirin/paracetamol/caffeine combination (Excedrin) is used for the treatment of headache and migraine [33], Carbidopa/levodopa/entacapone is used for the treatment of Parkinson's disease [34], and indacaterol/mometasone, used for the treatment of asthma [35]. Combination drug therapy is applied for many diseases such as: tuberculosis, leprosy, cancer, bacterial infections, malaria, and for many viral diseases like influenza, HIV/AIDS etc. [36]. Recently two combination drugs of Nitazoxanide/azithromycin [37] and another combination drug: lopinavir/oseltamivir/ritonavir are [38] being largely in use by medical practitioners to fight against SARS-CoV-2 infections. There are several advantages to the combination of drugs. They are increased action of drugs and efficiency, increase the efficiency of the therapeutic effect, reduced cost and side effects. However, combinations of drugs also include some disadvantages. Dose must be given in some fixed ratio otherwise mismatched pharmacokinetics may increase severe toxicity effects. Though several clinical trials are underway to identify drugs against SARS-CoV-2, but still currently there is availability of single approved drugs or vaccines. Urgent requirement of cure of current medical emergencies due to COVID-19 motivated us to investigate the possibility of inhibition of SARS-CoV-2 by using some repurposing of combination drugs: ivermectin and doxycycline.
In the present paper, we have described how the combination drug of ivermectin and doxycycline, can be used as a potential SARS-CoV-2 3CL pro inhibitor. 3CLpro is synonymous to another name which is Mpro (main protease). For several years ivermectin (C48H72O14) is used to treat many infectious diseases in mammals [39].
Ivermectin is a antiparasitic agent and doxycycline is a broad-spectrum tetracycline antibiotic.
Ivermectin have also antiviral activity against both RNA and DNA viruses [40]. Recently in April 2020, the in-vitro activity of Ivermectin against SARS-CoV-2 was reported [41]. Ivermectin's antiviral mechanism of action in COVID-19 may be block the activity of α/β1 receptors, which inhibiting viral protein transport in and out of the host nucleus [42]. Doxycycline was shown to have anti-SARS-CoV-2 activity in contaminated Vero E6 cells in vitro [43]. It's antiviral activity may be mediated by upregulation of the zinc finger antiviral protein, which binds to viral messenger RNAs and inhibits viral RNA translation [44]. Furthermore, doxycycline's anti-inflammatory effects were thought to add to its effectiveness in pulmonary inflammatory conditions such as asthma, cystic fibrosis, bronchiectasis, and lung damage. Both of these drugs are low cost and safe [45]. The aim of the present study that the combination of ivermectin and doxycycline can be proved to be effective against SARS-CoV-2 so that medical professionals may have alternative tool to treat patients. We have performed molecular docking and molecular dynamics (MD) simulations to understand the interaction mechanism of the proposed drugs for COVID-19.We hope that this work will provide other researchers with an important investigation way to identify new COVID-19 treatment.

Protein structure preparation
Coronavirus possesses a number of polyproteins (structural and nonstructural). Among them 3CL pro is a key CoV enzyme which plays an important role in mediating viral replication and transcription with the help of its glycoprotein. To rapidly discover the targeted drugs for clinical use, researchers focused on identifying drug leads that target 3CL pro protein of SARS-CoV-2 as it plays an important role for viral replication and transcription. In the present work, we have used one of 3CL pro proteases of CoV-2 virus in a complex with an inhibitor N3 (PDBID: 6LU7) [46,47] as the target protein. 6LU7 has been shown to be a promising target for designing COVID-19 drugs. We have chosen 6LU7 for checking the inhibiting and binding properties of it with the ivermectin and doxycycline drugs.
The structure of SARS-CoV-2 protease (PDB ID: 6LU7) was used as a receptor and retrieved from Protein Data Bank (http://www.rcsb.org/) [48,49] and are shown in Figure 1 (a). We have removed water and hydrogen from it. All the existing properties of the drugs are described in Table 1. For the preparation of protein, we have used Auto Dock and MG Tools of AutoDockVina software [50]. At first existing lead components, water molecules and ions have been removed from it. Later the process of cleaning has been done. We have calculated the Gasteiger charges of protein structures and after that polar hydrogen have been introduced. Then the non-polar bonds were merged and rotatable bonds were defined. Finally, by using Discovery studio 2020 [51] the intrinsic ligands were detached from the protein molecules and the final protein molecule was saved in the PDB format (Figure 2 a). For target protein by visualizing the dihedral angles ψ against φ of amino acid residues, Ramachandran plots have been drawn (Figure 2b). It predicted permissible and disfavored values of ψ and φ. Figure 2b shows Ramachandran plots for 6LU7; the plot specifies localization on chain residues, which reflect the consistency of the protein structure, implying effective and accurate docking capacity.

Figure 2. a) Target variable viral proteins (6LU7) SARS-CoV-2 protease enzyme as receptor and b)
Ramachandran plot for the receptor protein.

Ligand drug molecules preparations
Structures of the drug molecules were downloaded from Drug Bank in pdb format. Then these structure were fully optimized by using the Gaussian 09 program [52]. We have used the optimized structure for docking analysis as they provide better results than unoptimized one. The geometric optimization of all drug compounds were carried out using HartreeFock (HF) and STO 3G basis set.Gauss View 5 molecular visualization program was used for visualizing the optimized structure [53]. ADME-T properties of molecules were identified using Organic chemistry portal (http://www.organicchemistry.org/prog), a web based application for predicting in-silicoADME-T property. Protein-ligand interactive visualization and analysis was carried out in AutoDock 4.2 software on Windows 7 (64-bit).
For the present work, we have selected two potential ligand drugs: ivermectin (C48H72O14) and doxycycline (C22H24N2O8). Detail structures of these molecules were downloaded from Drug Bank in pdb format ( Figure 1 and Table 1). Different chemical, physical, drug likeness and pharmacokinetics properties obtained from SWISS ADME are shown in Table 1. Both the proposed drug molecules have molecular weight less than 875 gm/mol and topological polar surface area (TPSA) values less than 180 Å. 2 ( Table 1). All drug molecules have H-bond donors ≥6, H-bond acceptor ≥14 and have low synthetic accessibility count, this suggests that they can be synthesized easily. Though these drugs violate some drug likeness properties, still the availability of these drugs in the drug industry motivates us to consider these as potential inhibitors. The ligand file in pdbqt format is needed for molecular docking study with AutoDock Tools. AutoDock Tools 1.5.6 [54] have been used to save ligands in pdbqt format.

Methods: Molecular docking and Molecular dynamics simulations
To predict the target and drug interactions, molecular docking is commonly used in simulation. It minimizes the energy and calculates the binding energy of the interactions. In the molecular docking simulation, we normally make out the best pose of the ligand towards the receptor protein with the help of scoring functions [55]. Molecular docking can show the possibility of any biochemical reaction or whether a drug is docked with the receptor protein or not. The AutoDockvina with the best fitted parameters binding modes: 9, exhaustiveness: 8, applied maximum energy difference: 3 kcal/mol and Grid box center with x, y, and z coordinate of residue position of the protein is used for docking purpose [50]. Grid box was formed with centers of x, y, z coordinate of residue position of the receptor protein respectively. The value of the centers of x, y and z coordinates were considered as -10.729204 Å,  [51] have been used. After the analysis of individual docking, sequential docking is performed. For sequential docking, the grid box coordinates were set to the particular binding region of each drug with default grid spacing. In the procedure of sequential docking, the first ligand is docked and the complex is saved out as a single file, where the first ligand is considered part of the receptor.
Docking is then carried out on this complex with the second ligand. The structural dynamics of receptor and inhibitor interaction and thermodynamics stability of ligand: protein have been investigated with the help of Linux based platform "GROMACS 5.1 Package'' [56], Different thermodynamic parameters like temperature (T), density (D), potential energy (Epot), root mean square deviation (RMSD) for backbone, root mean square fluctuation (RMSF) for protein CαSolvent accessible surface area (SASA), intermolecular hydrogen bonds, interaction energies (∆G) of the protein and drug complex have been find out with GROMOS43A2 force fields [57]. For topology creation, we have used "PRODRG" server. PRODRG works with the concept of charge groups, which are defined as a group of bonded atoms with an integer charge. To assign atomic charges it recognizes the charge groups first. After topology creation of protein and ligand, according to the procedure followed for MD simulation, aqueous solution simulations have been performed using the water model: TIP3P. For solvation process protein in apo state, protein:ligand complexes were solvated in a cubic box, with a buffer distance of 10Å and volume as 893,000A 3 . For electrically neutralizing the system four Na + ions have been added. Then we minimize the energy in the vacuum. For energy minimization 50000 iterations have been taken. To check the stability of the system, we have been performed MD simulation for the period of 0 ps to 100000 ps. Number of particles (N), volume (V), and temperature (T) were constant under the 1 atmosphere pressure and 298K temperature. We have used Lennard-Jones and Coulomb short range interaction for the nonbonded interactions. All simulations were performed using a Berendsen thermostat and barostat [58] with the coupling time of 0.1 ps and 0.5 ps, respectively. Non-bonded interactions (electrostatic and LJ interaction) were calculated using a triple range scheme within a shorter range cutoff of 0.8 nm. Graphical tool Origin pro has been used to study the simulated results.
Where, ∆ is the molecular mechanical energy changes in gas phase and is the sum of covalent ∆ , electrostatic (∆ ), and van der Waals energy (∆ ) changes.
Covalent energy is the combination of bond angle and torsion and ∆ , is separated into its polar and nonpolar contributions.∆ , is solvation free energy change and -TΔS conformational energy change due to binding.For RMSD and RMSF multiple simulations were performed independently to validate the results obtained. MD simulation can simulate in ps/ns or further finer temporal stead-fastness [60]. The MD simulation force field plays an important role for estimating the forces within the molecule (intramolecular force) and between two molecules (intermolecular force). These intermolecular and intramolecular forces used to calculate the potential energy of the molecules. The total energy of the system is given as the sum of bonded and non-bonded energy and given as below: These equations show that the bonded energy is the combination of bond, angle and dihedral energies while nonbonded energy is the combination of hydrogen bond, electrostatic and van der waals energies (eq. 6, 7).

Computational facility
MD simulations and corresponding energy calculations have been computed using HP Intel Core i5 -1035G1 CPU and 8 GB of RAM with Intel UHD Graphics and a 512 GB SSD.

Individual docking of drugs against SARS-CoV-2 protease
In the present work, ivermectin and doxycycline drugs were docked to SARS-CoV-2 main protease (3CL pro ). Ivermectin and doxycycline drugs confirm the RO5,which means Lipinski's rule of fiveand other drug likeness rules etc. Hence, we have shown their strong application as potential drugs reaching the market (Table 1).
For the best docked 3 pose of ivermectin: protein complex, the obtained value of ki as 8.7 X 10 -6 M which proves the strong attraction of ivermectin towards protein ( Table 2)

Sequential docking of two drugs against SARS-CoV-2 protease
Individually, ivermectin and doxycycline drugs showed a good binding energy of -6.9 kcal/mol and -6.4 kcal/mol, respectively. The docked ligand molecules with the protease 3CL pro (6LU7) are shown in Figure 4a,b. The possible hydrophobic interactions and hydrogen bond between 3CL pro of the two considered drugs obtained with individual docking are presented in Table 1. We have performed sequential docking for checking the interaction of combinational drugs (two or more than two drugs  (Table 3).

MD simulation analysis
To analyze the stability of the studied structure, MD simulation of the complexes (ivermectin: 3CL pro , doxycycline: 3CL pro , ivermectin+doxycycline: 3CL pro ) have been studied for the period of 0ps to 100000 ps. For MD simulation, first we have to make all the structure energetically optimized (the potential energy should be minimum and negative with a maximum force value). Figure 5 represents energetically minimized protein and complex systems. We have obtained steady convergence of potential energy for all the cases. The comparison of the potential energy (Epot) of the stable structure of apo 6LU7 and in drugs: 3CL pro complex has been done carefully. In the apo state 3CL pro has Epot of -1.27x10 6 ±56.7 kJ/mol, while the complex ivermectin: 3CL pro , doxycycline: 3CL pro and ivermectin+doxycycline: 3CL pro has an average Epot of -0.6110×10 5 ±2.62 Kcal/mol, -0.60 × 10 5 ±3.34 Kcal/mol, -0.594×10 5 ±12.66 Kcal/mol respectively (Table 3). Now all the structures having their lowest Epot values are ready for MD simulation.  To stabilized different parameters (temperature (T), pressure (P), density (D), volume (V) etc.) within a time scale of 100 ps to 10000 ps, we have further check the optimized drugs: 3CL pro structures equilibrated by NVT and NPT ensembles. It is observed that over the period of 100 ps time trajectory the temperature of the complex rapidly reached the stable at 300 K (room temperature) value ( Figure 6).
This temperature stability is maintained throughout the process. The temperature, pressure and density values of the system were also observed to be very stable over the period of time trajectory 100 ps (SD 3, SD 4). This concludes that the system is well equilibrated and prepared for MD simulation.
The compactness of the system with respect to time of apo protein and protein: ligand complex can be measured with the help of radius of gyration (Rg) [62]. Normally for the stably folded protein structures the values of Rg keeps a relatively steady for full time scale [63]. Whereas the Rg values for the unfolded protein keeps changing for full time scale. Less compactness in the structures and high compactness with more stability exhibit a low and high Rg value respectively. In the present paper we have observed the apo protein (3CL pro ) has an average value ofRg as 2.225 nm (SD 5, Table 3). Almost similar variation is observed with the proposed drugs: 3CL pro complex (SD 5). This shows high compactness with more stability in the protein and drugs complexes (SD 5).
Further to validate the applicability of ivermectin, doxycycline and ivermectin+doxycycline ligands as proposed drug for COVID-19, we have simulated the SASA. SASA measures the area of exposure of the receptor to the solvents. The higher value of SASA indicates that the drug is more inserted into the water whereas, lower value represents that more drug is covered by the protein, which represents better complexation.In the present work, we have obtained the SASA mean value 22 nm 2 for apo protein (SD 6, Table 3). Similarly, for all the proposed drug and 3CL pro complex the mean value of SASA is 9 nm 2 . The low computed values of SASA observed for all drugs: protein complex shows that drug binding with the receptor protein increases the exposure of complexes to the protein (SD 6). Which validates the best complexation possibility. Intermolecular hydrogen bonding plays a significant role to get an idea about the binding strength between protein and drug. Ivermectin has a stable range of intermolecular hydrogen bonding with protein between 0 to 7 with an average value 3.5 in throughout the whole simulation process (Figure 7, Table 3). Doxycycline has a range of intermolecular hydrogen bonding with protein between 0 to 6 in throughout the whole simulation process with an average value of 3. However, the combination of both the drugs (ivermectin+doxycycline) has the highest stable range of intermolecular hydrogen bonding with protein between 0 to 12 with an average value 7 ( Table 3)  RMSD corresponds to any change in the conformational stability of the protein: drug complex and in the protein dynamics. RMSD of the free protein and protein: ligand complex have been simulated to 100000 psby using MD simulations. RMSD and RMSF have been measured by using the GROMACS module at an interval of 1000 ps. RMSD variation of apo 3CL pro lies in the range from 0.08 to 0.16Å. Ivermectin: 3CL pro , doxycycline: 3CL pro , ivermectin+doxycycline: 3CL pro complex, also ranges RMSD values from 0.08 to 0.16 Å ( Table 3). The RMSD value for complexes exactly matches with the apo protein. This provided a suitable basis for our study by the better stability with the probe drugs. Figure   For all amino acid residues with respect to Cα atom RMSF have been simulated. RMSF plot for 3CL pro in its apo state and ivermectin+doxycycline: 3CL pro complex have been shown in Figure 9, which depicts the fluctuations at the residue level. Residue fluctuation profile for both the cases shows a similar trend having an average RMSF value of 0.15 Å, which indicates that binding of both the drugs to the 3CL pro had no key effect on the flexibility of the protein and was quite stable. Jones interaction energy provides the total interaction energy. Figure 10 a (Table 3). Table 3 represents all the Coulombic interaction energy and Lennard Jones interaction energy for individual drugs: protein complex and combination of drugs: protein complex. The comparison suggests that for all the complex formation, short-range Lennard-Jones has shown stronger effect on binding affinity than the short range coulombic interaction energy.  [64][65][66]. Figure 11 represents the ΔG values for ivermectin: 3CL pro , doxycycline: 3CL pro and ivermectin+doxycycline: 3CL pro complex with respect to the time trajectory 0 ps to 10000 ps and inset of Figure 11 represents the ΔG values for ivermectin+doxycycline: Kcal/mol) ( Table 3, SD 7, SD 8, SD 9). The binding energy graph is going up (positive energy) for ivermectin+doxycycline: 3CL pro complex after 8000 ps (Inset of Figure 11). However, it is going down (negative energy) for ivermectin+doxycycline: 3CL pro complex after 10000 ps (Inset of Figure 11). For

Conclusions
In conclusion, two drugs (ivermectin and doxycycline) were tested as potential inhibitors for COVID-19 main protease 3CL pro via molecular docking. A strong inhibitory possibility of proposed drugs for SARS-CoV-2 protease 3CL pro was verified by Gastrointestinal absorption, pharmacokinetics, drug likeness, and medicinal chemistry properties by using ADME analysis. From docked compounds, we have proposed that ivermectin and doxycycline demonstrated high binding affinity to the 3CL pro and their combined docking increases the binding affinity on COVID-19 main protease. Strong binding affinity, lowest inhibition constant and existence of hydrogen bonded interaction established the better stability of ivermectin+doxycycline: 3CL pro complex structure. Further studies also conducted on these compounds using MD simulations in order to get more reliable data. Many thermodynamic parameters (Epot, T, V, D, Rg, SASA energy) obtained by MD simulation also validated the complexation between ivermectin+doxycycline and 3CL pro . The backbone of the complex and free 3CL pro illustrate similar RMSD and RMSF, which demonstrate the stability of the binding of drugs and protein. MD analyses have also confirmed the complexation between proposed drug and 3CL pro by the lower values of binding energy. All simulated results establish that combination of drugs is a stronger candidate as a potential inhibitor for SARS-CoV-2 than considering each drug separately. Our present in-silico study would provide a new approach to the researchers working in the field of new drug finding against SARS-CoV-2. However, a proper in-vivo and in-vitro rigorous research works are to be performed for the validation of our simulation work so that our recommended combination drug may be considered as a promising candidate for the drug design against COVID-19.

Declarations
Funding: Not applicable

Conflicts of interest/Competing interests:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.