Why Are Lopinavir and Ritonavir Effective against the Newly Emerged Coronavirus 2019? Atomistic Insights into the Inhibitory Mechanisms.

Since the emergence of a novel coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) was first reported from Wuhan, China, neither a specific vaccine nor an antiviral drug against SARS-CoV-2 has become available. However, a combination of two HIV-1 protease inhibitors, lopinavir and ritonavir, has been found to be effective against SARS-CoV, and both drugs could bind well to the SARS-CoV 3C-like protease (SARS-CoV 3CLpro). In this work, molecular complexation between each inhibitor and SARS-CoV-2 3CLpro was studied using all-atom molecular dynamics simulations, free energy calculations, and pair interaction energy analyses based on MM/PB(GB)SA and FMO-MP2/PCM/6-31G* methods. Both anti-HIV drugs interacted well with the residues at the active site of SARS-CoV-2 3CLpro. Ritonavir showed a somewhat higher number atomic contacts, a somewhat higher binding efficiency, and a somewhat higher number of key binding residues compared to lopinavir, which correspond with the slightly lower water accessibility at the 3CLpro active site. In addition, only ritonavir could interact with the oxyanion hole residues N142 and G143 via the formation of two hydrogen bonds. The interactions in terms of electrostatics, dispersion, and charge transfer played an important role in the drug binding. The obtained results demonstrated how repurposed anti-HIV drugs could be used to combat COVID-19.

In December 2019, there were many cases of patients reported to have a respiratory tract infection with severe pneumonia in Wuhan, China. It was found that these patients most likely had an epidemiological history related to a seafood market in that area of China.[1] However, a newly causative microbial infection cannot at first be identified in public databases. On the basis of whole genome sequencing, it was revealed that this microbial pathogen is a novel coronavirus, formally named 2019-nCoV, closely related to the bat severe acute respiratory syndrome (SARS)-like coronavirus, so-called SARS-CoV-2.[2,3] The World Health Organization (WHO) has officially confirmed the outbreak of 2019-nCoV on December 31, 2019, and eventually officially named it coronavirus disease 2019 or COVID-19. In general, coronaviruses are characterized as enveloped, positive-sense, single-stranded RNA viruses in the genus Coronavirus of the family Coronaviridae and can infect humans and several animals, including mammals and birds.[4−7] Nonetheless, some coronaviruses can potentially cause severe infection in patients such as the well-known outbreak of SARS-CoV in Guangdong, China,[8] and Middle East respiratory syndrome coronavirus (MERS-CoV) in many countries of the Middle East.[9] Likewise, COVID-19 has been confirmed to be transmitted from humans to humans and quickly spread in several countries throughout the world.[10]

SARS-CoV-2 is a betacoronavirus, like SARS-CoV and MERS-CoV, both of which have their origins in bats.[11] For the clinical symptoms, COVID-19 infection culminates in fatal pneumonia with the clinical presentation greatly resembling SARS-CoV infection.[1] Patients infected with SARS-CoV-2 might also develop acute respiratory distress syndrome, leading to a high rate of admission to intensive care units and ultimately death in severe cases.[7] After infection, patients presented mild to severe symptoms, including fever, cough, sore throat, rhinorrhea, severe pneumonia, and septic shock.[1,7]

To date, many companies and academic research groups around the world have focused on searching for and developing a specific vaccine or antiviral drug to prevent or control emerging SARS-CoV-2 infections (e.g., vaccine, monoclonal antibodies, and small-molecule drugs). However, these options need several months to years for their development. Because of the urgent need to alleviate the COVID-19 pandemic, the use of repurposed existing antiviral drugs approved for treatment of other viral infections such as human immunodeficiency virus (HIV), hepatitis B virus, hepatitis C virus, and influenza is somewhat promising,[12] based on previous successes of the therapeutic treatment with two relevant human coronaviruses, SARS-CoV and MERS-CoV. According to numerous previous studies,[1,7,13−16] the nonstructural protein of coronavirus, in particular, main proteases or 3C-like proteases (3CLpro), is considered an attractive drug target for the treatment of coronavirus infection. The role of this protease involves the proteolytic processing of the replicase polyprotein and is crucial for viral replication and maturation.[17] Moreover, 3CLpro has a similar common cleavage site among coronaviruses.[18] The sequence alignment of SARS-CoV-2 3CLpro (see Figure S1) shows that the SARS-CoV-2 proteinase is highly conserved compared to that of SARS-CoV with a 96.1% sequence identity.

