Hydroxychloroquine and Azithromycin Molecular Action against SARS-CoV-2 Viral Protein: A Molecular Dynamic Study

Abstract

For the past few months, the world has gone through hell with the emergence of the SARS-CoV-2 virus and the resulting pandemic. Faced with this disease, various therapeutic strategies have been developed to understand how to eradicate this virus. Here we present a molecular dynamics simulation study on the effect of a dual therapy (hydroxychloroquine and azithromycin) on the open and closed forms of a viral protein. We show in particular that hydroxychloroquine has no significant interaction with the viral receptor-binding domain RBD when it interacts with its host receptor. However, this molecule can, in the closed form of the virus, block the movement of these receptors and thus prevent the attachment of the virus to the host cell. The azithromycin molecule interacts very well with the open receptor but can also be inserted into the S₂ domain of the protein. It therefore presents two potential mechanisms of action against the virus, mainly on the closed state of the viral protein.

Keywords: Hydroxychloroquine; Azythromicin; SARS-CoV-2 structure; MD-Simulation; Docking

Introduction

In 2002, the first emergence of a pathogenic coronavirus revealed to the world the possibility of a pandemic. The severe acute respiratory syndrome coronavirus, or SARS-CoV, was responsible for very important breathing syndromes [1-4] with, however, a small amount of death around the world (8000 persons) while the mortality rate reached 10%. More recently (2012), a severe pneumonia appeared in Saudi Arabia due to a novel coronavirus [5,6]. This one, called MERSCoV for Middle East Respiratory Syndrome Coronovirus, still exists but concerns only the Arabic peninsula. While very localized on a small area, this MERS-Cov is highly dangerous since its mortality rate is about 35% (1 over 3 patients). These two cases are not the only ones since periodically, other coronaviruses, while less virulent are appearing [7,8,9].

These first viral apparitions should have alerted us to the possibility, in the long term, of the emergence of a more virulent coronavirus clearly difficult to manage. At the end of 2019, a new coronavirus called severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), or COVID- 19 disease, developed as a human pathogen in a Chinese city (Wuhan). Despite strong resolutions in China and over several states (locking up millions of people), global economic and touristic development is leading to a general spread of the disease. Although SARS-CoV-2 has many points in common with SARS-CoV and MERS-CoV, it appears to be transmitted more easily and much faster than SARS-CoV [10,11] and MERS-CoV [12,13]. To date, more than 64 million cases of COVID-19 have been confirmed and at least 1.400,000 deaths have been recorded worldwide by the World Health Organization which declared the first real pandemic of March 21, 2020. Although all the cases have not been identified, these figures lead to a mortality rate close to 2.1%, i.e. 4 times the rate of the seasonal flu. While the peak of the pandemic seems to be behind us, the second wave forces us to find a therapeutic method to treat this virus.

Many drugs have been clinically tested in numerous clinical projects (“discovery” for France, “recovery” for England”) but no one has determined a truly effective treatment against the SARS-CoV-2 virus. In parallel, a protocol, confirmed by many other studies [14- 17], seems to be the most appropriate for combating the virus and reducing the degree of contagiousness. This protocol, that should be administered as soon as the first symptoms appear, combines both the antimalarial drug Hydroxychloroquine (HCQ) and the antibiotic Azithromycin (AZM). The results of the various studies lead to a sharp decrease of the mortality rate (under 0.5%). However, many other analyzes question the results of the therapies. In order to better understand the role of this double drug treatment against this new virus, an analysis of its action is necessary at the molecular level [18].

