Interaction between favipiravir and hydroxychloroquine and their combined drug assessment: in silico investigations.

Hydroxychloroquine (HCQ) and favipiravir (FPV) are known to be effective antivirals, and there are reports about their use to fight the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) despite that these are not conclusive. The use of combined drugs is common in drug discovery, and thus, we investigated HCQ and FPV as a combined drug. The density functional theory method was used for the optimization of geometries, spectroscopic analysis and calculation of reactivity parameters. The quantum theory of atoms in molecules was applied to explain the nature of the hydrogen bonds and confirm the higher stability of the combined drug. We also evaluated the absorption, distribution, metabolism and excretion (ADME) parameters to assess their drug actions jointly using SwissADME. The preliminary findings of our theoretical study are promising for further investigations of more potent and selective antiviral drugs. Supplementary Information: The online version contains supplementary material available at 10.1007/s11696-021-01946-8. © Institute of Chemistry, Slovak Academy of Sciences 2021.

Introduction

In late 2019, a new type of coronavirus named SARS-CoV-2 appeared in Wuhan, China, and rapidly spread all over the World (Lai et al. 2020). To date, there have already been 164,929,440 persons being infected worldwide and about 3,419,966 persons lost their lives (COVID-19: Coronavirus Pandemic, Last updated 2021). One of the ways to combat the spread of this highly contagious virus has been to adopt confinement (Zhang et al. 2020), and this has led to huge economic lost (McKibbin and Fernando 2020) among several issues (Wong et al. 2020). Thus, there is an urgent need to identify suitable and reliable drugs (Le et al. 2020). One of the different strategies is to rely on existing drugs as the safety profiles of these drugs are already established (Colson et al. 2019; Gautret et al. 2020). Another strategy is to use combinations of drugs so as to improve the pharmacodynamics and pharmacokinetic efficacy (Flynn 2007). Although vaccination appears to be the most effective way to prevent an outbreak for a particular infectious disease, the development and the assessment of a vaccine require a lot of time. Nevertheless, there are many vaccines are available that have proven to be successful such as influenza, HPV and chickenpox vaccines (Fiore et al. 2009; Chang et al. 2009; Liesegang 2009).

FPV also known as T-705, as shown in Scheme 1, is an antiviral agent, and it is also known for its broad spectral activity against ribonucleic acid (RNA) viruses (Furuta et al. 2009). This drug was approved by Shenzan Health Commission for the possible treatment of COVID-19 patients (Wu et al. 2020). HCQ is another drug (Scheme 1) which has been reported for the possible treatment of COVID-19 patients (Gautret et al. 2020; Singh et al. 2020a). There are reports indicating that FPV and HCQ have some positive effects on the treatment of less severe COVID-19 patients (Singh et al. 2020b).

Scheme 1

Molecular structure of HCQ and FPV

QTAIM analysis has been a useful tool for the examination of intra- and intermolecular hydrogen bonding (HB) interactions (Bader 1990; Savin et al. 1997) although several objections were raised about the terminology preferred (Shahbazian 2018; Jabłoński 2019). At the base of this theory, three parameters of bond critical points (BCP), electron density ρ(r), Laplacian of electron density ∇2ρ(r) and electronic energy density H(r) are important to assess the nature of interaction (Bader 1990). According to Rozas et al. (Rozas et al. 2000), the nature of HB is defined as strong if ∇2ρ(r) and H(r) are negative. If ∇2ρ(r) is positive and H(r) is negative, then the HB interaction is intermediate or partially covalent. Finally, if ∇2ρ(r) and H(r) are positive, then the HB interactions are called as weak or non-covalent (Rozas et al. 2000).

The ability of electron donation and electron acceptance for a given molecular system can be assessed by analyzing highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO). The HOMO–LUMO energy gap is a well-known parameter to evaluate molecular charge transfers processes. They are also used for the identification of chemical reactivity, bioactivity, the stability and general reactivity of molecular systems (Fukui 1982; Parr et al. 1999; Samsonowicz et al. 2017).

Density functional theory (DFT) method has been widely applied for the examination of the structural, electronic and spectroscopic properties of different compounds (Hao et al. 2019; Astani et al. 2016). DFT is a useful tool for modeling drug design and the examination of molecular interactions since it is cost-effective and faster than the time-consuming experimental methods (Astani et al. 2016; Pierens et al. 2017; Khorram et al. 2017). In addition to ab initio method, the SwissADME web tool (Daina et al. 2017) is commonly used to have drug-related parameters.