A combination of the two approved drugs for HIV infection, lopinavir and ritonavir (KALETRA), has been reported to be active toward SARS and MERS.[14,19] Both anti-HIV drugs were initially purposed to inhibit 3CLpro of SARS-CoV and MERS-CoV, and they appeared to be related to clinical benefits of patients with SARS in a nonrandomized open-label trial.[17] Although ritonavir is a protease inhibitor, it is generally used to inhibit cytochrome P450 3A4 and markedly increases the plasma concentrations of other protease inhibitors.[20] Nevertheless, whether HIV protease inhibitors could effectively target SARS-CoV-2 3CLpro is under debate. This is based on the fact that HIV protease is from the aspartic protease family, whereas SARS-CoV-2 3CLpro belongs to the cysteine protease family. Previously, a theoretical study of the molecular interaction of lopinavir and ritonavir with 3CLpro of SARS-CoV suggested that these two drugs could bind well at the substrate-binding pocket of SARS-CoV 3CLpro.[15] To date, the three-dimensional structure of SARS-CoV-2 3CLpro in a complex with lopinavir and ritonavir has not been reported. Thus, in our study, we aimed to investigate the binding interactions of lopinavir and ritonavir with the SARS-CoV-2 proteinase using both molecular modeling and quantum chemical methods. It is our hope that this information can be useful for the future design or development of more specific inhibitors for the treatment of human coronaviruses.

Computational Details

System Preparation

The co-crystal structure of homodimeric SARS-CoV-2 3CLpro in complex with a peptidomimetic inhibitor [N-[(5-methylisoxazol-3-yl)carbonyl]alanyl-l-valyl-n-1-((1R,2Z)-4-(benzyloxy)-4-oxo-1-{[(3R)-2-oxopyrrolidin-3-yl]methyl}but-2-enyl)-l-leucinamide] was retrieved from the RSCB Protein Data Bank (PDB entry 6LU7).[21] Note that each moiety of the peptidomimetic inhibitor is associated with the -P5-P4-P3-P2-P1-P1′- positions of the 3CLpro cleavage site (Figure A). To generate the models of SARS-CoV-2 3CLpro with lopinavir and ritonavir bound to protomer A, the positions of the original peptidelike inhibitor were changed to the corresponding drug for each P position (Figure B,C). This modification was performed according to the similar orientation with the template inhibitor using the small molecules tool implemented in Accelrys Discovery Studio 2.5 (Accelrys Inc.). The protonation states of all ionizable amino acids were assigned at pH 7.4 using the H++ web server,[22] except for the catalytic residue H41 that was set as the neutral form with a protonated delta position (HID type of AMBER format) in accordance with the common reaction mechanism of cysteine protease.[23] The chemical structures of lopinavir (ZINC3951740) and ritonavir (ZINC3944422) were downloaded from the ZINC15 database.[24] The electrostatic potential (ESP) charges of both ligands were calculated with the HF/6-31G* level of theory using the Gaussian09 program.[25] The antechamber and parmchk modules of AMBER16 were used to generate the restrained ESP charges and missing parameters of the two drugs, respectively. The bonded and nonbonded parameters for the protein and ligand were treated with the AMBER ff14SB force field[26] and generalized AMBER force field version 2 (GAFF2),[27] accordingly. Missing hydrogen atoms were added using the LEaP module implemented in AMBER16. Afterward, the TIP3P water model[28] was used to solvate each system with a minimum distance of 10 Å between the protein surface and the solvation box edge, and the box dimensions were set to approximately 114 Å × 97 Å × 92 Å. The sodium ions were then randomly added to neutralize the simulated systems. The added hydrogen atoms and water molecules were minimized using 500 steps of steepest descent (SD) followed by 1500 steps of conjugated gradient (CG) methods before running the MD simulations, while the rest of the molecules were held fixed. Then, the protein and ligand were minimized by SD (500 iterations) and CG (1500 iterations) methods with constrained solvent molecules. Finally, the whole complex was fully minimized using the same procedure.