With regard to its genome sequence, SARS-CoV-2 belongs to the same beta-coronavirus family as SARS-CoV and MERS-CoV. These coronaviruses have a spherical envelope with a diameter close to 100 nm. The latter is composed of a N nucleoprotein surrounded by a lipid bilayer originating from the host cell. Three other proteins are then found on this surface, protein S (or spike protein), protein M (or membrane protein), and protein E (or envelope protein). Homotrimerization of S proteins [11] on the surface of the virion is the key step in viral infection. To decrease the viral progression, the drug must therefore target the S protein. However, the latter is separated into two domains which have a very specific role [19]. The S₁ domain, mainly formed by a Receptor Binding Domain (RBD), binds to the host cell receptor [20] (called angiotensin converting enzyme 2 (ACE2) [21,22] while that the S₂ domain is at the origin of the fusion of protein E with the host cell [19].

Therefore, the main objective of any viral treatment is to block the binding of the RBD domain to the ACE2 receptor in order to avoid the fusion of the virus with the host cell. Note that the drugs could prevent protein fusion with the host cell through structural changes.

Here we propose to determine what are the main sites of interaction of the HCQ and the AZM drug molecules on structural S protein. Our work, based on several molecular dynamics simulations will be separated into two parts. We will first determine the binding of each drug when RBD is associated with the ACE2 receptor. Then, in a second step, we will show that these molecules can also have another target to on the state close to the S protein.

Method

Hydroxychloroquine (HCQ) is (RS)-2-[[4-[(7-chloroquinoline- 4-yl)amino]pentyl](ethyl)amino] ethanol. Its 3D structure has been obtained through the Pub-Chem CID: 3652 file. Azithromycin (AZM) is (2R, 3S, 4R, 5 100 R, 8R, 10R, 11R, 12S, 13S, 14R)-11- [(2S, 3R, 4S, 6R)-4-(dimethylamino)-3-hydroxy-6-methyloxan-2-yl] oxy-2-ethyl-3, 4, 10-trihydroxy-13- [(2R, 4R, 5S, 6S)-5-hydroxyl- 4-methoxy-4,6-dimethyloxan-2-yl] oxy-3, 5, 6, 8, 10, 12, 14 - heptamethyl-1-oxa-6-azacyclopentadecan-15-one. Its 3D structure has been obtained through the Pub-Chem CID: 447043 file.

Our strategy was organized in two different steps. First of all, we focused on the simulation of the Receptor Binding Domain (RBD) bound to the ACE2 receptor of the host cell with the HCQ or the AZM molecules. The goal of these first calculations is to observe whether the molecules of the drug can interact directly with the viral protein when they are attached to its host cell. To simulate such arrangement of proteins, we use the 6M0J pdb structure. The relaxed crystal structures approaching the living organism as well as the effect of glycosylation on the stability of the structure were studied.

Then, in a second step, the full conformation of the spike glycoprotein trimer SARS-CoV-2 was simulated in presence of HCQ or AZM molecules. Its structure was obtained from pdb file #6VXX. It has a resolution of 2.80Å as determined from electron microscopy. It is composed of 3 chains intercalated with different respective domains such as the NTD (N-terminal domain) and the RBD (receptor binding domain) which belong to the S₁ part of the protein.

Classical MD simulations were performed by constructing the full molecular force field for HCQ and AZM using the SwissParam Force Field Toolkit package [23,24].

For the protein, the molecular force field was constructed according to the CHARMM-GUI procedure in order to appropriately relax the different parts of the protein [25,26]. N-glycosylation of the proteins, when necessary, was achieved using the CHARMMGUI Glycolipid Modeler.[27] The structure of these proteins was first minimized and then progressively equilibrated (1ns) and run (40ns) using MD simulations (NAMD 2.12 package) for a total of 41 ns in saline solution media [28]. All the structures of proteins and molecules are shown in Figure 1a-d.

Figure 1: Molecular structure of a) hydroxychloroquine, b) azithromycin, c) RBD (red)+ACE2 (blue) complex and d) 6VXX protein (note that the RBD, NTD and S2 parts 11of the protein were depicted in van der Waals, licorice and ribbon+lines modes for the three monomers shown in different colors).

    

    

    Figure 1: Molecular structure of a) hydroxychloroquine, b) azithromycin, c)
RBD (red)+ACE2 (blue) complex and d) 6VXX protein (note that the RBD,
NTD and S2 parts 11of the protein were depicted in van der Waals, licorice
and ribbon+lines modes for the three monomers shown in different colors).
    

Close

The systems (protein or protein+drug) were then combined and finally solvated in a water box large enough to prevent interaction of each entity with its neighbor in adjacent cells when periodic boundary conditions are used. To mimic a salt medium, NaCl ions (at a concentration of 0.15M) were added to the water model (transferable intramolecular potential with 3 points, TIP3P). CHARMM36 forcefield optimization parameters are used in all simulations [29]. The complete systems contained 245805, 245757, 301911 and 301896 atoms, for HCQ/RBD, AZM/RBD, HCQ/6VXX and AZM/6VXX systems, respectively.

All simulations were performed under constant temperature and pressure with a 13Å cutoff for unbonded forces. The temperature is fixed at 310K (Langevin dynamics) [30] and the pressure at 1atm (Langevin piston), [31] respectively. Long range electrostatic forces are evaluated using the classical Particle Mesh Ewald (PME) method [32] with a grid spacing of 1.2Å, and fourth-order spline interpolation. The simulations were carried out in three stages: first, a minimization stage (steepest descent) of 5000 steps was carried out to remove the strongest atomic hindrance from the system. Then, an equilibration phase was imposed on the system where the movement of the protein were constrained in order to stabilize the drug, the solvent and the ions around the protein. The harmonic force constants on the protein backbone and the sidechain were chosen to be respectively equal to 400kJ mol^-1 nm^-2 and 40kJ mol^-1nm^-2, for a total duration of 1ns. Finally, the production step was performed without position restrains for a total of at least 60ns (time step of 1fs). Each simulation uses periodic boundary conditions in the three directions of space and the list of neighbors is refreshed every 12-time steps. No constraint is imposed during the production phase of MD simulations. Note that for system equilibration, all bond lengths involving hydrogen atoms were fixed using the SHAKE algorithm [33]. Consequently, our results are obtained with all the atoms left free in the simulation box.

To determine the sites of interaction in such huge systems, several simulations were run with different starting configurations. These configurations were obtained through a docking procedure between the viral protein and the drugs. AutoDock Vina (a fast and accurate evolution of AutoDock) was used as the molecular docking engine, being able to process a massive number of ligand positions in a relatively short time. The different starting configurations used in the molecular dynamics simulations were thus chosen with regard to the best scoring functions obtained thanks to the optimization algorithm. When the drug interacted with the 6M0J, at least 7 different simulations were run for each drug. For the 6VXX system, only 3 simulations were performed due to the huge number of atoms in the system. Note that for each system, we also performed an additional simulation with a drug position located outside the protein surface in order to study the path to the active site of the protein.

Results

Relaxations of the proteins and the molecules in saline solution

