key: cord-0829630-6s0i6vpb
authors: Mishra, Deepak; Raman Maurya, Radha; Kumar, Kamlesh; Munjal, Nupur S.; Bahadur, Vijay; Sharma, Sandeep; Singh, Prashant; Bahadur, Indra
title: Structurally modified compounds of hydroxychloroquine, remdesivir and tetrahydrocannabinol against main protease of SARS-CoV-2, a possible hope for COVID-19: Docking and Molecular Dynamics Simulation studies
date: 2021-04-16
journal: J Mol Liq
DOI: 10.1016/j.molliq.2021.116185
sha: 570c6d5304bc63ede603aac9c38f6e3d6980e9f3
doc_id: 829630
cord_uid: 6s0i6vpb

Now a days, more than 200 countries faces the health crisis due to epidemiological disease COVID-19 caused by SARS-CoV-2 virus. It will cause a very high impact on world’s economy and global health sector. Earlier the structure of main protease (Mpro) protein was deposited in the RCSB protein repository. Hydroxychloroquine (HCQ) and remdesivir were found to effective in treatment of COVID-19 patients. Here we have performed docking and molecule dynamic (MD) simulation study of HCQ and remdesivir with Mpro protein which gave promising results to inhibit Mpro protein in SARS-CoV-2. On the basis of results obtained we designed structurally modified 18 novel derivatives of HCQ, remdesivir and tetrahydrocannabinol (THC) and carried out docking studies of all the derivatives. From the docking studies six molecules DK4, DK7, DK10, DK16, DK17 and DK19 gave promising results and can be use as inhibitor for Mpro of SARS-CoV-2 to control COVID-19 very effectively. Further, molecular dynamics simulation of one derivative of HCQ and one derivative of tetrahydrocannabinol showing excellent docking score was performed along with the respective parent molecules. The two derivatives gave excellent docking score and higher stability than the parent molecule as validated with molecular dynamics (MD) simulation for the binding affinities towards Mpro of SARS-CoV-2 thus represented as strong inhibitors at very low concentration.

In the Wuhan city, the capital of Hubei province (China), since late December of 2019, firstly a case of pneumonia like viral infection attracted the attention of scientists, academicians and pharmaceuticals worldwide [1] [2] [3] [4] . It was an outbreak of pneumonia like symptoms with abnormal lungs, due to the anonymous reasons with an epidemiological among humans [2] . Earlier the terms like "the new coronavirus" and "Wuhan coronavirus" were used in a very common manner by January 2020. On the 11 th February 2020, this virus strain, earlier termed as '2019-nCoV & Wuhan virus' was officially announced a taxonomic name "severe acute respiratory syndrome coronavirus-2" (SARS-CoV-2) [5] due to the resemblance of the RNA genome of the previously identified coronavirus (SARS-CoV). On the same day, finally the world health organization (WHO) referred the name of this disease as COVID-19 [6] . This virus represents a pandemic threat to global public health and has already infected millions of the people and leads to more than hundred thousands of deaths globally in more than 200 countries. These numbers are increasing rapidly as of this manuscript writing. If the rapid spread of this virus-SARS-CoV-2 will be not controlled, it will probably bring a foremost challenge for global health system and world's economy [7] [8] [9] .