Figure 1

(A) Three-dimensional structure of the peptidelike inhibitor binding to the active site of the SARS-CoV-2 3CLpro homodimer (PDB entry 6LU7) in one asymmetric unit (A, yellow; B, cyan). Protomers are shown as ribbons, and the inhibitor is shown as an orange ball and stick model. Chemical structures of (B) lopinavir and (C) ritonavir, where the atomic labels are also given.

Molecular Dynamics Simulations

Each system was simulated under the periodic boundary condition with the isothermal–isobaric (NPT) scheme, as previously described.[29−32] Briefly, a cutoff distance for nonbonded interactions was set to 10 Å, while the particle mesh Ewald summation method[33] was employed to treat the electrostatic interactions. The SHAKE algorithm[34] was used to constrain all covalent bonds involving hydrogen atoms. A 2 fs simulation time step was used throughout the MD simulation. The temperature and pressure were controlled by the Langevin thermostat[35] with a collision frequency of 2 ps–1 and the Berendsen barostat[36] with a pressure-relaxation time of 1 ps. For the heating step, each simulated system was gradually heated from 10 to 310 K for 200 ps with a harmonic positional restraint of 30.0 kcal mol–1 Å–2 applied to the Cα atoms of the protein. To equilibrate the system, each complex was then subjected to four steps of restrained MD simulations at 310 K with harmonic restraints of 30, 20, 10, and 5 kcal mol–1 Å–2 for 1300 ps in total and another 200 ps without any restraint. Subsequently, the entire system was simulated under the NPT ensemble (310 K and 1 atm) until reaching 100 ns. The MD trajectories were saved every 10 ps. The calculations in terms of the root-mean-square displacement (RMSD), the number of protein–drug hydrogen bonds (# H-bonds), and the number of atomic contacts (# atom contacts) with the drug molecule were used to verify the stability of all investigated models as well as the structural features of each complex. The H-bond interactions were calculated using two criteria: (i) distance between the hydrogen donor (HD) and acceptor (HA) of ≤3.5 Å and (ii) HD–H···HA angle of ≥150°. The # atom contacts were counted from the number of 3CLpro atoms within 3.5 Å of each drug. To evaluate the binding affinity and the key binding residues involved in ligand binding of the protein–ligand complexes, the total binding free energy (ΔGbind) and per-residue decomposition free energy (ΔGbindresidue) calculations based on molecular mechanics/Poisson–Boltzmann (generalized Born) surface area [MM/PB(GB)SA] methods were performed on 100 MD snapshots extracted from the last 20 ns of the MD production phase. It should be noted that the protein–ligand interactions and the binding free energies were carried out using the CPPTRAJ[37] and MMPBSA.py[38] modules of AMBER16, respectively. Furthermore, the intermolecular interactions between the drug and its binding residues at the atomic level were described by pair interaction energy decomposition analysis (PIEDA) based on the fragment molecular orbital (FMO) method at the second-order Møller–Plesset perturbation energy and the 6-31G* level of theory (FMO-MP2/6-31G*) using GAMESS.[39] The solvation effect was also treated with a polarizable continuum model (PCM). The paired interaction energy (ΔEPIEDA) involved in ligand binding was evaluated by a summation of the electrostatic (EES), dispersion (EDI), charge transfer (ECT+mix), and charge exchange (EEX) energies, as well as PCM solvation (GsolPCM), using the following equation:[40−42]

Results and Discussion

System Stability

The stability of each simulated model was determined by calculating the all-atom RMSD for the drug–3CLpro complex, the number of intermolecular H-bonds (# H-bonds), and the number of atom contacts (# atom contacts) with the drug molecule versus the simulation time. As shown in Figure , the RMSD values of both lopinavir and ritonavir systems continuously increased in the first 40 ns and then maintained at a fluctuation of ∼2.0 Å until the end of the simulation time, as supported by a high fluctuation of # H-bonds and # atom contacts in the first 40 ns. It should be noted that the # H-bonds and # atom contacts of the ritonavir system (# H-bonds of 3 ± 1 and # atom contacts of 19 ± 4 over the last 20 ns) were higher than those of the lopinavir model (# H-bonds of 2 ± 1 and # atom contacts of 14 ± 4), suggesting that ritonavir was more susceptible to SARS-CoV-2 3CLpro (discussed in more detail below). In this work, the MD trajectories from 80 to 100 ns were thus extracted for further analysis in terms of (i) the binding affinity between lopinavir or ritonavir and SARS-CoV-2 3CLpro, (ii) key binding residues involved in drug binding, (iii) protein–drug H-bonding, and (iv) water accessibility at the enzyme active site.

Figure 2

All-atom RMSD, # H-bonds, and # atom contacts of lopinavir (left) and ritonavir (right) in complex with SARS-CoV-2 3CLpro plotted along the 100 ns MD simulation.

Predicted Inhibitory Efficiency

The susceptibility of lopinavir and ritonavir to SARS-CoV-2 3CLpro was estimated using the MM/PB(GB)SA approach on 100 snapshots extracted from the last 20 ns of simulation. Note that the MM/PBSA and MM/GBSA results in Table provided a similar trend in the binding free energy predictions.

Table 1

ΔGbind Values (kilocalories per mole) of Lopinavir and Ritonavir in Complex with SARS-CoV-2 3CLpro Calculated by the MM/PB(GB)SA Methoda

energy component	lopinavir	ritonavir
Gas Term
ΔEvdW	–50.29 ± 0.62	–74.09 ± 0.50
ΔEele	–19.68 ± 0.50	–23.35 ± 0.56
ΔEMM	–69.97 ± 0.87	–97.44 ± 0.60
–TΔS	22.47 ± 1.79	29.46 ± 1.77
Solvation Term
PBSA
ΔGsol(PBSA)ele	41.64 ± 0.53	59.47 ± 0.46
ΔGsol(PBSA)nonpolar	–5.03 ± 0.04	–6.42 ± 0.02
ΔGsol(PBSA)	36.61 ± 0.51	53.04 ± 0.45
GBSA
ΔGsol(GBSA)ele	40.11 ± 0.40	48.72 ± 0.37
ΔGsol(GBSA)nonpolar	–6.43 ± 0.07	–8.53 ± 0.04
ΔGsol(GBSA)	33.67 ± 0.37	40.19 ± 0.37
Binding Free Energy
ΔGbind(MM/PBSA)	–10.89 ± 1.89	–14.93 ± 1.83
ΔGbind(MM/GBSA)	–13.83 ± 1.91	–27.78 ± 1.82

Data are shown as means ± the standard error of the mean (SEM).

According to the molecular mechanics (ΔEMM) calculations, we found that van der Waals (vdW) interaction is the main force inducing molecular complexation with SARS-CoV-2 3CLpro of both lopinavir (ΔEvdW of −50.29 ± 0.62 kcal/mol) and ritonavir (−74.09 ± 0.50 kcal/mol) and is ∼2–3-fold stronger than the electrostatic attraction (ΔEele values of −19.68 ± 0.50 and −23.35 ± 0.56 kcal/mol, respectively). This finding was in a good agreement with the ΔEMM results of (i) these two anti-HIV drugs in complex with SARS-CoV 3CLpro[15] and (ii) darunavir and amprenavir binding to the HIV-1 protease.[43] Taken together with the solvation effect and entropic term, the binding affinity (ΔGbind) of ritonavir with SARS-CoV-2 3CLpro (−14.93 ± 1.83 and −27.78 ± 1.82 kcal/mol taken from the MM/PBSA and MM/GBSA methods, respectively) was relatively higher than that of another drug (−10.89 ± 1.89 and −13.83 ± 1.91 kcal/mol, respectively). The predicted results suggested that ritonavir may have an ∼1.4–2.0-fold greater inhibitory efficiency than lopinavir on the focused target; however, it should be noted that this drug can clinically also boost the lopinavir efficacy in the fight against COVID-19.

Key Residues for the Repurposed HIV Drug to Combat COVID-19

To evaluate the key residues involved in anti-HIV drug binding at the active site of SARS-CoV-2 3CLpro, the ΔGbindresidue calculation based on the MM/GBSA method was performed. The total energy contribution from each residue associated with drug binding of both complexes is plotted in Figure , and the drug orientation in the enzyme active site is illustrated in the right panel, in which the contributing amino acids are colored according to their ΔGbindresidue values. It is important to note that only residues exhibiting an energy stabilization of ≤−1.0 kcal/mol were taken into account.

Figure 3

ΔGbindresidue values of lopinavir (top) and ritonavir (bottom) in complex with SARS-CoV-2 3CLpro. The contributing residues involved in ligand binding are colored according to their ΔGbindresidue values, where the highest to lowest free energies are shaded from gray to blue, respectively.

The obtained results revealed that there were four (M49, M165, P168, and Q189) and nine residues (L27, H41, M49, F140, N142, G143, H164, M165, and E166) that were important for the binding of lopinavir and ritonavir, respectively. The more strongly contributing residues found in the SARS-CoV-2 3CLpro–ritonavir complex were in accordance with the results mentioned above (Figure and Table ). Similarly, it was reported that the SARS-CoV 3CLpro residues (i) M49 and M165 and (ii) E166 strongly stabilized PJ207, a quinazoline-containing compound, mainly through hydrophobic interaction and H-bond formation, respectively.[44] In line with this finding, our ΔGbindresidue calculation showed that the phenyl ring at the P1 position of lopinavir and the thiazole ring at the P2 position of ritonavir were inserted into the S2 pocket of SARS-CoV-2 3CLpro, forming hydrophobic contacts with residues M49 and M165 (Figure ), and the O1 atom of ritonavir was stabilized by the -NH group of E166 via H-bond formation (Figure , discussed below). It is of note that only ritonavir could interact with the catalytic residue H41[45] and the oxyanion hole residues N142 and G143[46] (ΔGbindresidue of −1.4 kcal/mol for each residue), suggesting that the molecular structure of ritonavir fitted well within the active site of SARS-CoV-2 3CLpro. The obtained information was supported by the well-aligned structures of ritonavir over the last 20 ns simulation (Figure ), whereas a large fluctuation of the P2′ site was detected in lopinavir.

Figure 6

Pair interaction energy decomposition analysis (PIEDA) of lopinavir (top) and ritonavir (bottom) interacting with individual residues in the binding pocket of SARS-CoV-2 3CLpro based on FMO-MP2/PCM/6-31G* calculation.

Figure 4

Superimposed structures over the 20 snapshots of (A) lopinavir and (B) ritonavir in complex with SARS-CoV-2 3CLpro derived from the last 20 ns of MD simulations.

In terms of the contribution from the electrostatic (ΔEele + ΔGpolar, black line) and vdW (ΔEvdW + ΔGnonpolar, red line) energies from each residue, it can be seen from Figure that the main energy contribution for stabilizing the HIV-1 protease inhibitors lopinavir and ritonavir was the vdW energy (≲−4.0 kcal/mol), especially for residues M49, M165, and Q189, whereas the electrostatic contribution was observed in the range of approximately −1.0 to 2.0 kcal/mol as related to the ΔEMM calculations (see Table ). However, there is a common concern about the drawbacks of the molecular mechanics calculations, in particular electronic properties. The FMO-MP2/PCM/6-31G* calculation was additionally performed on the last snapshot of each complex to investigate the protein–drug interactions at the enzyme active site. The advantage of this method is known to include the electron correlation, which trustworthily describes the hydrophobic effect and hydrogen bond interaction between the ligand and individual residue.[40,47,48]Figure highlights that residue Q189 dramatically stabilized lopinavir (−37.3 kcal/mol) mainly through electrostatic, dispersion, and charge transfer contributions (−28.1, −18.3, and −8.1 kcal/mol, respectively), whereas residues H41, A46, M49, E166, L167, L187, A191, and A193 had moderately stabilizing effects with energy values in the range of approximately −5 to −10 kcal/mol via dispersion and/or electrostatic interactions. In contrast, there were four essential amino acids (i.e., N142, M165, E166, and Q189) contributing to ritonavir binding with an energy stabilization of approximately −12 to −20 kcal/mol, whereas residues G143, S144, C145, and D187 showed less energy stabilization of approximately −5 to −7 kcal/mol, most likely through dispersion and electrostatic interactions, except for residue E166, for which the impact of the solvation effect was predominantly found on the two drugs binding to SARS-CoV-2 3CLpro.

Figure 5

Electrostatic (ΔEele + ΔGpolar, black line) and vdW (ΔEvdW + ΔGnonpolar, red line) energy contributions from each residue of SARS-CoV-2 3CLpro to the binding of lopinavir (top) and ritonavir (bottom).

As compared to the recently reported potent α-ketoamide inhibitor (13b) bound to SARS-CoV-2 3CLpro,[49] the anti-HIV drugs lopinavir and ritonavir lack (i) the P3 moiety, (ii) an α-ketoamide group located between positions P1 and P1′, and (iii) polar moieties on the phenyl ring in the P1 region. However, the phenyl ring at the P1′ site of both HIV-1 protease inhibitors and 13b remains the same, and the P2 (larger and more hydrophilic than 13b) and P2′ (found only in HIV-1 protease inhibitors) regions of lopinavir and ritonavir strongly interacted with SARS-CoV-2 3CLpro at the S2 and S2′ pockets, respectively (Figure ). Accordingly, on the basis of this evidence, we proposed a rational design of novel protease inhibitor(s) derived from lopinavir and ritonavir to enhance the binding efficiency with SARS-CoV-2 3CLpro as follows: (i) changing the phenyl ring in the P1 region to the γ-lactam ring to increase the extent of H-bond formation with residues F140, H163, and E166, (ii) introducing P3 polar moieties (e.g., pyridone ring) to form H-bonds with the negatively charged E166, and (iii) adding the reactive α-ketoamide group located between positions P1 and P1′ to covalently bond with the catalytic residue C145, as depicted in Figure .

Figure 7

Rational drug design of the SARS-CoV-2 3CLpro inhibitors. Note that the green dashed line indicates H-bond formation.

Protein–Drug Hydrogen Bonding

Although the protein–drug interactions were primarily contributed by the vdW interactions, the formation of H-bonds between both anti-HIV drugs and their surrounding residues in the active site of SARS-CoV-2 3CLpro could be an important factor in the inhibition of this targeted enzyme. To monitor such interaction, the percentages of H-bond occupations are plotted in Figure , whereas the representative structures determined from the last MD snapshot are depicted in Figure . As expected, H-bond formation (three to four H-bonds) between the drugs and the surrounding residues inside the active site of SARS-CoV-2 3CLpro was observed. In the case of lopinavir, there were three H-bonds: (i) O1(P2 site)···H–NE2(Q189) at 70.6%, (ii) O(T190)···H–N1(P2 site) at 14.8%, and (iii) OE1(Q189)···H–N4(P1 site) at 83.9%. Meanwhile, four H-bonds were detected in the ritonavir complex: (i) O1(P2 site)···H–N(E166) at 65.6%, (ii) O(H164)···H–N4(P1 site) at 72.6%, (iii) OD1(N142)···H–N5(P1′ site) at 17.3%, and (iv) O3(P1′ site)···H–N(G143) at 85.6%. It can be noticed that the importance of residues E166 and Q189 also agreed well with the reported binding of such HIV-1 protease inhibitors to SARS-CoV.[15] Notably, ritonavir showed a slightly higher level of H-bond formation with 3CLpro, and its P1′ site was stabilized by the oxyanion hole residues N142 and G143.

Figure 8

Percentage of H-bond occupation of SARS-CoV-2 3CLpro that contributed to the binding of lopinavir (top) and ritonavir (bottom).

Figure 9

Binding patterns of (A) lopinavir and (B) ritonavir in complex with SARS-CoV-2 3CLpro demonstrated from the last MD snapshot. Black dashed lines represent H-bonds.

Solvent Accessibility in the 3CLpro Active Site

As the active site of SARS-CoV-2 3CLpro is located near the solvent-exposed area, water molecules might play a role in stabilizing protein–ligand interactions. To characterize the effect of water accessibility at the 3CLpro active site, the solvent-accessible surface area (SASA) (Figure ) calculations were performed on the residues within 4 Å of each anti-HIV drug. Moreover, the radial distribution function (RDF) (Figure ) of water oxygens toward all heteroatoms of lopinavir and ritonavir was used to verify the number of water molecules approaching the ligands.

Figure 10

(A) SARS-CoV-2 3CLpro homodimer, in which chain A with a drug bound and chain B without a drug bound are shown in shades of yellow and blue, respectively. Note that the amino acid residues within 4 Å (stick model) of the ligand (ball and stick representation) were used for SASA calculations. (B) SASA plots along the simulation time of the two studied systems. (C) Average SASAs of lopinavir and ritonavir systems.

Figure 11

Radial distribution functions, g(r), of the water oxygen atom and integration numbers, n(r), up to the first minimum around the heteroatoms of (A) lopinavir and (B) ritonavir (see Figure for definitions) in complex with SARS-CoV-2 3CLpro. H-Bond water networks of the drug, bridging water, and SARS-CoV-2 3CLpro residue are shown as black dashed lines.

Our complex model of the SARS-CoV-2 3CLpro homodimer contained the anti-HIV drugs binding only to protomer A. It can be seen from Figure that the SASAs for protomer B (blue) were 883.81 ± 75.24 and 918.15 ± 110.03 Å2 for the lopinavir and ritonavir systems, respectively. Upon molecular complexation in protomer A (yellow), the SASAs of both lopinavir (500.21 ± 67.38 Å2) and ritonavir (451.60 ± 61.62 Å2) systems dramatically decreased in a manner similar to those of other studies.[50−52] By considering the RDF results, we found that the first peak at ∼3 Å of the lopinavir system (up to ∼2) was slightly higher than that of the ritonavir system (up to ∼1.4), indicating that accessible water molecules were more pronounced in the SARS-CoV-2 3CLpro–lopinavir complex. As the H-bond water network, protein···water···drug, also plays a role in stabilizing the protein–ligand complex,[53] we further investigated the H-bond water bridges around the heteroatoms of both ligands using the same structural criteria of H-bonds described elsewhere. The obtained results demonstrated that there was only one bridging water coordinating with (i) lopinavir (O2) and E166 as well as (ii) ritonavir (O4) and N142 with a very low percentage of occupancy (∼5–13%), suggesting that the bridging water did not make a significant contribution to mediating the H-bond network of protein–ligand binding.

Conclusions

In this work, the binding pattern and susceptibility of the two HIV-1 protease inhibitors lopinavir and ritonavir in complex with SARS-CoV-2 3CLpro were fully revealed by all-atom MD simulations, binding free energy estimation, and PIEDA based on the MM/PB(GB)SA and FMO-MP2/PCM/6-31G* calculations, respectively. According to ΔGbind prediction, the susceptibility against SARS-CoV-2 3CLpro of ritonavir was somewhat higher than that of lopinavir, supported by energy stabilization from individual residues that resulted from both methods: (i) M49, M165, P168, and Q189 from MM/GBSA for lopinavir and L27, H41, M49, F140, N142, G143, H164, M165, and E166 from MM/GBSA for ritonavir and (ii) H41, A46, M49, E166, L167, L187, Q189, A191, and A193 from FMO-MP2/PCM/6-31G* for lopinavir and N142, G143, S144, C145, M165, E166, D187, and Q189 from FMO-MP2/PCM/6-31G* for ritonavir. In addition, the oxyanion hole residues N142 and G143 were found to interact with ritonavir via hydrogen bonds. From the FMO-MP2/PCM/6-31G* data, the electrostatics, dispersion, and charge transfer were considered as the important interactions for drug binding. The obtained results demonstrated how repurposed anti-HIV drugs could be used to combat COVID-19 and how fundamental knowledge at the atomic level could also be helpful for the further design or development of more specific inhibitors in treating human coronaviruses.