The 6M0J and 6VXX structures were relaxed in their respective solvated plus ionized solvent box via MD simulations. As with half of the proteins, the glycosylation of the protein should also be studied [34]. Indeed, this could have an importance on how the virus adapted during the infection of the host cell [35]. Its impact on S protein has been examined for Sars-CoV with mainly N- glycosylation type versus O-glycosylation [36]. Different mutagenesis analyzes [37] have identified N-glycosylation sites. For 6M0J, the latter were determined but only four sites (3 on ACE2 and 1 on RBD chain) were accessible by GlyProt [38] (90, 322, 546 for ACE2 and 343 for RBD). For 6VXX, we can find 22 possible N-linked glycosylation sites per monomer and only 13 are conserved in the pdb structure. In this context and amongst all the possibilities of glycosylation (sialylation, fucosylation, …), we have limited our studies to branches of N-glycosylation formed by association of N-acetyl-glycosamine (BGLC) and Mannose (BMAN and AMAN) groups such as --ASN-BGLC-BGLC-BMAN-(AMAN)₂.

For each protein the progressive modification of the structure was followed by calculations of its root mean square deviation. The results are shown in Figure 2a-b and the superimposition of relaxed glycosylated and non-glycosylated proteins are shown in Figure 2c-f.

Figure 2: a) RMSDs of 6M0J structures (with non-glycosylated sites (black curve) and with glycosylated ones (red curve)). b) Same but for 6VXX structures. c-d) Superimposed crystallized (blue) and relaxed (red) structures for 6M0J. (c) and 6VXX (d) proteins, respectively. e-f) Superimposed crystallized (blue) and relaxed (red) structures for glycolyzed 6M0J. (e) and 6VXX (f) ) Proteins, respectively.

    

    

    Figure 2: a) RMSDs of 6M0J structures (with non-glycosylated sites (black
curve) and with glycosylated ones (red curve)). b) Same but for 6VXX
structures. c-d) Superimposed crystallized (blue) and relaxed (red) structures
for 6M0J. (c) and 6VXX (d) proteins, respectively. e-f) Superimposed
crystallized (blue) and relaxed (red) structures for glycolyzed 6M0J. (e) and
6VXX (f) ) Proteins, respectively.
    

Close

As shown in Figure 2a-b, the simulations ended as soon as the RMSD and the total energy converged for a reasonable simulation time thanks to the system size. The change in structure is similar, although slightly higher for glycosylated proteins compared to nonglycosylated proteins. These slight differences mainly come from the glycan groups which are free to move around in the solvent and can assume random positions during their interactions with water molecules. The 6M0J RMSDs converge towards a value close to 2.5Å while the 6VXX ones tend to 4.0Å. The difference in molecular weight between the two respective proteins could explain this RMSD difference. To compare the final structures obtained after the relaxation of the proteins, we describe the superimposed structure between the crystallized and relaxed structures in Figure 2c-f. As expected, the proteins did not show strong changes in their backbone and no significant folding or unfolding.

HCQ and AZM were also relaxed under the same conditions of solvation in order to start the simulations using both drugs and proteins with optimized structures. According to the analyzes of their RMSD, the modification of their structure did not undergo a value greater than 2.0Å, which is very low. We can use these structures to study their interaction with the different proteins.

Interactions of the RBD-ACE2 complex with the therapeutics

Our initial investigation was devoted to the action of drugs on the Receptor Binding Domain (RBD) bound to the ACE2 receptor exhibiting N-glycosylation. Indeed, this part of the Sars-Cov-V2 virus is dedicated to binding the virus to the host cell. Most studies aim to develop specific treatments the goal of which is to block the interaction between the virus (RBD) and the host cell (ACE2).

Due to the large size of the system, docking simulations were first performed in order to determine the most relevant interaction sites between the drugs and the protein. Seven configurations were then carried out in the case of the HCQ/RBD systems with scoring functions varying from -7.78 to -6.73. For AZM/RBD systems, nine systems were obtained with stable scoring functions, i.e. ranging from -9.94 to -7.52. Based on these different configurations, molecular dynamics simulations in full solvent were performed in order to test the stability of the drug under biological conditions. For these simulations, a 1ns equilibration phase was performed before running some 10ns production simulations. The energy between the drug and the protein was then determined to obtain the most stable site of interaction by molecular dynamics simulation.

As shown in Figure 3, the different sites of interaction determined by simulations can belong to several parts of the protein. RBD (red part of 6M0J) is however at the origin of the best pair interaction energies for each molecule. We can also observe, from these data, that the AZM molecules tend to interact strongly with the 6M0J protein compared to HCQ. We did not observe any specific residue in the adsorption configurations.

Figure 3: The three best stable configurations obtained through molecular dynamics simulations. a,b,c) for AZM; d,e,f) for HCQ. For each configuration, the mean interaction energy and the protein residues interaction with the molecule are given to precise the differences between each system.

    

    

    Figure 3: The three best stable configurations obtained through molecular
dynamics simulations. a,b,c) for AZM; d,e,f) for HCQ. For each configuration,
the mean interaction energy and the protein residues interaction with the
molecule are given to precise the differences between each system.
    

Close