Coronaviruses (CoVs) are comparatively larger in size than other viruses and shared a >79% sequence identity to SARS-CoV [10] . Coronaviruses (CoVs) are enveloped single stranded positive sense RNA viruses and their family can be subdivided into four groups-alpha-, beta-, gamma and delta-CoV [11, 12] . Before the identification of SARS-CoV-2, only six of the previously known CoVs were cause mild to severe illnesses in humans. The alpha coronavirous (HCoV-229E and HCoV-NL63), cause milder respiratory disease in the upper respiratory tract. However the beta form of CoVs (HCoV-OC43, HKU1, and SARS-CoV are responsible for severe acute respiratory syndrome [12, 13] . While it's another form MERS-CoV (Middle east resipartory syndrome coronavirus) is responsible for infection in the lower respiratory tract [13] . The coronavirus has its appearance like crown shape ( Figure 1 ) and its structure is quite complex which contains four structural proteins [14] spike proteins (S), envelope protein (E), membrane protein (M) and nucleocapsid proteins (N) and several unidentified nonstructural protein. Recent studies shows that the spike protein and chymotrypsin-like protease (3CL pro ) also called as main protease (M pro ) plays an important role in mediating viral replication and transcription [15] [16] [17] [18] and hence is an attractive target for drug development for SARS-CoV-2. The 3CL pro catalyze the cleavage of polyprotein at 11 distinct sites and produces a various nonstructural protein which is responsible for viral replications. Therefore by the inhibiting the activity of this enzyme the researchers can cure the viral replication [17, 18] . The coronavirus can infect both animals as well as humans. Few animals like bat, hosts various types of coronaviruses appears to be immune to coronavirus prompted sickness [19] . Latest establishment shows that SARS-CoV-2 infects the lower respiratory system and results in viral pneumonia. Along with viral pneumonia SARS-CoV-2 may also attacks at kidney, heart, liver, gastrointestinal and central nervous system causing multiple organ failure [20, 21] . Current advancement about its transmission suggests that this new coronavirus strain, SARS-CoV-2 is more contagious than another coronavirus family [22] . Therefore structure-based virtual computational studies will be used to investigate the best drugs and is a fundamental approach to design highly potent inhibitor which can cure the viral replication.

The crystal structure of main protease (M pro ) of SARS-CoV-2 is determined and deposited in RCSB with PDB ID: 6LU7 [23] Figure 2 . According to our recent literature studies, most of the potential drugs contains heterocyclic moieties and was found to be active against SARS-CoV-2 M pro . Recently antiviral drug Remdesivir [24, 25] and antimalerial drug hydroxychloroquine (HCQ)/chloroquine [26] [27] [28] [29] were found to inhibit effectively to control SARS-CoV-2 in vitro and are currently being used in the treatment of patients with SARS-CoV-2 infection ( Figure 3 ). HCQ is less toxic than chloroquine and also found to be a potent inhibitor against the SARS-CoV-2 in vitro [30] . Herein, in order to manage/control this pandemic situation we have suggested few modified chloroquine, remdesivir and tetrahydrocannabinol (THC) and studied their inhibition abilities towards M pro using molecular docking and were further validated by molecular dynamics (MD) simulation study. In order to understand the properties of molecular assemblies in terms of their structure and the microscopic interactions between them computer simulations plays a significant role. Molecular dynamics simulation provides a route to dynamical properties of a system; transport coefficients, time-dependent responses to perturbations, spectra, rheological properties [31] .

Molecular dynamics simulation comprised of numerical, step-by-step solution of the classical equations of motion that for a simple atomic system can be written as equation 1.

For this we need to calculate the forces f i, acting on the atoms, and these are derived from a potential energy μ (r N ), where r N = (r 1 , r 2 , r 3 …..r N ) representing the complete set of 3N atomic coordinates.

Molecular dynamics simulations in drug discovery play an important role as they provide insight into the protein motion. The static models provided by homology modelling, NMR, X-ray crystallography gives insights about the macromolecular structure while drug binding and molecular recognition are very dynamic processes that can be studied using molecular dynamics simulations [32] . MD simulations are often used to explore the conformation space of molecules, especially the large molecules such as proteins. In MD simulation the energy surface is studied using the newton's law of motion for system. In the docked structures MD simulation is useful to confirm whether the major contacts found in docking are maintained during the MD simulation hence presenting the more reliable results. In other words the static view presented by docking is verified by the MD simulation.

Generally by combining the two in-silico techniques, docking and MD simulation will provide the more reliable results. 

The step-wise methodology followed for in-silico interaction pattern study of hydroxychloroquine (HCQ), remdesivir, and tetrahydrocannabinol (THC) and their derivatives for characterizing the potent drug molecules targeting M pro /3CL pro protein of SARS-CoV-2 is described as follow.

A dataset of eighteen derivatives of hydroxychloroquine, remdesivir, and tetrahydrocannabinol were designed ( Figure 4 ). The initial structures of all the derivatives were drawn in Gauss view 3 .0 and all the molecular structures of derivatives along with the parent drugs were optimized at the level of parametrization method 6 (PM6) semi-empirical method in Gaussian 09, quantum chemistry software [33] [34] [35] . M pro /3CL pro protein of SARS-CoV-2 has a catalytic dyad comprises of Cys-His and substrate binding site is located in the cleft between domain I and II [37] . [23] . In the current study chain A was taken for macromolecule preparation. Chimera 1.13.1 software was used for carrying out the minimization of protein structure implementing Amber ff99Sb force field [38] .

The M pro active site mainly contains sub-sites commonly named as S1, S2, S1', S2', S4. The amino acid residue Phe140, Leu141, His163, Met165, Gln166 and His172 are present in S1 site, while Met49 and Asp187, Glu189 are present in bulky hydrophobic subsite S2 [39] . The amino acid residue Thr25, Leu27, Cys38, Pro-39, Val42, and Cys145 present in S1' subsite and play an important role in catalysis [18] . The Thr-26, Asn-28, Tyr-118, Asn-119, and Gly-143 are found in S2,'while the flexible S4 site contains Leu167, Gln192 [39, 40] . The active site was predicted by CastP software active site prediction tool [41] .

Dataset of optimized molecular structures of hydroxychloroquine, remdesivir, THC and its derivatives were subjected to molecular docking against the M pro of SARS-CoV-2 in the above defined active site using AutoDock 4.2.2 [42] and AutoDock Vina [43] . AutoDock implies 3dimemsional potential grids to calculate the potential energy of docking molecules. Grid 

To study the stability of protein-ligand complex i.e., The molecular dynamics (MD) simulation study [44] , we have used the Gromacs 5.1.2 software which is widely used. The MD simulation of the hydroxychloroquine and tetrahydrocannabinol and their derivative DK7 and DK16 respectively those were showing the best binding affinity from the docking study against M pro was carried out. All the four systems were employed for 100 ns time scale simulation. Protein topology was generated by GROMOS 96 54a7 forcefield [45] in Gromacs while the ligands topology was obtained using automated topology builder (ATB) server [46] . Solvation of proteinligand complex was carried out in a cubic box of size 9.73 X 9.73 X 9.73 nm using SPC water model. Systems were neutralized by the addition of four Clions and then the energy minimization of the systems was carried out using the steepest descent algorithm. Particle-Mesh Ewald method [47] was used for the calculation of electrostatic interactions of the systems. A normal cut-off of 0.9 nm was implied for Van der Waals interaction. Constraints imposed by the covalent bonds including hydrogen atoms of the system were calculated by a linear constraint solver for molecular simulations (LINCS) algorithm. NVT (Number of particles, Volume and Temperature)

and NPT (Number of particles, Pressure and Temperature) simulations were run on the system for 1 ns using Berendsen thermostat method [48] to control the temperature at 300 K with a coupling time of 0.1 ps and pressure was restrained at 1 bar [49] . The co-ordinates were saved in 0.002 fs for all the systems [50] . Different analytical tools such as gmx rms, gmxrmsf, gmx gyration, gmxh_bond, and gmxsasa were implied for the calculation of RMSD, RMSF, Rg, Hydrogen bonds and SASA respectively. The RMSD and RMSF trajectories were also resolved with the help of visualization software like Visual Molecular Dynamics (VMD) [51] and Chimera 1.13.1 [38] visualization tools and the graphs were plotted using the Origin software.

Analysis of eigenvectors and eigenvalues was performed using Principal component analysis (PCA) or Essential dynamics method in Gromacs for obtaining the projection of first two Principal Components (PCs). Protein function is majorly defined by the associated motion of the protein which was obtained by PCA [52] . All the rotational and translational movements were removed and co-variance matrix was formed. The atomic coordinates from the positional covariance C and its eigenvectors were used in the following equation for the calculation of elements of the positional covariance matrix C:

The i th C α atom was represented by the cartesian coordinate q i and C α atom was denoted by N in the protein-ligand bound systems. Complete translational and rotational movements using 'Least square method' of the equilibrated trajectory was removed that was further superimposed on a reference structure. For the prediction of rest of eigenvectors and eigenvalues λ i , all the matrices were diagonalized using orthogonal coordinate transformation matrix Λ.

In the equation, eigenvectors corresponding to the direction of motion relative to <q i > were represented in the columns and each eigenvector associated with the eigenvalue represents the total mean-square fluctuation of the system ahead the corresponding eigenvector. gmxcovar and gmxanaeig tools of Gromacs were used for the calculation of eigenvector and eigenvalues for the obtained last 20 ns trajectories.

Since chloroquine based antimalarial drug HCQ and remdesivir, the antiviral drug is in use for the treatment of COVID-19. Few reports suggest that cannabis and their derivatives have many therapeutic applications including antiviral activity [53] [54] [55] [56] [57] . Herein, we have done structural modification in the HCQ/chloroquine, remdesivir and THC ( Figure 4) We have also compared the free binding energy and inhibition activity of DK7 -7.56 kcal mol -1 with its isomer DK8 (-6.78 kcal mol -1 ), surprisingly, it was a more than three-fold change in inhibition constant against M pro . We further investigate the effect of more polar -COOH functional group of pyrrolidine DK10. The molecular docking data reveals that the activity of DK10 (-7.13 kcal mol -1 ) is (c) more than that of DK8 but less than that of DK7. We have also compared the effect of introducing one more OH group in pyrrolidine ring (adjacent and syn to -OH group in DK7), in DK5 (free binding energy -6.01 kcal mol -1 , and it was found to be less active than DK7. It was also observed that the replacement of pyrrolidine moiety with planar 1,2,3 triazole derivative (DK20 and DK21, see table   1 ) has no significant effect on the activity of the molecule against M pro of SARS-CoV-2.

Further, molecular dynamics (MD) simulation which is a computerized method for studying the physical movements of the atoms while binding [58] was carried out for HCQ-M pro and DK7-M pro complexes. MD simulation of the docked complexes for 100 nano second (ns) was carried out and equilibration state for both the systems was obtained. The trajectories were observed for RMSF, Rg, SASA, and H-bonds for last 20 ns equilibrated trajectories including the binding free energy analysis also for last 20 ns.

Dynamic stability of all the systems is described by the RMSD, it calculates the changes in the 

Protein-ligand stability is majorly defined by the hydrogen bonds as they contribute to transient interactions in the complex formation providing stability to the docked complex. In this study, we 

In the molecular dynamics simulation, the principal component analysis predicts the correlated motion of the complexes, in the calculation of the Principal Components (PCs) the overall motions in the proteins is described by a few key eigenvectors. Hence in the study, we had plotted the eigenvalues with the eigenvectors for the 10 eigenvectors in (Figure 10a ). The first ten eigenvectors 

In our second studies we have designed a conjugate of hydroxychloroquine with remdesivir.

Several report suggest that remdesivir is effectively in use with HCQ for the treatment of COVID-19

patient. So, we planned to design a hybrid of HCQ and remdesivir. These hybrid molecules (DK2, DK3, DK6, DK9 and DK12) were docked with active site of the M pro of SARS-CoV-2. The docking score obtained from the molecular docking studies shows, among these conjugates DK3, DK6, DK9 and DK12 are little better but not that much active as compared to HCQ alone while DK2 (free binding energy -4.22 kcal.mol -1 ), was found to have lower free binding energy and poor inhibition with respect to HCQ against M pro of SARS-CoV-2 (Table 1) .

Since tetrahydrocannabinol is useful in the cure of various disease [53] [54] [55] [56] [57] . Hence, taking THC as a parent moiety (DK15, Figure 11 DK17, DK18 and DK19) it was found that all derivatives are more active than the parent molecule (DK15). Among all these the best result was obtained by DK16, DK17 and DK19 (-7.64 kcal mol -1 , -7.62 kcal mol -1 , -7.42 kcal mol -1 respectively). It has also been found that DK16 is more potent inhibitor for M pro of SARS-CoV-2 than the HCQ's derivative DK7. The molecular docking of DK16

suggests that the carboxylic acid group present at parent moiety makes a hydrogen bond with Lys102 and oxygen atom of tetrahydropyran ring is involved in hydrogen bond formation with Gln110. While the other amino acid residue Val104, Ile106, Gln107, Thr111, Asn151, Ile152, Asp153, Ser 158, Phe294 are involved in hydrophobic interaction with the different residue of the moiety ( Figure 12 ). 

Hence molecular dynamics simulation for 100 ns of the docked complex of M pro with THC and DK16

was carried out to further validate the results obtained from the docking studies. In the study, equilibration state for both the systems were obtained and last 20 ns equilibrated trajectories were observed for RMSF, Rg and SASA H-bonds including binding free energy analysis for last 20 ns.

Average showed average RMSF value less than the THC representing that ligand complex is a stable complex as compared to the parent drug (Figure 13a and 13b) . 

Number 

Average 

In 

The binding ability of the ligands to the protein is assessed with the help of binding free energy calculation which is the sum of non-bonded interaction energies of the complex. MMPBSA tool [59] that is directed by the Gromacs, was adopted for the calculation of binding free energy. All the binding free energies of the complexes are summarized in Table 2 , showing the stability of the complex formation with the derivatives as compared to the parent molecule. 

Since COVID-19 is epidemiological and created a globally challenged health crisis, till date there is no any particular vaccine or drug have been developed to control it. Initially it was found that hydroxychloroquine and remdesivir are little bit active and are using as a drug for management of COVID-19. Keeping it in mind we have designed derivatives of these drugs and studied about their inhibition and binding with M pro of SARS-CoV-2. In our docking studies we have concluded that DK4, DK7, DK10, DK16, DK17 and DK19 can be use as inhibitor for M pro of SARS-CoV-2 to control COVID-19 very effectively. Also, MD simulation for HCQ, DK7, THC, and DK16 those were highest scoring molecules obtained in the docking studies was carried out. RMSD, RMSF, hbonds, Rg, SASA, PCA parameters were calculated to validate the stability of the protein ligand complexes and found that DK16 showed higher stability as compared to the parent molecule as well as from DK7. DK7 molecule showed comparable stability as the HCQ with M pro protein of SARS-CoV-2. It is anticipated that both these derivatives will be potent drug candidates and can be further validated experimentally.

Titled "Structurally modified compounds of hydroxychloroquine, remdesivir and tetrahydrocannabinol against main protease of SARS-CoV-2, a possible hope for COVID-19: Docking and Molecular Dynamics Simulation studies.

"I, PROF. Indra Bahadur (the Corresponding Author), declare that this manuscript is original, has not been published before and is not currently being considered for publication elsewhere. I further confirm that the order of authors listed in the manuscript has been approved by all of us.

There is also no conflict of interest in any way.

Sincerely,

February 16, 2021

Pneumonia of unknown etiology in Wuhan, China: potential for international spread via commercial air travel

The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health-the latest 2019 novel coronavirus outbreak in Wuhan, China

A new coronavirus associated with human respiratory disease in China

A pneumonia outbreak associated with a new coronavirus of probable bat origin

The species severe acute respiratory syndrome-related coronavirus: classifying 2019-nCoV and naming it SARS-CoV-2

Evidence of the COVID-19 virus targeting the CNS: tissue distribution, host-virus interaction, and proposed neurotropic mechanisms

Severe acute respiratory syndrome-related coronavirus: The species and its viruses-a statement of the Coronavirus Study Group

Will novel virus go pandemic or be contained?

Coronavirus is now expected to curb global economic growth by 0

Network-based drug repurposing for novel coronavirus 2019-nCoV/SARS-CoV-2

Coronaviruses post-SARS: update on replication and pathogenesis

Coronaviruses: an overview of their replication and pathogenesis

Coronavirus entry and release in polarized epithelial cells: a review

Analysis of therapeutictargets for SARS-CoV-2 and discovery of potential drugs bycomputational methods

Therapeutic drugs targeting 2019-nCoV main protease byhigh-throughput screening

Nelfinavir was predicted to be a potential inhibitor of 2019-nCovmain protease by an integrative approach combining homology modelling,molecular docking and binding free energy calculation

Structure of M pro from SARSCoV-2 and discovery of its inhibitors

Structures of two coronavirus main proteases: Implications for substratebinding and antiviral drug design

Global patterns in coronavirus diversity

Epidemiology, genetic recombination, and pathogenesis of coronaviruses

A novel coronavirus from patients with pneumonia in China

An updated estimation of the risk of transmission of the novel coronavirus (2019-nCov)

Structural and evolutionary analysis indicate that theSARS-CoV-2 Mpro is a challenging target for smallmolecule inhibitordesign

Currentpharmacological treatments for COVID-19: What's next?

Favipiravir, an anti-influenza drug againstlife-threatening RNA virus infections

A systematic review on the efficacy and safety of chloroquine for thetreatment of COVID-19

Screening ofchloroquine, hydroxychloroquine and its derivatives for their bindingaffinity to multiple SARS-CoV-2 protein drug targets

Hydroxychloroquine and azithromycinas atreatment of COVID-19: Results of an open-label non-randomizedclinical trial

Remdesivir and chloroquine effectively inhibit the recently emerged novel coronavirus (2019-nCoV) in vitro

Hydroxychloroquine, a less toxic derivative of chloroquine, is effective in inhibiting SARS-CoV-2 infection in vitro

Molecular dynamics simulation for all

Molecular dynamics simulations and drug discovery

Gaussian 09 (Reversion B. 01)

Molecular structure, vibrational assignments and non-linear optical properties of 4, 4'Dimethylaminocyanobiphenyl (DMACB) by DFT and ab Initio HF Calculations

Validation of semiempirical methods for modeling of corrinoid systems

Structure of Mpro from COVID-19 virus and discovery of its inhibitors

Structure of coronavirus main proteinase reveals combination of a chymotrypsin fold with an extra α-helical domain

UCSF Chimera-a visualization system for exploratory research and analysis

Structure of mainprotease fromhuman coronavirus NL63: Insights for wide spectrumanti-coronavirus drug design

Identification of potential molecules against COVID-19 main proteasethrough structure-guided virtual screening approach

CASTp 3.0: computed atlas of surface topography of proteins

Automated Docking of Flexible Ligands: Applications of AutoDock

AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization and multithreading

GROMACS: Fast, Flexible, and Free

Gunsteren Definition and testing of the GROMOS force-field versions 54A7 and 54B7

An Automated Force Field Topology Builder (ATB) and Repository: Version 1.0

Accuracy and efficiency of the particle Mesh Ewald method

Molecular dynamics with coupling to an external bath

Molecular dynamics simulations at constant pressure and/or temperature

A leap-frog algorithm for stochasticdynamics

VMD: Visual molecular dynamics

Points of significance: Principal component analysis

Possible therapeutic applications of cannabis in the neuropsychopharmacology field

Medicinal cannabis for psychiatric disorders: a clinically-focused systematic review

History of cannabis and its preparations in saga, science, and sobriquet

The religious and medicinal uses of Cannabis in China, India and Tibet

Potential of cannabidiol for the treatment of viral hepatitis

High performance molecular simulations through multi-level parallelism from laptopsto supercomputers

Open Source Drug Discovery Consortium, & Lynn, A.g_mmpbsa--A GROMACS tool for high-throughput MM-PBSA calculations

The authors are thankful to Scfbio, IIT Delhi for providing the facility of Gaussian software to optimize the structures and also supporting in carrying out molecular dynamics simulation studies.