In view of the above resume and in continuation with our interests to study drugs using theoretical methods (Rhyman et al. 2018; Alver et al. 2019; Parlak and Alver 2017), we were interested to investigate FPV and HCQ working jointly as a drug. Within this aim, we set objectives, namely: (i) to use DFT and QTAIM to understand the nature of the interactions between the two drugs and (ii) to use SwissADME to have insights in the drug capabilities of the combined drugs. We hereby report the findings of our research.

Computational details

All ab initio calculations were performed using Gaussian 09 (Frisch et al. 2009) within the DFT framework with the B3LYP functional in conjunction with the 6-31G(d) basis set for all atoms. The polarizable continuum model was used to include the solvation effect (Tomasi et al. 2005). Molecular structures and initial geometrical configurations were built and visualized with the GaussView program (Dennington et al. 2009). Firstly, isolated FPV and HCQ were optimized in water without imposing any structural limitations. Geometry optimizations were followed by vibrational frequency calculations to ensure that the resultant structures converge to a ground state or a minimum energy state. An analysis of the lowest energy configurations natural bond orbital (NBO) charge (Reed et al. 1988) was conducted to identify the possible interaction edges between FPV and HCQ. After having decided the possible interaction edges for FPV and HCQ, the combined drugs were optimized until no imaginary frequency was obtained.

We followed two steps to determine the most possible geometry of the interacted drugs. Initially, binding energies (Eb) of the possible drug couples are calculated using Eq. (1) (Zhai et al. 2019):

In Eq. (1), Ecombined, EFPV and EHCQ are the optimized electronic energies of the investigated systems. It is known that a negative value of Eb reflects a favorable process. Consequently, a positive value of Eb means an unfavorable process. After the optimizations, the structures with positive Eb values were not considered for further investigations. It is also known that at the sites where the interactions take place, basis functions overlap and this leads to an error known as basis set superposition error (BSSE) causing some alterations in the Eb energies of the investigated molecular systems (Simon et al. 1996; Zheng et al. 2016). Thus, in the present study, the BSSE errors were included by using the counterpoise correction method (Boys and Bernardi 1970). After eliminating the unlikely candidates of drug couples, a topological analysis of electron density was carried out with QTAIM analyses for the structures having negative Eb values. QTAIM calculations were performed by the Multiwfn program (Lu and Chen 2012). SwissADME web tool was used to predict the drug-related parameters (SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. 2017).

Results and discussion

The coordinates of the optimized structures and the dipole moments FPV, HCQ and interacted drugs (C1-C5) are collected in Tables SI1-7 (SI: Supplementary information). The optimized structures (five possibilities) of the interacted drugs are displayed in Fig. 1 with the Eb energy. The calculated Eb values are found to be positive for C1, C4 and C5 and thus unfavorable. Therefore, in this section, the discussion is limited to C2 and C3 structures since they have exothermic binding energies values of -4.45 kcal/mol and -3.41 kcal/mol.

Fig. 1

Optimized structures with Eb energies

Analysis of C2 interaction

The molecular topography map and some important topological parameters of C2 are given in Fig. 2. In this representation, BCPs are shown by orange spheres. It can be observed from Fig. 2 that QTAIM calculations provide two possible intermolecular HB interactions, namely site I and site II. At the site I, interaction occurs via O─H…O group. For this edge, the ∇2ρ(r) was calculated as positive and the value of H(r) was found to be negative, and thus, these indicate the partial covalent interaction with EHB value of -8.02 kcal/mol. At site II, interaction occurs via N─H…O groups with EHB value of -6.37 kcal/mol. It is seen that the HB character at site I is stronger than site II leading to a larger electron density at BCP.

Fig. 2

Molecular topography map of C2

In C2, the effect of interaction can be best followed from the hydroxyl edge of HCQ drug since it is easy to identify the characteristic OH vibration by the infrared (IR) spectroscopic method. For a single isolated HCQ, the OH stretching vibration is observed at 3737 cm−1 with an IR intensity of 26.8 km/mol (Fig. 3). After interaction with FPV, this band shifts down to 3541 cm−1 and IR intensity increases to 747.2 km/mol. This red-shift indicates a weakening in the bond character of OH. As a result, the force constant of OH bond of HCQ decreases from 8.77 to 7.89 mdyn/Å upon interaction with FPV. The reduction in the force constant causes a change in the OH bond length of HCQ from 0.971 Å to 0.982 Å after interaction with FPV (Fig. 3).

Fig. 3

Vibrational properties and bond characteristics at the interaction site of C2