We finally carried out a final simulation where the drugs were placed far from the adsorption sites of the protein (15Å). The complete systems were gradually relaxed in order to study the diffusion of the drug to this part of the protein. The two molecules (protein and drug) being relaxed in independent simulations, the equilibration phase here was equal to 1ns, and the system was running for a production phase for 38ns. We extracted during this phase the pair interaction between drug and protein over time. Figure 4 shows different simulation times where we can observe the position of the drug and also the state of the protein. At the beginning of the simulation, the molecule is left 15.0Å from the junction between the RBD and ACE2. The drugs diffuse in the solvent due to different pair interactions and thermal agitation and come closer to the protein. As can be seen in the middle of the simulation, each drug has found a site where it can be stabilized. Figure 5a, which exclusively shows the pairwise interaction of AZM (in black) and HCQ (in red) with the protein, clearly demonstrates that the drug could be adsorbed on the protein. However, while the AZM molecule can arrange itself to find a better adsorption site on the protein, HCQ quickly desorbs and return to the solvent bath (Figure 4f), where it is confirmed by a pair interaction equal to 0 with the protein (Figure 5a). On the contrary, during the simulations, AZM never desorbed from the protein and found a stable site which is more stable than the first one (Figure 4c) and of the same order as those obtained in Figure 3.

Figure 4: Figure 4a-c) Sketch of AZM positions during the simulation (at t=0 ns, t=4 ns and t=19 ns, respectively). d-f) Sketch of HCQ positions during the simulation (at t=0 ns, t=7.4 ns and t=19 ns, respectively). The RBD is depicted using red ribbon while ACE2 is in blue ribbon The chains showed in color correspond to the N-glycan parts of the protein.

    

    

    Figure 4: Figure 4a-c) Sketch of AZM positions during the simulation (at t=0
ns, t=4 ns and t=19 ns, respectively). d-f) Sketch of HCQ positions during
the simulation (at t=0 ns, t=7.4 ns and t=19 ns, respectively). The RBD is
depicted using red ribbon while ACE2 is in blue ribbon The chains showed in
color correspond to the N-glycan parts of the protein.
    

Close

Figure 5: a) Pair interaction (smooth data in line representation), and b) RMSD of the RBD+ACE2 protein in presence of AZM (black) or HCQ (red) drugs.

    

    

    Figure 5: a) Pair interaction (smooth data in line representation), and b)
RMSD of the RBD+ACE2 protein in presence of AZM (black) or HCQ (red)
drugs.
    

Close

The RMSDs of the protein were also plotted during these two production phases to see the impact of the drug on the stabilized protein structure (Figure 5b). Regardless of the drug considered, the RMSD followed the same trend when approaching the drug. It should be noted however that a series of upward peaks is observed around t=14.4 ns, which corresponds to the arrivals of the HCQ on the protein. These small successive deformations could be at the origin of the HCQ desorption since for AZM case, we rather observed some small downward peaks. It should be noted that we have not observed any direct influence of the N-glycosylation functions toward the adsorption or the desorption of the HCQ and AZM molecules, although they are very labile when they are subjected to thermal agitation. They could however interfere with the accessible sites determined by the previous calculations by steric hindrance.

Interactions of the 6VXX trimer with the drugs

HCQ and AZM were then relaxed close to optimized 6VXX protein. Due to the large size of the protein, we only put the drug in two different situations. First, the drugs were placed in a cavity of the protein where the molecule could take place. This cavity is located between two RBDs of the trimer and was determined by docking calculations implemented in QuikVina-w [39]. The purpose of these first calculations is to determine if the drug could accommodate on the viral protein. Then, in order to simulate a more realistic diffusion of the drugs, we let them freely approach the molecules. Due to the size of the protein, and the time simulations, no additional configuration was performed as the production runs were longer than with RBDACE2 alone to converge.

When the drugs are placed in the cavity near the two RBD sites (Figure 6a-f), the interaction with the protein is initially slightly attractive (Figure 7a). Then, progressively the molecules tend to translate and find the way to fall into a well of significant attractive potential well. For AZM, this well is located near a N-Terminal Domain (NTD) shown in CPK representation in Figure 6. Two states of equivalent energies were explored by AZM which are around the same area as the final configuration observed in Figure 6c. The mean pair interaction energy obtained during the simulation time at this site is equal to -32±5 kcal/mol. In contrast, the HCQ molecule falls into its lower potential well (-37±3 kcal/mol) by translating between the two RBDs (Figure 6d-f). Then it eventually moves into an area where its pair interaction energy fluctuates slightly around its average value. The modification of the protein structure is directly impacted as soon as the drugs fall into their potential valley. Indeed, although slightly modified during the first ten ns of simulations, it sharply increases to reach, after 42 ns, a value of 4Å for AZM (6Å for HCQ) which is not stabilized yet (Figure 7b). We observed that the main deformation of the 6VXX is located at the bottom of the protein, near the intrusion zone. The progressive unfolding of this zone can be observed in the sketches shown in figure 6. Although very stable when placed alone in the solvent, this zone is particularly sensitive to the presence of the drug during the simulations. The estimate of the RMSF (root mean square fluctuations) of the protein confirms these observations since the RBD (NTD) zones have a maximum RMSF equal to 3.2Å (4.7Å) while the intrusion zone RMSF can reach a value of 8.3 Å.

