ORIGINAL RESEARCH article
Sec. Theoretical and Computational Chemistry
Quantum-Chemistry Based Design of Halobenzene Derivatives With Augmented Affinities for the HIV-1 Viral G4/C16 Base-Pair
- 1Sorbonne Université, Laboratoire de Chimie Théorique, UMR7616 CNRS, Paris, France
- 2UR EGP, Centre d'Analyses et de Recherche, Faculté des Sciences, Université Saint-Joseph de Beyrouth, Beirut, Lebanon
- 3SABNP, Univ. Evry, INSERM U1204, Université Paris-Saclay, Evry, France
- 4Institut Universitaire de France, Paris, France
- 5Department of Biomedical Engineering, The University of Texas at Austin, Austin, TX, United States
The HIV-1 integrase (IN) is a major target for the design of novel anti-HIV inhibitors. Among these, three inhibitors which embody a halobenzene ring derivative (HR) in their structures are presently used in clinics. High-resolution X-ray crystallography of the complexes of the IN-viral DNA transient complex bound to each of the three inhibitors showed in all cases the HR ring to interact within a confined zone of the viral DNA, limited to the highly conserved 5′CpA 3′/5′TpG 3′ step. The extension of its extracyclic CX bond is electron-depleted, owing to the existence of the “sigma-hole.” It interacts favorably with the electron-rich rings of base G4. We have sought to increase the affinity of HR derivatives for the G4/C16 base pair. We thus designed thirteen novel derivatives and computed their Quantum Chemistry (QC) intermolecular interaction energies (ΔE) with this base-pair. Most compounds had ΔE values significantly more favorable than those of the HR of the most potent halobenzene drug presently used in clinics, Dolutegravir. This should enable the improvement in a modular piece-wise fashion, the affinities of halogenated inhibitors for viral DNA (vDNA). In view of large scale polarizable molecular dynamics simulations on the entirety of the IN-vDNA-inhibitor complexes, validations of the SIBFA polarizable method are also reported, in which the evolution of each ΔE(SIBFA) contribution is compared to its QC counterpart along this series of derivatives.
The HIV-1 integrase (IN) catalyzes the transfer of a viral DNA (vDNA) strand into the genome of the host cell (Lesbats et al., 2016). It is also involved in reverse transcription (Hironori et al., 2009), nuclear import (Mouscadet et al., 2007), and HIV-1 particle maturation (Kessl et al., 2016). It has no counterpart in human cells and thus constitutes an emerging target for the design of novel anti-retroviral inhibitors (Liao et al., 2010).
Three integrase inhibitors have been approved by the FDA in anti-HIV therapies, Raltegravir (RAL) (Summa et al., 2008), Elvitegravir (EVG) (Shimura et al., 2008), and Dolutegravir (DTG) (Underwood et al., 2012). All three act as integrase strand transfer inhibitors (INSTIs) (Ammar et al., 2016). They embody two distinct structural motives, namely a large diketo acid pharmacophore, and a halobenzene derivative. A major advance toward the design of novel derivatives is enabled by high-resolution X-ray crystal structures of the ternary complexes of IN, vDNA, and each of the three INSTI's (Hare et al., 2010a, 2011). These show all three drug complexes to be two-pronged. The keto oxygen and a coplanar neighboring oxygen both coordinate two IN catalytic Mg (II) cations, structural water molecules, and, either directly or through water, IN residues. The halobenzyl moiety is confined in a narrow cleft, binding to the G4 and C16 bases of the highly conserved 5′CpA 3′/5′TpG 3′ step on the viral DNA ends (Hare et al., 2010b). The emergence of IN mutations weakening IN-drug interactions is a major limitation for INSTI-based therapies (Wainberg et al., 2011). Such mutations were reported to occur solely on IN and not on vDNA: this thus leaves open the possibility that additional enhancements of INSTI-vDNA binding should not be adversely impacted by IN mutations.
The present study was motivated by two findings. First, reports from one of our Laboratories showed the DTG > EVG > RAL ranking of affinities for the IN-vDNA complex (denoted as the intasome, INT) to be paralleled by their corresponding affinities for the sole vDNA (El Khoury et al., 2017). Could, thus, increases of the INSTI-INT binding affinities be attempted upon focusing on the sole “ternary” complexes of G4, C16, and a halobenzene ring? Second, recent spectrometric and computational studies of the binding of INSTIs to viral DNA extremities showed the ranking of affinities to be governed by the enthalpy (ΔH) component of the binding free energies (ΔG), the entropy component (TΔS) being similar for all three complexes (El Khoury et al., 2019), a “signature” for intercalative binding (Chaires, 2006). Furthermore, the ΔH ranking of the three inhibitor affinities was itself paralleled by the corresponding ranking of the ab initio quantum chemistry (QC) intermolecular interaction energies, ΔE(QC), of their halobenzene rings with the sole G4/C16 base pair. Upon focusing on the halobenzene ring of the best bound compound, namely DTG, could, then, ΔE (QC) Energy Decomposition Analysis (EDA) along with electronic structure considerations offer insight for affinity-enhancing chemical substitutions?
The full structures of RAL, EVG and DTG are presented in Figure 1. There is a zone of electron depletion along the extension of the C-F bond of DTG, denoted as the “sigma-hole” (Murray et al., 2008). This bond is para to the C-C bond connecting the ring to the central diketo acid group (Figure 1). In the crystal structure of the DTG-INT complex, it points toward the electron-rich bicyclic ring of G4. This could constitute a key stabilizing feature of the DTG-G4 complex.
Figure 1. Molecular structure of Raltegravir (RAL), Elvitegravir (EVG), and Dolutegravir (DTG) represented with the Chimera software. The coordinates used for these representations are taken from the Protein Data Bank crystal structures of the Prototype Foamy Virus (PFV) intasome in complex with magnesium and Raltegravir (PDB code: 3OYA—Hare et al., 2010a), Elvitegravir (PDB code: 3L2U—Hare et al., 2010b) as well as with Dolutegravir (PBD code: 3S3M—Hare et al., 2011). We emphasized with a circle the difluoro halobenzyl ring of DTG, where the substitutions are going to be made to conceive new inhibitors. Carbon atoms are colored tan, hydrogen atoms are colored white, fluorine atoms are colored light green, chlorine atoms are colored dark green, nitrogen atoms are colored blue and oxygen atoms are colored red.
We have considered and analyzed several analogs of the DTG halobenzene ring (HR), 13 of which will be reported in this paper. It is noted that with the exception of compound A1, all compounds have an extracyclic -NHCH3 or -NH2 proton donor replacing in para the extracyclic CF bond. All additional substitutions were done ortho or para to the C-C connector by electron-withdrawing groups. Each substituent can impact ΔE(QC) by a combination of different factors:
- on account of its electron-withdrawing character, a favorable increase of the electrostatic contribution, EC, of ΔE (QC), of the para substituent with G4;
- a favorable increase of the polarization contribution, Epol, of ΔE (QC), due to the contribution of its polarizability to the polarization energy of the ligand;
- when in meta, it could contribute to additional electrostatic interactions with sites belonging to G4 or C16, as reflected by EC and Epol. There is a cone of electron-rich density around the CX bonds of halogens. Halogen substituents could enhance EC if such a cone was in the vicinity of electron-deficient sites of C16.
- these three factors could be counteracted in some cases by increases of the short-range repulsion, EX, and, in the energy balances, by a larger solvation penalty than DTG, since most designed derivatives bear a more polar character than DTG.
EDA unravels the relative magnitudes of the individual ΔE (QC) contributions along the series investigated, and how each contribution can be impacted by each substituent. Optimizing the localized HR-G4-C16 interactions could contribute to a modular design of INSTIs as a preliminary to long-duration molecular dynamics (MD) of the entirety of the INSTI with the entirety of the INT. The binding site includes highly polar protein residues, two Mg (II) cations, and structural, highly polarizable water molecules. Polarizable, multipolar Molecular Mechanics/Dynamics approaches, such as SIBFA (Gresh et al., 2007) or AMOEBA (Ponder et al., 2010) should be adapted to these simulations. Such procedures were demonstrated to reliably account for the impact of the sigma hole and the dual character of the CX bond of HR on electrostatics (El Hage et al., 2013). Very large macromolecular complexes, such as the INT-INSTI ones, are now amenable to long-duration MD upon resorting to the massively parallel code Tinker-HP (Lagardère et al., 2018). Prior to these, and in the context of the present study, it was thus essential to evaluate the accuracy of one of these procedures, SIBFA. Specifically, how well will the evolutions of ΔE(QC) and of each of its contributions, be paralleled by the SIBFA contributions along the series of the thirteen HR ligands?
As a complement to EDAs, we will also report the contours of their electrostatic potential maps.
Quantum Chemistry Calculations
Energy Decomposition Analysis
The decomposition of the ab initio SCF interaction energy is done using the Reduced Variational Space (RVS) analysis (Stevens and Fink, 1987), where the intermolecular interaction energy is separated into four contributions: Coulomb (EC) and short-range exchange-repulsion (EX) in first order (E1) and polarization (Epol) and charge-transfer (Ect) in second order (E2). Finally, the dispersion contribution is assessed as the difference between the BSSE-corrected B97-D3 intermolecular interaction energies and the Hartree-Fock (HF) ones. The basis set superposition error (BSSE) is taken into consideration in the final energy values. The GAMESS software with the cc-pVTZ (-f) basis set (Schmidt et al., 1993) were used in this analysis. We denote below as Epol (VR), a “variational” value of Epol, obtained as the difference between ΔE(RVS) and the sum of E1 and Ect. This is to be contrasted to Epol(RVS) at the outcome of the RVS procedure, the sum of all individual ligand polarization energies computed in a process when the occupied molecular orbitals (MO) of this ligand is relaxed toward its own virtual MO's, the other ligands being frozen.
Energy decomposition analyses were also performed at the correlated levels and the ω-B97-D3 functional (Grimme, 2006) resorting to the Absolutely Localized Molecular Orbitals method (ALMOEDA) (Azar and Head-Gordon, 2011) using the Q-Chem software (Shao et al., 2015). ΔE is decomposed into a “frozen density” component (FRZ), namely the sum of the Coulomb and short-range contribution, and a polarization (Epol) and a charge transfer contribution (Ect) in second order.
The intermolecular interaction energies (ΔE) of the complexes formed by the substituted rings with G4 and C16 were computed at the correlated level using the dispersion-corrected functionals B97-D3 and B3LYP-D3 (Goerigk and Grimme, 2011), with the cc-pVTZ basis set (-f) (Feller, 1996) and the Gaussian software G09 (Frisch et al., 2009) software. The obtained values were corrected for BSSE (Simon and Duran, 1996). The values of the corresponding correlation and dispersion contributions were computed as the differences between the dispersion-corrected correlated ΔE and uncorrelated Hartree-Fock (HF) ΔE values.
G09 was used as well for energy-minimization of the complexes, using a starting structure determined by X-ray crystallography and taken from Protein Data Base site (PDB code: 3S3M, Hare et al., 2011). The C-C bond connecting the HR to the diketoacid (DKA) ring was replaced by a CH bond, since the DKA ring was removed.
The halobenzene rings were relaxed except for the hydrogen atom of this CH bond. This choice was made to avoid wanderings over the G4/C16 bases prevented in the complete complexes by HR anchoring to the rest of the drug. We also chose to hold both guanine and cytosine frozen in their X-ray geometry to account for their anchoring in the DNA backbone. We did not allow for conformational relaxation around the glycosidic bonds of G4 and C16 since we considered that this viral DNA base pair is held in an experimental position “tailored” for DTG: any additional positional relaxation in complexes with derivatives with bulkier ligands could a priori be expected to further optimize, rather than penalize, the binding of such derivatives, with possibly an even more favorable outcome than the one from the present study. This could clearly, only be clarified at the outcome of long-duration MD simulations on the ternary DNA-IN-ligand complexes, enabling to explore the widths and flexibility of the energy basins in the G4/C16/HR zone. Such studies will be reported in due course.
A Continuum desolvation energy ΔGsolv, computed following the Polarizable Continuum Model (PCM) procedure (Mennucci et al., 2002) was also calculated for each HR ring. It was considered as an upper bound to its actual PCM desolvation energy prior to its complexation.
In the context of the SIBFA procedure, the intermolecular interaction energy (ΔEtot) is computed as the sum of five contributions: electrostatic multipolar (EMTP), short-range repulsion (Erep), polarization (Epol), charge transfer (ECT), and dispersion (Edisp). EMTP is computed with distributed multipoles (up to quadrupoles) derived from the QC molecular orbitals precomputed for each individual molecule, derived from the Stone analysis (Stone and Alderton, 1985) and distributed on the atoms using a procedure developed by Vigné-Maeder and Claverie (1988). EMTP is augmented with a penetration term (Piquemal et al., 2007). The anisotropic polarizabilities are distributed on the centroids of the localized orbitals (heteroatom lone pairs and bond barycenters) using a procedure due to Garmer and Stevens (1989). Erep and ECT, the two short-range contributions, are computed using representations of the molecular orbitals localized on the chemical bonds and on localized lone-pairs. Edisp is computed as an expansion into 1/R6, 1/R8, and 1/R10 and embodies an explicit exchange—dispersion term (Creuzet et al., 1991).
Contours of Electrostatic Potentials
The contours of molecular electrostatic potentials (MEPs) of the HRs derived from their electronic densities were displayed by the Gaussview software, implemented in the Gaussian software (Frisch et al., 2009). The colored zones are based on their electronic densities.
Results and Discussion
All DTG derivatives have at least one extracyclic halogen atom; F, Cl, or Br. The DTG coordinates are derived from the crystal structure of the viral intasome with Dolutegravir and two magnesium ions (PDB code: 3S3M, Hare et al., 2011). The first derivative considered, A1, has a chlorine substituent replacing the second DTG fluorine, which is para to the CH bond connecting it to the diketo moiety. This did not bring a significant ΔEtot increase with respect to DTG (see Table 1 below). Several alternative substituents could be considered. However, for the present study we deemed it more instructive to replace upfront the para-F group, which binds G4 by its sigma-hole prolonging the CF bond, by another electron-deficient group, namely a proton donor. We could thus leverage these C-X bonds (X halogen or electron-deficient group) interactions and increase their magnitudes by depleting this donor electron-density with selected electron-attracting substituents. A natural choice bore on substituents such as -NHCH3 or –NH2. Both groups have the additional advantage of acting as electron-donating substituents: this is in contrast to fluorine, which is electron-withdrawing. This could favor the interaction of the electron-deficient sites of C16 with the halobenzene ring and/or with the electron-rich cone around the halogen substituent. It is another means of leveraging the “Janus-like” character of the halobenzene derivatives (El Hage et al., 2014, 2015). The –NHCH3 substituted derivatives were split into two groups, having two or one electron-withdrawing substituents, respectively. The search for –NH2 substituted derivatives was limited to three derivatives, with one or two fluorine substituents. The four groups are indexed as Group 1, 2, 3, and 4, and the derivatives in each group are indexed with capital letters from A to E. All derivatives are represented in Figure 2.
Table 1. Intermolecular interaction energies (kcal/mol) of halogenated derivatives with the G4/C16 vDNA base pair.
Figure 2. Molecular structures of the halobenzene ring of DTG (Group 1) and the newly designed ones (Groups 1, 2, 3, and 4) represented with the Chimera software. Carbon atoms are colored tan, hydrogen atoms are colored white, fluorine atoms are colored light green, chlorine atoms are colored dark green, bromine atoms are colored dark red and nitrogen atoms are colored blue. The squares indicate the new atoms or moieties chosen after substitution.
Para-Electron Attracting Substituted Derivatives
This group contains the compound A1 which has a similar structure to DTG, but with a chlorine atom in para instead of a fluorine one.
Para–NHCH3 Substituted Derivatives
Five derivatives are considered. They all have, similar to DTG, fluorine in ortho, and a third, electron-withdrawing substituent:
A2, B2, and C2 have in the second ortho position: Br, Cl, and –CF3, respectively. D2 has chlorine in meta, and E2 has a cyano substituent in the second ortho position.
Para–NHCH3 Substituted Derivatives
A3, B3, and C3 have as ortho substituent chlorine, bromine and fluorine, respectively, while D3 has fluorine in meta instead of ortho.
Para–NH2 Substituted Derivatives
A4 and B4 have a fluorine in ortho and in meta, respectively, while C4 has the two fluorines in ortho and meta.
Figure 3 gives a representation of the complexes of the G4/C16 vDNA base pair with representative compounds of series 1 to 4. We have displayed the parent compound, DTG from the first group, C2 and E2 from the second, B3 from the third and C4 from the fourth group. As seen below, the inhibitory rings face the ring of C16, interacting with the latter via π-π stacking, while the moiety in para position face the electron-rich rings of G4. Therefore, the electron donor groups in C2, E2, B3, and C4 as well as the electron-depleted region in the prolongation of the sigma-hole in DTG interacts with G4 via electrostatic interactions.
Figure 3. Representation of the G4/C16 complexes with representative HRs. The structures of G4, C16 and DTG are taken from the Protein Data Bank crystal structure of the Prototype Foamy Virus (PFV) intasome in complex with magnesium and Dolutegravir (PDB code: 3S3M—Hare et al., 2011). The newly designed compounds are represented with the Chimera software. The color code for the atoms remains the same as in Figure 1, with the oxygen atoms colored in light red.
Table 1 lists the ΔE (QC) values at the B97-D3, B3LYP-D3, and ω-B97D levels together with the ΔE (SIBFA) values. It also lists the continuum solvation energies, ΔGsolv, of each isolated derivative prior to complexation, and the resulting values of partial B97-D3 energy balances after ΔGsolv subtraction from ΔE (B97-D3). These ΔGsolv values should be considered as upper bounds to the actual PCM solvation energies of the halobenzene ring of DTG and its derivatives, because in the “complete” inhibitor, the exposure of the polar group is less on account of a lesser accessibility and possible intramolecular interactions. The MD calculation of the actual solvation free energies of DTG and its derivatives in a box of water molecules will be considered in the next step of our studies in parallel with the MD simulations of their complexes with the intasome (El Darazi et al., work in progress).
The listed QC ΔE values take into account the BSSE correction as the HF and at the correlated levels. It has values in the 1.5–2 kcal/mol range at the HF level, and 3.5–4.6 kcal/mol range at the correlated levels, depending upon the functional, but changing little for a given functional with the ligand.
The evolution of ΔE as a function of the compound number is plotted in Figure 3 for each functional.
All three ΔE (QC) curves run parallel. Table 1 and Figure 4 show that the ΔE (QC) values of all compounds in series 2–4 have significantly larger magnitudes than those of either DTG or its chlorinated derivative A1. Taking the B97-D3 results as an example, the differences are in the 4–8.3 kcal/mol range. All three QC procedures concur into having five derivatives which stand out from the rest. Four belong to series 2, namely A2, B2, C2, and E2, and one belongs to series 3, namely B3. In series 2, all four derivatives have their two substituents on both ortho sides of the connecting -CH bond, or, equivalently, on both sides meta to the –NHCH3 para substituent. The two best-bound compounds are C2 and E2, with either a –CF3 or a cyano substituent –CN in meta to the –NHCH3 group, respectively. C2 has a larger ΔGsolv value than E2, namely −10.1 compared to −7.8 kcal/mol. As mentioned above, this 2.3 kcal/mol difference should represent an upper bound to the actual desolvation energy difference between the two derivatives. It could reflect a greater “hydrophilicity” of the C2 than the E2 ring. Including such a difference in partial relative energy balances would give rise to a 2 kcal/mol preference for the cyano derivative compared to the trifluoromethyl one.
Figure 4. Evolution of the intermolecular interaction energies (ΔE) of the ternary complex; G4/C16/substituted ring as a function of the halobenzene ring of DTG and the newly conceived rings (Compounds A1 to C4). ΔE values are calculated via quantum chemistry correlated calculations with the B3LYP-D3 functional (ΔE B3LYP-D3-raw before BSSE correction and ΔE B3LYP-D3-corrected after BSSE correction) and the B97-D3 functional (ΔE B97-D3-raw before BSSE correction and ΔE B97-D3-corrected after BSSE correction). ΔE obtained via quantum chemistry energy decomposition analyses methods (EDA) are also presented, precisely the Absolutely Localized Molecular Orbitals method using two functionals; B97D3 (ΔE ALMOEDA-B97D3) and ω-B97D (ΔE ALMOEDA/ω-B97D). ΔE values are obtained as well via molecular mechanics calculations with the polarizable force field SIBFA (ΔEtot SIBFA).
In series 3, the best-bound derivative is B3, with a bromine substituent. Owing to the large hydrophilic character of this atom compared to fluorine in DTG, the values of ΔGsolv are correspondingly larger in series 2 and 3 than that of DTG, in the range 2.8–4.0 kcal/mol. This leaves out relative energy preferences for A2–C2, E2, and B3 over DTG in the range 4–6.7 kcal/mol.
The ΔE values in series 4 with an –NH2 para substituent have smaller magnitudes than in series 2 and 3. These are 4–5 kcal/mol more favorable than DTG, but these are virtually fully compensated for by correspondingly larger ΔGsolv values, leaving partial energy balances <0.5 kcal/mol more favorable than DTG, which is inconclusive in terms of augmented affinities.
As for the BSSE values, we wish to note that the average values of the BSSE amounts to 4.14 kcal/mol for the B97D-3 functional and to 3.41 kcal/mol for the W97D functional, both used for the ALMOEDA energetic calculations. With each, the relative variations among the compounds are around 0.5 kcal/mol. Whereas, at the HF level, the average value of BSSE is in the range of 1.5–2 kcal/mol.
In Figures 5–8 we report the evolution of the individual QC contributions along the series. We denote by E1 the “frozen” energy contribution from ALMOEDA and the sum of the Coulomb, EC, and exchange, EX, contributions from the RVS analyses. E1 (SIBFA) is correspondingly the sum of its penetration-augmented electrostatic and short-range repulsion contributions. We denote by E2 the sum of the polarization and charge-transfer contributions in all approaches.
Figure 5. Evolution of the first order energy (E1) of the ternary complex; G4/C16/substituted ring as a function of the halobenzene ring of DTG and the newly conceived rings (Compounds A1 to C4). E1 values are calculated via quantum chemistry Energy Decomposition Analysis methods (EDA) precisely the Reduced Variational Space method (E1 RVS) and Absolutely Localized Molecular Orbitals method using two functionals; B97D3 (E1 ALMOEDA/B97D3) and ω97D (E1 ALMOEDA/ω-B97D). E1 values are obtained as well via molecular mechanics calculations with the polarizable force field SIBFA (E1 SIBFA).
In the perspective of long-duration polarizable MD on the complexes of vDNA and of the Integrase-vDNA assembly, it is essential to evaluate how well could a procedure such as SIBFA account for the ΔE (QC) trends. This alone could lend credence to comparative energy balances and prospective free energy calculations bearing on complexes out of reach of QC calculations, but amenable to this, and related procedures such as those between “improved” DTG derivatives with v-DNA, let alone with INT, which total several thousands of atoms. Table 1 and Figure 4 show that ΔEtot (SIBFA) is fully able to recover the trends from ΔE(QC). C2 and E2 are found to be the two best-bound compounds, while B2, B3, and A2 come next with a small margin. Such agreements are encouraging, considering that no extra calibration effort was done on the DTG derivatives.
Figure 5 reports the evolutions of E1 at the uncorrelated RVS and correlated B97D3 and ω-B97D levels and those of E1 (SIBFA). Figure 6 reports the corresponding evolutions of E2. Figure 7 compares the evolution of uncorrelated ΔE (RVS) and ΔE (SIBFA) without the dispersion contribution Edisp. Figure 8 reports the evolution of a contribution denoted as “Edisp/Ecorr,” the gain in ΔE upon passing from the uncorrelated RVS ΔE to the correlated B3LYP-D3, B97-D3, and ω-B97D levels, along with Edisp (SIBFA). There is no explicit ALMOEDA “dispersion” contribution, as it is included in the van der Waals kernel for both E1 and E2.
Figure 6. Evolution of the second order energy (E2) of the ternary complex; G4/C16/substituted ring as a function of the halobenzene ring of DTG and the newly conceived rings (Compounds A1 to C4). E2 values are calculated via quantum chemistry Energy Decomposition Analysis methods (EDA), precisely the Reduced Variational Space method (E2 RVS) and Absolutely Localized Molecular Orbitals method using two functionals; B97D3 (E2 ALMOEDA/B97D3) and ω-B97D (E2 ALMOEDA/ω-B97D). E2 values were obtained as well via molecular mechanics calculations with the polarizable force field SIBFA (E2 SIBFA).
Figure 7. Evolution of the interaction energy of the ternary complex; G4/C16/substituted ring as a function of the halobenzene ring of DTG and the newly conceived rings (Compounds A1 to C4). ΔE values are calculated via quantum chemistry Energy Decomposition Analysis method (EDA), precisely the Reduced Variational Space method (ΔE RVS) and via molecular mechanics calculations with the polarizable force field SIBFA (ΔE SIBFA).
Figure 8. Evolution of the dispersion energy (Edisp) of the ternary complex; G4/C16/substituted ring as a function of the halobenzene ring of DTG and the newly conceived rings (Compounds A1 to C4). Edisp values correspond to the ones calculated via quantum chemistry correlated calculations using two functionals; B97D3 (Edisp B97D3) and B3LYPD3 (Edisp B3LYPD3) and via molecular mechanics calculations with the polarizable force field SIBFA (Edisp SIBFA).
Throughout we will denote by ΔE (SIBFA) and ΔEtot (SIBFA) the SIBFA intermolecular interaction energies without and with the dispersion contribution. ΔE (RVS) denotes the BSSE-corrected RVS intermolecular interaction energy. ΔE (QC) denotes the correlated QC intermolecular interaction energy with the B3LYP-D3, B97-D3, or ω-B97D functionals. The QC-derived “dispersion” energy is the difference between ΔE (QC) and ΔE (HF), the latter being derived from a Hartree-Fock computation without removing the BSSE.
Figure 5 shows that the trends in ΔE (QC) and ΔEtot (SIBFA) giving distinct preferences for derivatives A2, B2, C2, E2, and B3 are retrieved by E1, while the E2 curves (Figure 6) are much shallower. Figure 7 shows that similar to ΔE (QC) and ΔEtot (SIBFA), both ΔE (RVS) and ΔE (SIBFA) curves have minima with derivatives A2–C2, E2, and B3. ΔE (RVS) has an additional minimum with derivative A4 from series 4. The corresponding ΔE (SIBFA) minimum is higher in energy. There is a lesser correspondence between Edisp (SIBFA) and those found from the B97D3 and B3LYP-D3 calculations (Figure 8). For compounds A1-E2, it has values intermediate between those from both functionals. For A3 and B3, it is close to the B97D3 values, with differences of 1.2 and 0.2 kcal/mol, respectively. For C3 and D3, it is close to the B3LYP-D3 values, with differences of 0.0 and 0.4 kcal/mol, respectively. For all three compounds in series 4 having the para –NH2 substituent, Edisp (SIBFA) has smaller values than either functional. In fact, the trends are not consistent between the two functionals. Thus, while Edisp (B3LYP-D3) is larger by 2.7 kcal/mol than Edisp (B97-D3), both functionals give very close Edisp values in both B3 and C3 complexes.
Even though it appears satisfactory, the agreement between ΔEtot (SIBFA) and ΔE (QC) could be further improved in the near future. There are ongoing SIBFA refinements in the context of reconstruction of a new library of protein, DNA, and ligand constitutive fragments. They resort to multipoles and polarizabilities from correlated calculations, and a rescaling of the individual energy contributions on the basis of correlated Symmetry Adapted Perturbation Theory SAPT-DFT (Misquitta and Szalewicz, 2002; Podeszwa et al., 2006) energy decomposition analyses. They will be reported elsewhere.
Molecular Electrostatic Potentials (MEP)
It is instructive to represent the impact of some of the reported substitutions on the MEP contours. Figure 9 presents such contours around five representative compounds: DTG; C2 and E2; the meta-substituted trifluoro- and cyano derivatives of group 2, respectively; B3, the meta-substituted bromine derivative of series 3; and C4, the difluorine-substituted derivative of group 4.
Figure 9. Contours of the Molecular Electrostatic Potential (MEP) around: DTG, C2, E2, B3, and C4. The distribution of the electronic charges is 10−3 electron Bohr-3 isodensity surface. The colored bar at the top of the figure indicates the variation of the electron density between −2.094e−2 (corresponding to the color red) and 2.094e−2 (corresponding to the color blue).
There is a much wider extension of the “blue” zone of positive potential around the para –NHCH3 or –NH2 substituents than around the para-fluorine substituent of DTG. This explains the better efficiency observed in the case of the newly conceived compounds in comparison with DTG. In fact, the lower electron density of NHCH3 and NH2 confirms the better electrostatic interaction of these rings with the electron rich polycyclic structure of G4. For the compounds of group 2, the –CF3 substituted derivative C2 generates a much wider and delocalized “red” zone of negative potential than the cyano-substituted derivative E2. The extension of positive potential around the -NHCH3 group is also larger for C2 than for E2. To a large extent, the first feature should explain the 2.1 kcal/mol larger ΔGsolv energy of C2 than E2, impacting the energy balance in favor of E2 while both derivatives had close ΔE (QC) and ΔEtot (SIBFA) values. Derivative B3 also displays a wide zone of negative electrostatic potential around the –Br substituent, and this zone is significantly more extended than around the two –F substituents of derivative C4. The latter presents a low electron density in the middle of its ring due to the presence of its two-electron attracting fluorines, giving it the least favorable attraction with the DNA viral base pairs.
Conclusions and Perspectives
We have designed a series of derivatives of the halobenzene ring of dolutegravir (DTG), the most potent inhibitor of the HIV-1 integrase to date. These derivatives target the highly conserved G4/C16 base pair of viral DNA. The -para fluorine ring of DTG is replaced by a –NHCH3 or a –NH2 substituent. All derivatives reported in this study had more favorable intermolecular interaction energies with this base pair than DTG, as computed by quantum chemistry, and a polarizable molecular mechanics procedure, SIBFA. The two compounds with the highest G4/C16 affinities had a para –NHCH3 substituent and two substituents meta to it: fluorine and trifluromethyl, denoted as C2, and fluorine and cyano, denoted as E2. ΔEtot (SIBFA) displayed trends consistent with ΔE (QC). This study, along with our previous ones (El Hage et al., 2014, 2015; El Khoury et al., 2019) shows that it is possible to design, in a piece-wise fashion, novel derivatives targeting a well-defined, highly conserved, subset of the recognition site of a large macromolecular target. It also constitutes an essential validation step prior to large-scale polarizable molecular dynamics simulations of the complex of the entirety of the drug with the entirety of the target.
Data Availability Statement
The datasets generated for this study are available on request to the corresponding author.
PD, LK, and KH performed computations. PD, RM, ZH, J-PP, and NG analyzed results. PD, J-PP, and NG wrote the paper.
This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Program (grant agreement no. 810367), project EMC2.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We wish to thank the Grand Equipement National de Calcul Intensif (GENCI): Institut du Développement et des Ressources en Informatique Scientifique (IDRIS), Centre Informatique de l′Enseignement Supérieur (CINES), France (Project No. x2009-075009), and the Centre Régional Informatique et d'Applications Numériques de Normandie, Rouen, France (Project 1998053). We thank also for funding this work the Research Council of Saint-Joseph University of Beirut, Lebanon (Project FS71), the Lebanese National Council for Scientific Research, CNRS-L (Projects CNRS-FS80 and CNRS-FS 116), and the French Institute-Lebanon and the French-Lebanese Program CEDRE (Project 35327UJ).
Ammar, F. F., Hobaika, Z., Abdel-Azeim, S., Zargarian, L., Maroun, R. G., and Fermandjian, S. (2016). A targeted DNA substrate mechanism for the inhibition of HIV-1 integrase by inhibitors with antiretroviral activity. FEBS Press 6, 234–250. doi: 10.1002/2211-5463.12025
Azar, R. J., and Head-Gordon, M. (2011). An energy decomposition analysis for intermolecular interactions from an absolutely localized molecular orbital reference at the coupled-cluster singles and doubles level. J. Chem. Phys. 132:024103. doi: 10.1063/1.3674992
Creuzet, S., Langlet, J., and Gresh, N. (1991). Adjustment of the SIBFA method for potential maps to study hydrogen bonding vibrational frequencies. J. Chem. Phys. 88, 2399–2409. doi: 10.1051/jcp/1991882399
El Hage, K., Piquemal, J. P., Hobaika, Z., Maroun, R. G., and Gresh, N. (2013). Could an anisotropic molecular mechanics/dynamics potential account for sigma hole effects in the complexes of halogenated compounds? 34, 1125–1135. doi: 10.1002/jcc.23242
El Hage, K., Piquemal, J. P., Hobaika, Z., Maroun, R. G., and Gresh, N. (2014). Substituent-modulated affinities of halobenzene derivatives to the HIV-1 integrase recognition site. Analyses of the interaction energies by parallel quantum chemical and polarizable molecular mechanics. J. Phys. Chem. A 118, 9772–9782. doi: 10.1021/jp5079899
El Hage, K., Piquemal, J. P., Hobaika, Z., Maroun, R. G., and Gresh, N. (2015). Could the “Janus-like” properties of the halobenzene CX bond (X- Cl, Br) be leveraged to enhance molecular recognition? 36, 210–221. doi: 10.1002/jcc.23786
El Khoury, L., El Hage, K., Piquemal, J. P., Fermandjian, S., Maroun, R. G., Gresh, N., et al. (2019). Spectrometric and computational studies of the binding of HIV-1 integrase inhibitors to viral DNA extremities. Peer J Phys. Chem. 1:e6. doi: 10.7717/peerj-pchem.6
El Khoury, L., Piquemal, J. P., Fermandjian, S., Maroun, R. G., Gresh, N., and Hobaika, Z. (2017). The inhibition process of HIV-1 integrase by diketoacids molecules: understanding the factors governing the better efficiency of dolutegravir. Biochem. Biophys. Res. Commun. 488, 433–438. doi: 10.1016/j.bbrc.2017.05.001
Goerigk, L., and Grimme, S. (2011). A thorough benchmark of density functional methods for general main group thermochemistry, kinetics, and noncovalent interactions†. Phys. Chem. Chem. Phys. 13, 6670–6688. doi: 10.1039/c0cp02984j
Gresh, N., Cisneros, G. A., Darden, T. A., and Piquemal, J. P. (2007). Anisotropic, polarizable molecular mechanics studies of inter- and intramolecular interactions and ligand-macromolecule complexes. a bottom-up strategy. J. Chem. Theory Comput. 3, 1960–1986. doi: 10.1021/ct700134r
Hare, S., Smith, S. J., Métifiot, M., Jaxa-Chamiec, A., Pommier, Y., Hughes, S. H., et al. (2011). Structural and functional analyses of the second-generation integrase strand transfer inhibitor dolutegravir (S/GSK1349572). Mol. Pharmacol. 80, 565–572. doi: 10.1124/mol.111.073189
Hare, S., Vos, A. M., Clayton, R. F., Thuring, J. W., Cummings, M. D., and Cherepanov, P. (2010a). Molecular mechanisms of retroviral integrase inhibition and the evolution of viral resistance. PNAS 107, 20057–20062. doi: 10.1073/pnas.1010246107
Hironori, N., Hayashi, T., Takahashi, T., Miyano, M., Kannagi, M., and Masuda, T. (2009). Augmentation of reverse transcription by integrase through an interaction with host factor, SIP1/Gemin2 is critical for HIV-1 infection. PloS ONE. 4:e7825. doi: 10.1371/journal.pone.0007825
Kessl, J. J., Kutluay, S. B., Townsend, D., Rebensburg, S., Slaughter, A., Larue, R. C., et al. (2016). HIV-1 integrase binds the viral RNA genome and is essential during virion morphogenesis. Cell 166, 1257–1268. doi: 10.1016/j.cell.2016.07.044
Lagardère, L., Jolly, L. H., Lipparini, F., Aviat, F., Stamm, B., Jing, Z. F., et al. (2018). Tinker-HP: a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields. Chem. Sci. 9, 956–972. doi: 10.1039/C7SC04531J
Mennucci, B., Tomasi, J., Cammi, R., Cheeseman, J. R., Frisch, M. J., Devlin, F. J., et al. (2002). Polarizable Continuum Model (PCM) calculations of solvent effects on optical rotations of chiral molecules. J. Phys. Chem. A 106, 6102–6113. doi: 10.1021/jp020124t
Misquitta, A., and Szalewicz, K. (2002). Intermolecular forces from asymptotically corrected density functional description of monomers. Chem. Phys. Lett. 357, 301–306. doi: 10.1016/S0009-2614(02)00533-X
Piquemal, J. P., Chevreau, H., and Gresh, N. (2007). Toward a separate reproduction of the contributions to the hartree-fock and DFT intermolecular interaction energies by polarizable molecular mechanics with the SIBFA potential. J. Chem. Theory Comput. 3, 824–837. doi: 10.1021/ct7000182
Podeszwa, R., Bukowski, R., and Szalewicz, K. (2006). Density-fitting method in symmetry-adapted perturbation theory based on kohn–sham description of monomers. J. Chem. Theory Comput. 2, 400–412. doi: 10.1021/ct050304h
Schmidt, M. W., Baldridge, K. K., Boatz, J. A., Elbert, S. T., Gordon, M. S., Jensen, J. H., et al. (1993). General atomic and molecular electronic structure system. J. Comput. Chem. 14, 1347–1363. doi: 10.1002/jcc.540141112
Shao, Y., Gan, Z., Epifanovsky, E., Gilbert, A. T. B., Wormit, M., Kussmann, J., et al. (2015). Advances in molecular quantum chemistry contained in the Q-Chem 4 program package. Mol. Phys. 113, 184–215. doi: 10.1080/00268976.2014.952696
Shimura, K., Kodama, E., Sakagami, Y., Matsuzaki, Y., Watanabe, W., Yamataka, K., et al. (2008). Broad antiretroviral activity and resistance profile of the novel human immunodeficiency virus integrase inhibitor elvitegravir (JTK-303/GS-9137). J. Virol. 82, 764–774. doi: 10.1128/JVI.01534-07
Stevens, W. J., and Fink, W. H. (1987). Frozen fragment reduced variational space analysis of hydrogen bonding interactions. Application to the water dimer. Chem Phys Lett. 139, 15–22. doi: 10.1016/0009-2614(87)80143-4
Summa, V., Petrocchi, A., Bonelli, F., Crescenzi, B., Donghi, M., Ferrara, M., et al. (2008). Discovery of raltegravir, a potent, selective orally bioavailable HIV-integrase inhibitor for the treatment of HIV-AIDS infection. J. Med. Chem. 51, 5843–5855. doi: 10.1021/jm800245z
Underwood, M. R., Johns, B. A., Sato, A., Martin, J. N., Deeks, S. G., and Fujiwara, T. (2012). The activity of the integrase inhibitor dolutegravir against HIV-1 variants isolated from raltegravir-treated adults. J. Acquired Immunodefic. Synd. 61, 297–301. doi: 10.1097/QAI.0b013e31826bfd02
Vigné-Maeder, F., and Claverie, P. (1988). The exact multicenter multipolar part of a molecular charge distribution and its simplified representations. J. Chem. Phys. 88, 4934–4948. doi: 10.1063/1.454705
Keywords: HIV-1 integrase, quantum chemistry, polarizable force fields, viral DNA, Dolutegravir, INSTIs
Citation: El Darazi P, El Khoury L, El Hage K, Maroun RG, Hobaika Z, Piquemal J-P and Gresh N (2020) Quantum-Chemistry Based Design of Halobenzene Derivatives With Augmented Affinities for the HIV-1 Viral G4/C16 Base-Pair. Front. Chem. 8:440. doi: 10.3389/fchem.2020.00440
Received: 25 November 2019; Accepted: 27 April 2020;
Published: 19 June 2020.
Edited by:Nino Russo, University of Calabria, Italy
Reviewed by:Tiziana Marino, University of Calabria, Italy
Pedro Alexandrino Fernandes, University of Porto, Portugal
Copyright © 2020 El Darazi, El Khoury, El Hage, Maroun, Hobaika, Piquemal and Gresh. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
†Present address: Léa El Khoury, Department of Pharmaceutical Sciences, University of California, Irvine, Irvine, CA, United States