HOMO–LUMO distributions over the surface of a given molecular structure provide useful information for the identification of molecular activity. For example, the more negative LUMO energy values are considered as the descriptors of chemically more active molecules (Fukui 1982; Ayala and Scuseria 1999; Kawakami et al. 2013). It was previously reported that molecular compounds with more stabilized LUMO orbitals have higher biological activities (Kumar et al. 2018). It was also reported that the molecular sites where the LUMO orbitals are distributed play an important role for the biological activities (Kawakami et al. 2013; Kumar et al. 2018). Furthermore, ionization potential (I) estimated as -EHOMO and electron affinity (A) estimated as -ELUMO are also considered as general reactivity identifiers (Samsonowicz et al. 2018). HOMO–LUMO plots and calculated energy values for C2 are given in Fig. 4. The HOMO orbital is located on HCQ, and the LUMO orbital is distributed over FPV. These results indicate that while HCQ behaves as an electron donor, FPV behaves as an electron acceptor. ELUMO values of single isolated FPV and HCQ were calculated as -2.27 eV and -1.26 eV. After the interaction, ELUMO of C2 becomes more negative with a value of -2.33 eV suggesting an increase in the chemical activity and the electron affinity of the C2 system compared to FPV and HCQ separately. There is also an increase in the activity of FPV for the interacted C2 system since LUMO is located over the FPV following the interaction. Eg value of C2 (3.30 eV) is found to be smaller than FPV (4.57 eV) and HCQ (4.37 eV), and this also shows that the chemical reactivity of C2 is higher than FPV and HCQ alone.

Fig. 4

Frontier molecular orbitals of C2

Analysis of C3 interaction

The molecular topography map along with some important topological parameters of C3 is provided in Fig. 5. The QTAIM calculations for C3 indicate two possible intermolecular HB interactions labeled as site I and site II. At site I, interaction occurs via N─H…O groups. For this edge, it is found that ∇2ρ(r) > 0 and H(r) < 0 and these refer to a partial covalent interaction with EHB value of -8.04 kcal/mol. At site II (O─H…N), both ∇2ρ(r) and H(r) are found to be positive indicating a non-covalent interaction with EHB value of -4.92 kcal/mol. Electron density at BCP was found larger for site I interaction (Fig. 5).

Fig. 5

Molecular topography map of C3

In C3, the OH stretching vibration of HCQ is red shifted to 3658 cm−1 as a result of HB interaction and IR intensity increases to 348 km/mol (Fig. 6). The force constant of OH bond of HCQ decreases to 8.43 mdyn/Å upon interaction with FPV. Consequently, the weakening of the force constant causes an increase in the OH bond length of HCQ from 0.971 to 0.976 Å following the interaction with PFV. The HOMO orbital is located on HCQ and the LUMO orbital is in FPV. The ELUMO value of C3 is found to be more negative (-2.34 eV) compared to FPV (-2.27 eV) and HCQ (-1.26 eV) separately. As shown from Fig. 7, Eg value of C3 (3.29 eV) is also found to be smaller than drugs (4.57 eV and 4.37 eV). Eg and ELUMO values are found to be slightly smaller than C2. These indicate an increase in activity following the interaction.

Fig. 6

Vibrational properties and bond characteristics at the interaction site of C3

Fig. 7

Frontier molecular orbitals of C3

SwissADME analysis

We predicted the physicochemical properties and drug-likeness of FPV, HCQ and their combined structure. Physicochemical properties give a global description of the structure of the compound including molecular weight (g/mol), molecular refractivity, topological polar surface area (Å2), number of rotatable bonds, heavy atoms and hydrogen bond acceptors and donors (Table 1).

Table 1

Physicochemical properties and lipophilicity of FPV, HCQ and FPV–HCQ

Properties	FPV	HCQ	FPV–HCQ
Molecular weight	157.10	335.87	492.97
Heavy atom	11	23	34
Arom. heavy atom	0	0	0
Fraction Csp3	0.00	0.56	0.43
Rotatable bond	1	9	10
H–Bond acceptor	5	3	8
H–Bond donor	2	2	4
Molar refractivity	45.10	102.41	147.51
Polar surface area	63.32	35.50	98.82
Lipophilicity
MLOGP	–1.84	1.59	–0.39
WLOGP	–2.89	1.42	–1.46
XLOGP3	–2.89	0.70	–2.20