Figure 6: a-c) Initial and final AZM positions during the simulation (at t=0 ns and t=42 ns, respectively), and protein residues in interaction with AZM at the final stage of the simulation. d-f) Initial and final HCQ positions during the simulation (at t=0 ns and t=42 ns, respectively) and protein residues in interaction with HCQ at the final stage of the simulation. The RBD, NTD and S₂ parts of the protein were depicted in van der Waals, licorice and ribbon+lines modes for the three monomers shown in different colors.

    

    

    Figure 6: a-c) Initial and final AZM positions during the simulation (at t=0
ns and t=42 ns, respectively), and protein residues in interaction with AZM
at the final stage of the simulation. d-f) Initial and final HCQ positions during
the simulation (at t=0 ns and t=42 ns, respectively) and protein residues in
interaction with HCQ at the final stage of the simulation. The RBD, NTD
and S2 parts of the protein were depicted in van der Waals, licorice and
ribbon+lines modes for the three monomers shown in different colors.
    

Close

Figure 7: a) Pair interaction (smooth data in line representation), and b) RMSD of the 6VXX protein in presence of AZM (black) or HCQ (red) drugs.

    

    

    Figure 7: a) Pair interaction (smooth data in line representation), and b)
RMSD of the 6VXX protein in presence of AZM (black) or HCQ (red) drugs.
    

Close

When the drugs are now placed far from the 6VXX structure (close to 15Å) (Figure 8a-f), the pair interactions with the protein are at first zero then gradually become attractive, ending up with a large stabilized plateau in the last 15ns (Figure 9a). For AZM, this well is located at the bottom of the trimer shown in ribbon+line representations in Figure 8. The approach of the AZM is quite slow. For 15ns, it did not interact with the protein. Then it finds an attractive interaction site and reaches its final stable state translating along the surface of the protein (as shown in Figure 6a-c). The mean pair interaction energy obtained after this simulation is equal to -30+/-5 kcal/mol. On the other hand, the HCQ molecule finds more rapidly its lower potential well (-40+/- 3kcal/mol) by diffusing in the solvent for 8ns then adsorbed on the protein in its lowest energy state from t=20 ns until the end of the simulations (Figure 8d-e). The final adsorption state of the HCQ is between the RBD and the NTD. The modification of the protein structure is of the same order as previously. However, we can notice that the RMSDs tend to the same stable value after 42ns of simulations for AZM and for HCQ (around 6Å) as shown in Figure 9b. The main deformation of the protein subjected to AZM is again at the bottom of the protein, close to the intrusion zone. Furthermore, AZM is directly responsible for this deformation since it adsorbs precisely in this zone at the end of the simulation and progressively unfold the protein. On the other hand, we observed that the bottom of protein remains stable during the interaction with HCQ while the RSMD exhibits the same feature. The S₂ zone (purple ribbon) of the protein is here impacted by the presence of the drug HCQ.

Figure 8: a-c) Initial and final AZM positions during the simulation (at t=0 ns, t=42 ns) and protein residues in interaction with AZM at the final stage of the simulation. d-f) Initial and final of HCQ positions during the simulation (at t=0 ns, and t=42 ns, respectively) and protein residues in interaction with HCQ at the final stage of the simulation. The RBD, NTD and S₂ parts of the protein were depicted in van der Waals, licorice and ribbon+lines modes for the three monomers shown in different colors.

    

    

    Figure 8: a-c) Initial and final AZM positions during the simulation (at t=0 ns,
t=42 ns) and protein residues in interaction with AZM at the final stage of the
simulation. d-f) Initial and final of HCQ positions during the simulation (at t=0
ns, and t=42 ns, respectively) and protein residues in interaction with HCQ at
the final stage of the simulation. The RBD, NTD and S2 parts of the protein
were depicted in van der Waals, licorice and ribbon+lines modes for the three
monomers shown in different colors.
    

Close

Figure 9: a) Pair interaction (smooth data in line representation), and b) RMSD of the 6VXX protein in presence of AZM (black) or HCQ (red) drugs.

    

    

    Figure 9: a) Pair interaction (smooth data in line representation), and b)
RMSD of the 6VXX protein in presence of AZM (black) or HCQ (red) drugs.
    

Close

Discussion

Recently, it was well shown, using combined docking and MD simulations, that both HCQ and AZM could block the binding of the protein virus to its cell receptor and more particularly, to gangliosides as soon as they are placed in the corresponding cavities [40]. AZM was found to occupy the binding domain of the spike protein and thus could limit virus binding to the cell receptor. HCQ was also responsible for blocking the virus binding through attachment to ganglioside and RBD, with the same mechanism as AZM but with a different targeted site [41]. It has also been pointed out that AZM and HCQ can act synergistically to improve therapy against viral disease. Their results thus confirm the strong binding of these pharmaceutical molecules against the viral protein [41,42].

Our MD were carried out simultaneously with first the influence of drugs (AZM or HCQ) on RBD linked to its cell receptor (ACE2). Our results, performed under physiological conditions, suggest that AZM is strongly docked to RBD when it is attached to ACE2 domain while HCQ is less bound to this protein conformation and even returned to the solvent. Our simulations showed that AZM could bind to the 6M0J while, for HCQ, diffusion to the stable adsorption site seemed more difficult to achieve. The modification of the protein is of the same order for the two different drugs. On the contrary, when the full S₁/S₂ parts of the protein are simulated in presence of AZM or HCQ, the deformation of the protein is greater regardless of the conditions of the simulations. Indeed, we observed a progressive unfolding of the lower part of the protein or of its S₂ domain. The pair interaction of HCQ with the RBD domain in the whole protein is of the same order as the AZM molecule with the bottom of the protein. In this particular case, the 6VXX presents a RBD domain in a closed state whereas one could estimate that, in the 6M0J case, the RBD is in its open state. The position and orientation of the RBD thus appear to be a fundamental means of understanding the key interaction of the HCQ and AZM with the protein.

Binding to virus receptor like SARS coronavirus (SARS-CoV), MERS coronavirus (MERS-CoV) and SARS-Cov-2 aims to recognize the receptor in the host cell. This RBD constantly evolves between two conformations, i. e. standing up (open state) and lying down (closed state). In the open state, the interaction between the S₁ and S₂ parts of viral proteins are weakened. These parts can thus easily dissociate and allow the fusion of the virus inside the host cell [19,21,43,44]. RBD thus plays a fundamental role in the attachment of the virus inside the host cell. NTD, located on the trimer side showed no conformational changes during the fusion process but may play a role in viral attachment.

Based on our results, we can suggest that AZM has two potential targets in SARS-Cov-2. First, it can adsorb at the bottom of the closed state of the viral protein (near the fusion zone) or, second, (nicely reported in literature [41,42]) be adsorbed close to the binding domain between the open state of RBD and the cell receptor. In contrast, we found only one sensitive adsorption site for HCQ in the closed state of the protein. Its interaction with two RBDs and one NTD could prevent the evolution of the viral protein from evolving from its closed to its open form. The literature concerning the use of HCQ in the Sars-COV2 disease shows some controversy over the efficacy of the molecules. [45,46] However, when HCQ is given early in the infection, the majority of clinical studies have reported positive effects [47,48]. Our results, although made only from numerical data, may explain why HCQ might be more effective in the early stages of the disease since our data suggest that its interaction is only effective on the closed state of viral protein. It may also explain why not all treatments using HCQ alone were as relevant as biotherapy [14,15]. Indeed, the AZM does not have the same action and could therefore serve to eradicate the virus before its fusion with the host cell, either by blocking the process of binding or by limiting the dissociation of the S₁ and S₂ part of the protein. This would help lower the viral load in the infected body and limite the degree of transmission of SARSCov- 2, as also demonstrated in the study by Fantini et al., [41]. In this latest study, it was suggested that the two drugs act together in specific targets of the SARS-Cov-2 protein, namely the RBD and cofactors of cell attachment. Taking at the earlier stage of the disease the double cocktail of HCQ+AZM would be an interesting solution to prevent the development of this virus.

Conclusion

Our molecular dynamics simulations aimed to study the interaction of hydroxychloroquine and azithromycin molecules with the SARS-Cov-2 protein. To do this, we first simulated the interaction of these molecules with the receptor binding domain and its host receptor in order to understand the role of each molecule on this important site. We demonstrated that AZM was strongly attracted to RBD in each simulation while HCQ had to be close to its adsorption site to fall inside. When the viral protein is in closed state, HCQ interacts with two RBDs in closed state and one NTD while AZM only adsorbes to the bottom of the viral protein, near the fusion domain. The deformations of the viral protein are always on the same order suggesting the equivalent role of drugs toward the virus. According to our data, HCQ can still interact with the closed state of the virus while AZM can act in both the closed and open states, explaining why biotherapy can be used to be effective.

Acknowledgement

Calculations were performed at the supercomputer regional facility Mesocentre of the University of Franche-Comté with the assistance of K. Mazouzi. This work was granted access to the HPC resources of IDRIS, Jean Zay supercomputer, under the allocation 2020 - DARI AP010711661 made by GENCI. We would like to express our gratitude to the IDRIS team (S Requena, PF Lavallée, R Lacroix and S Van Crienkingen), which was able to be very reactive to our request in a very tense pandemic climate, without whom this work would not have been possible.

References