The bioavailability radar of the compounds displayed a rapid evaluation of drug-likeness. It includes the following six physicochemical properties: (1) Lipophilicity (XLOGP3 between − 0.7 and + 5.0), (2) size (molecular weight between 150 and 500 g/mol), (3) polarity (the total polar surface area between 20 and 130 Å2), (4) solubility (log S not higher than 6), (5) saturation (fraction Csp3 not less than 0.25) and (6) flexibility (the number of rotatable bonds not more than 9).

The bioavailability radar of FPV, HCQ and upon interaction for six physicochemical properties is shown in Fig. 8. The pink area represents the optimal range of these properties, and the red line represents the properties of the compounds. In Fig. 8, the red lines of the combined compound are in the range of the pink area except for the flexibility value (slightly). Therefore, we can conclude that the combined drugs are predicted orally as bioavailable.

Fig. 8

Bioavailability radar of FPV (a), HCQ (b) and (c) FPV–HCQ

Furthermore, drug-likeness was established based on the physicochemical properties to find oral drug candidates. There are five different rule-based filters which are defined as follows:

(1) Lipinski’s filter includes molecular weight ≤ 500, MLOGP (lipophilicity) ≤ 4.15, hydrogen bond acceptors ≤ 10 and hydrogen bond donors ≤ 5.

(2) Ghose’s filter includes 160 ≤ molecular weight ≤ 480, − 0.4 ≤ WLOGP (lipophilicity) ≤ 5.6, 40 ≤ the molar refractivity ≤ 130 and 20 ≤ number of atoms ≤ 70.

(3) Veber’s filter includes the number of rotatable bonds ≤ 10 and the total polar surface area ≤ 140.

(4) Egan’s filter includes WLOGP (Lipophilicity) ≤ 5.88 and the total polar surface area ≤ 131.6.

(5) Muegge’s filter includes 200 ≤ molecular weight ≤ 600, − 2 ≤ XLOGP3 (lipophilicity) ≤ 5, the total polar surface area ≤ 150, the number of rings ≤ 7, the number of carbon > 4, the number of heteroatoms > 1, the number of rotatable bonds ≤ 15, the hydrogen bond acceptors ≤ 10 and the hydrogen bond donors ≤ 5.

The result of the drug-likeness evaluation of FPV, HCQ and upon interaction is shown in Table 2. The combined drug is in agreement with the Lipinski’s, Veber’s and Egan’s rules. Turning to Ghose’s four rules, it is in agreement with one rule excluding MW > 480, MR > 130 and WLOGP < -0.4. The combined drug passes to Muegge’s nine rules excluding one lipophilicity value (XLOGP3 < -2). It can be concluded that this interaction is predicted orally bioavailable. These preliminary results provide the lead for the design of more potent and selective antiviral drugs.

Table 2

Drug-likeness evaluation of FPV, HCQ and FPV–HCQ

Rule-based filters	FPV	HCQ	FPV–HCQ
Lipinski	Yes	Yes	Yes
Ghose	No; 3 violations	Yes	No; 3 violations
Veber	Yes	Yes	Yes
Egan	Yes	Yes	Yes
Muegge	No; 2 violation	Yes	No; 1 violations

Conclusions

The interaction mechanism, structural, reactivity and spectroscopic examinations of the combined FPV and HCQ drugs were investigated with DFT and QTAIM approaches. The drug-likeness evaluations of FPV–HCQ were also performed by ADME parameters.

It was observed that following the interaction LUMOs were located over the FPV and ELUMO values of the interacted systems became more negative indicating an increase in the antiviral activity and electron affinity. The charge transfer was found from HCQ to FPV. It is shown that for the drug interactions occurring through electron acceptance, low LUMO energy values suggest that FPV–HCQ couple shows higher activities via LUMOs distributing over FPV when compared to single FPV. It is also observed that HB interactions with bond energies almost -8 kcal/mol lead to a red-shift for the hydroxyl stretching vibration of HCQ. Further, the IR intensity and bond length of the hydroxyl group increase as a result of HB interactions.

C3 interaction has slightly more negative LUMO energy and smaller gap energy than C2, whereas C2 is more energetically stable than C3. Overall, both structures are effective, suitable for this combination and increase the activity of FPV. Further, FPV–HCQ interacted system fulfills all requirements of Lipinski, Veber and Egan rules, and it is predicted to be orally bioavailable.

These preliminary results may provide the lead for the design of more potent and selective COVID-19 drugs or vaccines. It is also worth noting that the possible combination of FPV and HCQ has the capacity to increase the effectiveness based on electron acceptance properties. Therefore, the combination might open new considerations and strategies against the COVID-19 and possible challenges in the future.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 28 kb)