1. Zaki SR, Goldsmith CS. SARS coronavirus infection: pathology and pathogenesis of an emerging virus disease, Coronaviruses with Special Emphasis on First Insights Concerning SARS. 2005; 87-99.
2. Drosten C, Günther S, Preiser W, van der Werf S, Brodt HR, Becker S, et al. Identification of a novel coronavirus in patients with severe acute respiratory syndrome, The New England journal of medicine. 2003; 348: 1967-1976.
3. Zhong NS, Zheng BJ, Li Poon YM, Xie ZH, Chan KH, Li PH, et al. Epidemiology and cause of Severe Acute Respiratory Syndrome (SARS) in Guangdong, People’s Republic of China, in February. 2003. Lancet (London, England). 2003; 362; 1353-1358.
4. Peiris JS, Lai ST, Poon LL, Guan Y, Yam LY, Lim W, et al. Coronavirus as a possible cause of severe acute respiratory syndrome, Lancet (London, England). 2003; 36: 1319-1325.
5. Zaki AM, van Boheemen S, Bestebroer TM, Osterhaus AD, Fouchier RA. Isolation of a novel coronavirus from a man with pneumonia in Saudi Arabia. The New England journal of medicine. 2012; 367; 1814-1820.
6. de Groot RJ, Baker SC, Baric RS, Brown CS, Drosten C, Enjuanes L, et al. Middle East respiratory syndrome coronavirus (MERS-CoV): announcement of the Coronavirus Study Group, J Virol. 2013; 87: 7790-7792.
7. van der Hoek L, Pyrc K, Jebbink MF, Vermeulen-Oost W, Berkhout RJ, Wolthers KC, et al. Identification of a new human coronavirus, Nature medicine. 2004; 10: 368-373.
8. Woo PC, Lau SK, Chu CM, Chan KH, Tsoi HW, Huang Y, et al. Characterization and complete genome sequence of a novel coronavirus, coronavirus HKU1, from patients with pneumonia. J Virol. 2005; 79; 884-895.
9. Woo PC, Lau SK, Yip CC, Huang Y, Yuen KY. More and More Coronaviruses: Human Coronavirus HKU1, Viruses. 2009; 1: 57-71.
10. Li W, Moore MJ, Vasilieva N, Sui J, Wong SK, Berne MA, et al. Angiotensinconverting enzyme 2 is a functional receptor for the SARS coronavirus, Nature. 2003; 426: 450-454.
11. Wrapp D, Wang N, Corbett KS, Goldsmith JA, Hsieh CL, Abiona O, et al. Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation, Science (New York, NY). 2020; 367: 1260-1263.
12. Zhou P, Yang XL, Wang XG, Hu B, Zhang L, Zhang W, et al. A pneumonia outbreak associated with a new coronavirus of probable bat origin, Nature. 2020; 579: 270-273.
13. Paraskevis D, Kostaki EG, Magiorkinis G, Panayiotakopoulos G, Sourvinos G, Tsiodras S. Full- genome evolutionary analysis of the novel corona virus (2019-nCoV) rejects the hypothesis of emergence as a result of a recent recombination event, Infection, genetics and evolution : journal of molecular epidemiology and evolutionary genetics in infectious diseases. 2020; 79: 104212.
14. Gautret P, Lagier JC, Parola P, Hoang VT, Meddeb L, Mailhe M, et al. Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial. International journal of antimicrobial agents. 2020; 382: 105949-105949.
15. Das S, Bhowmick S, Tiwari S, Sen S. An Updated Systematic Review of the Therapeutic Role of Hydroxychloroquine in Coronavirus Disease-19 (COVID-19), Clin Drug Investig. 2020; 1-11.
16. Erickson TB, Chai PR, Boyer EW. Chloroquine, hydroxychloroquine and COVID-19, Toxicol Commun. 2020; 4; 40-42.
17. Rodrigo C, Fernando SD, Rajapakse S. Clinical evidence for repurposing chloroquine and hydroxychloroquine as antiviral agents: a systematic review, Clin Microbiol Infect. 2020; S1198-1743X (1120): 30293-30297.
18. Baildya N, Ghosh NN, Chattopadhyay AP. Inhibitory activity of hydroxychloroquine on COVID-19 main protease: An insight from MDsimulation studies, Journal of Molecular Structure. 2020; 1219: 128595.
19. Gui M, Song W, Zhou H, Xu J, Chen S, Xiang Y, et al. Cryo-electron microscopy structures of the SARS-CoV spike glycoprotein reveal a prerequisite conformational state for receptor binding, Cell Research. 27; 2017: 119-129.
20. Walls AC, Tortorici MA, Frenz B, Snijder J, Li W, Rey FA, et al. Glycan shield and epitope masking of a coronavirus spike protein observed by cryo-electron microscopy, Nature structural & molecular biology. 2016; 23: 899-905.
21. Song W, Gui M, Wang X, Xiang Y. Cryo-EM structure of the SARS coronavirus spike glycoprotein in complex with its host cell receptor ACE2. 2018; 14: e1007236.
22. Li W, Moore MJ, Vasilieva N, Sui J, Wong SK, Berne MA, et al. Angiotensinconverting enzyme 2 is a functional receptor for the SARS coronavirus, Nature. 2003; 426: 450-454.
23. Zoete V, Cuendet MA, Grosdidier A, Michielin O. SwissParam: A fast force field generation tool for small organic molecules, Journal of Computational Chemistry. 2011; 32: 2359-2368.
24. Mayne CG, Saam J, Schulten K, Tajkhorshid E, Gumbart JC. Rapid parameterization of small molecules using the Force Field Toolkit, J Comput Chem. 2013; 34: 2757-2770.
25. Jo S, Cheng X, Islam SM, Huang L, Rui H, Zhu A, et al. CHARMM-GUI PDB Manipulator for Advanced Modeling and Simulations of Proteins Containing Nonstandard Residues. Karabencheva-Christova T, editors. In: Advances in Protein Chemistry and Structural Biology, Academic Press. 2014; 235-265.
26. Jo S, Kim T, Iyer VG, Im W. CHARMM-GUI: A web-based graphical user interface for CHARMM, Journal of Computational Chemistry. 2008; 29: 1859- 1865.
27. Lee J, Patel DS, Ståhle J, Park SJ, Kern NR, Kim S, et al. CHARMM-GUI Membrane Builder for Complex Biological Membrane Simulations with Glycolipids and Lipoglycans, J Chem Theory Comput. 2019; 15: 775-786.
28. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, et al. Schulten, Scalable molecular dynamics with NAMD, Journal of Computational Chemistry. 2005; 26: 1781-1802.
29. Best R, Zhu X, Shim J, Lopes P, Mittal J, Feig M, et al. Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone Φ, Ψ and Side- Chain Χ¹ and Χ² Dihedral Angles, Journal of Chemical Theory and Computation. 2012; 8; 3257-3273.
30. Jambrina PG, Aldegunde J. Computational Tools for the Study of Biomolecules. Martín M, Eden MR, Chemmangattuvalappil NG, editors. in: Computer Aided Chemical Engineering, Elsevier. 2016. 583-648.
31. Feller SE, Zhang Y, Pastor RW, Brooks BR. Constant pressure molecular dynamics simulation: The Langevin piston method. The Journal of Chemical Physics. 1995; 103: 4613-4621.
32. Darden T, York D, Pedersen L. Particle mesh Ewald: An Nlog (N) method for Ewald sums in large systems. The Journal of Chemical Physics. 1993; 98: 10089-10092.
33. Ryckaert J, Ciccotti G, Berendsen H. Numerical-Integration of Cartesian Equations of Motion of a System with Constraints - Molecular-Dynamics of N-Alkanes, Journal of Computational Physics. 1977; 23: 327-341.
34. Fung TS, Liu DX. Post-translational modifications of coronavirus proteins: roles and function, Future Virology. 2018; 13; 405-430.
35. Sugrue RJ. Viruses and Glycosylation: an overview. Sugrue RJ, editors. In: Glycovirology Protocols, Humana Press, Totowa, NJ. 2007; 1-13.
36. Goettig P. Effects of Glycosylation on the Enzymatic Activity and Mechanisms of Proteases, International Journal of Molecular Sciences. 2016; 17: 1969.
37. Zheng J, Yamada Y, Fung TS, Huang M, Chia R, Liu DX. Identification of N-linked glycosylation sites in the spike protein and their functional impact on the replication and infectivity of coronavirus infectious bronchitis virus in cell culture, Virology. 2018; 513: 65-74.
38. Böhm M, Bohne-Lang A, Frank M, Loss A, Rojas-Macias MA, et al. Glycosciences. DB: an annotated data collection linking glycomics and proteomics data (2018 update). Nucleic Acids Research. 2018; 47: D1195-D1201.
39. Hassan NM, Alhossary AA, Mu Y, Kwoh CK. Protein-Ligand Blind Docking Using QuickVina-W With Inter-Process Spatio-Temporal Integration. Scientific Reports. 2017; 7; 15451.
40. Matrosovich M, Herrler G, Klenk HD. Sialic Acid Receptors of Viruses. Gerardy-Schahn R, Delannoy P, von Itzstein M, editors. In: Top current Chemistry, Springer International Publishing. Cham. 2015; 1-28.
41. Fantini J, Chahinian H, Yahi N. Synergistic antiviral effect of hydroxychloroquine and azithromycin in combination against SARS-CoV-2: What molecular dynamics studies of virus-host interactions reveal. Int J Antimicrob Agents. 2020; 106020.
42. Fantini J, Di Scala C, Chahinian H, Yahi N. Structural and molecular modelling studies reveal a new mechanism of action of chloroquine and hydroxychloroquine against SARS-CoV-2 infection, International Journal of Antimicrobial Agents. 2020; 55: 105960.
43. Yuan Y, Cao D, Zhang Y, Ma J, Qi J, Wang Q, et al. Cryo- EM structures of MERS-CoV and SARS-CoV spike glycoproteins reveal the dynamic receptor binding domains. Nat Commun. 2017; 8: 15092.
44. Shang J, Wan Y, Liu C. Structure of mouse coronavirus spike protein complexed with receptor reveals mechanism for viral entry. PLOS Pathogens. 2020; 16: e1008392.
45. Alexander PE, Debono VB, Mammen MJ, Iorio A, Aryal K, Deng D, et al. COVID-19 coronavirus research has overall low methodological quality thus far: case in point for chloroquine/hydroxychloroquine, Journal of Clinical Epidemiology. 2020; 123: 120-126.
46. White NJ, Watson JA, Hoglund RM, Chan XHS, Cheah PY, Tarning J. COVID-19 prevention and treatment: A critical analysis of chloroquine and hydroxychloroquine clinical pharmacology, PLOS Medicine. 2020; 17: e1003252.
47. Sulaiman T, Mohana A, Alawdah L, Mahmoud N, Hassanein M, Wani T, et al. The Effect of Early Hydroxychloroquine-based Therapy in COVID-19 Patients in Ambulatory Care Settings: A Nationwide Prospective Cohort Study, medRxiv. 2020; 41-43.
48. Soto-Becerra P, Culquichicón C, Hurtado-Roca Y, Araujo-Castillo RV. Realworld effectiveness of hydroxychloroquine, azithromycin, and ivermectin among hospitalized COVID-19 patients: results of a target trial emulation using observational data from a nationwide healthcare system in Peru. medRxiv. 2020.

Download PDF

Citation: Picaud F and Herlem G. Hydroxychloroquine and Azithromycin Molecular Action against SARS-CoV-2 Viral Protein: A Molecular Dynamic Study. Austin J Nanomed Nanotechnol. 2021; 9(1): 1061.

Home

Journal Scope

Editorial Board

Instruction for Authors

Submit Your Article

Contact Us