Abstract

This study is part of the search for solutions to the resistance of the parasite to antimalarials. The aim of this work is to analyze the relative stability and reactivity of adducts in order to understand the resistance of the parasite to artemisinin. The DFT/B3LYP/GenECPs method associated with bases 6-31G, 6-31G ** and LANL2DZ is used to perform these calculations. On one hand analysis of the thermodynamic quantities of formation and reaction revealed that β adduct is the most stable among the studied adducts. In other words, the latter is the majority product of the interaction between heme and artemisinin. On the other hand, the adduct γ is the least stable. The values of the pameters such as dipole moments, intramolecular hydrogen bond lengths and thermodynamic parameters (formation and reaction quantities) have shown that it can exist a correlation between the following three parameters those are the relative stability of adducts (α, β, δ), their dipole moment and the binding intramolecular hydrogen. The analysis of the Frontier Molecular Orbitals also showed that the β adduct is susceptible to electrophilic attack, unlike the γ adduct, which is rather favorable to nucleophilic attack and can be the most excited component. Moreover, the determination of the alkylation rate will be considered in order to specify the order of complexation of the meso positions α, β and δ of the heme.

1. Introduction

Artemisinin is an antimalarial which is the basis of a lot of antimalarial drugs currently used in first-line treatment ^{1, 2, 3, 4, 5}. This antimalarial alkylates the heme on its meso positions α, β, γ and δ to form adducts ^{6, 7}. This alkylation is weak in γ and strong in α, β and δ positions. In 2001, WHO (World Health Organization) considered artemisinin to be "the world's biggest hope in the struggle against malaria" ⁸. Unfortunately, it was a short hope because since 2007, artemisinin-resistant parasites have appeared in Southeast Asia ^{9, 10}. This molecule does not completely eliminate all malaria parasites, hence the need to combine it with other antimalarials ^{11, 12}. But the first cases of clinical failures with artemisinin-based combination therapies (ACTs) have been identified in the same area (Southeast Asia) ^{13, 14, 15}. The problem of parasite resistance to antimalarials remains despite the introduction of new molecules. This study is a side of the research for solutions to the resistance of the parasite to antimalarials. It assumes that the instability of the adducts can explain the resistance of the parasite to artemisinin. The aim of this work is to analyze the relative stability of heme-artemisinin adducts (Figure 1) in order to understand the ineffectiveness of this antimalarial. The theoretical approach of adduct reactivity is based on the theory of Frontier Molecular Orbitals (FMO). The calculations are performed in gaseous phase and in solution using the DFT/B3LYP/GenECPs method associated with bases 6-31+G (d,p), 6-31G and LanL2DZ.

2. Materials and Methods

Theoretical studies are carried out with Gaussian 09 calculation software ¹⁶. The different structures were optimized using density functional theory (DFT) with B3LYP / GenECPs as the level of calculation ^{17, 18}. The atoms of carbon, nitrogen and oxygen (C, N and O) are optimized at B3LYP/6-31 + G (d, p) theory level. The atoms of hydrogen (H) and iron (Fe) are optimized respectively at the B3LYP/6-31G and at B3LYP/LanL2DZ level. LanL2DZ is a pseudo-potential base ^{19, 20}. All these calculations are made by using the same method (DFT). An artemisinin derivative was added at one of the meso positions of the heme unit to form adducts in order to improve the representation of the coordination environment found in the parasite vacuole. As the alkylation reaction takes place in an acid environment, two solvents were used for the solution calculation, namely water and acetic acid. The Conductor-like Polarizable Continuum Model (CPCM) is the solvation model used for solution calculations ²¹. The complete optimizations were carried out without any symmetry constraint. The harmonic vibration frequencies have been calculated to confirm that the optimized geometry corresponds exactly to a local minimum which does not possess a negative frequency.

The knowledge of the variations of energy contributions to the internal energy at 0 K and at 298.15 K under normal atmospheric pressure (P = 1atm.) between products and reagents contributes to the energetic characterization of a chemical reaction ^{22, 23}. For a given energy parameter X, its variation is determined according to the following relation (1):

(1)

The studied energetic parameters are the variations of the electronic energy of reaction (internal energy of this reaction), the variation of the enthalpy and the variation of the free enthalpy . To get access to the internal energy, the contributions of different movements such as translation, rotation and vibration of the nuclei must therefore be taken into account in order to evaluate them. In the ideal gas approximation, the rotational and translational contributions are given by relation (2):

(2)

The contribution ZPVE (Zero Point Vibrational Energy) which means the energy of the lowest vibrational level due to the 3N-6 normal modes of vibration (3N-5 for the linear molecules) of frequencies of N nuclei at 0 K, is defined by the relation below (3):

(3)

To obtain the energy corresponding to 298 K, it is necessary to take into account the additional energy due to the population of the vibration levels during the temperature rising from 0 to 298 K, defined by the relation (4) below:

(4)

At the end, the variation of internal energy at 298 K can be written like relation (5):

(5)

The formation’s enthalpy of reaction at 298 K corresponds to the internal energy variation which has been corrected from the term Δ (PV), ie (–RT) because the variation in the number of moles is equal to (-1). This variation is given by the following expression (6):

(6)

The estimation of the free reaction enthalpy is then obtained by equations (7) and (8).

(7)

(8)

HOMO’s orbitals usually act as electron donors and LUMO’s ones as electron acceptors. They can give a reasonable qualitative prediction of the excited state properties of a molecule and the capability of a molecule to transport electron. They are also quantum chemical parameters those determine molecular reactivity ²⁴. The energy gap noted ΔE is the energy difference between the energy levels of HOMO (the highest occupied molecular orbital) and LUMO (lowest vacant molecular orbital). This difference (ΔE) serves as a measure of the molecule’s excitability. Thus, the lower the energy gap is, the more the molecule can interact with the environment. A large gap HOMO-LUMO implies a high stability for the molecule and low reactivity in terms of chemical reactions. Similarly, a small difference implies a high reactivity of the molecule ^{25, 26, 27}.

3. Results and Discussion

Table 1 contains the enthalpies of formation (Δ_fH^298K) in atomic unit, the reaction enthalpies (Δ_rH^298K) and the free reaction enthalpies (Δ_rG^298K) in Kcal.mol^-1, the lengths of intramolecular hydrogen bonds (d) in angström and the dipolar moments (μ) in Debye of adducts α, β, γ and δ.As with tables and equations, figures should be set in one column if possible unless two-column display is essential. The resolution of graphics and image should be adequate to reveal the important detail in the figure.

The comparison of the enthalpies of formation mentioned in Table 1 indicates that the alkylation of heme on the meso position β gives the lowest energy and the γ position leads to the highest energetic value. α and δ adducts’ energies are very closed to the same energetic value so they have a similar relative stability. The enthalpies of formation of these adducts increase in acetic acid and decrease in water which means that these adducts are more stable in water than in the acidic area. β adduct remains the most stable molecule and the least stable is γ adduct because they possess respectively the lowest and the highest values of energy. These results permit us to establish the following descending order of stability in acidic environment: β > δ > α > γ. In water, α adduct is more stable than δ adduct and it allow us to get the following order of decreasing stability: β > α > δ > γ. The instability of the adduct γ can be explained by a steric gene (trouble) due to the presence of carboxyl groups. Robert et al. ^{6, 7} showed that the alkylation rate of heme is high in α, β and δ positions, and low in the γ’s one. However, their work does not give any comment on the relative stability of these adducts. The results of our work allow us to establish the order of relative stability of these molecules. The alkylation of heme in the β position leads to the most stable molecule. We can argue that the major product of the interaction between heme and artemisinin is β adduct.

The dipole moment (μ) indicates the stability of a molecule in water, in particular in aqueous solution. Thus, a strong dipole moment will result low solubility in organic solvents and high solubility in water ^{28, 29}. Like enthalpies of formation, the dipole moment increases in the acid environment. α adduct has the highest dipole moment. The relative stability of α, β and δ adducts increases with the dipole moment that is to say more the value of the dipole moment of the adduct is high, less the adduct is stable. Thus, α adduct with the highest dipole moment value (6.1887 and 6.9855 Debye respectively in water and in acid), is the least stable among the studied adducts. The most stable is β adduct because it has the smallest value of dipole moment (3.9948 and 4.1574 Debye respectively in water and in acid). There is a correlation between the relative stability of these adducts and their dipole moment.

The free enthalpy values and calculated enthalpies of reaction (Δ_rH^298K and Δ_rG^298K) are all positive (Table 1). The alkylation at 298.15 K is an endothermic and non-spontaneous reaction at this temperature. In other words, the alkylation reaction takes place at high temperatures.

Concerning H bond lengths (distance d), the practice consist to consider a contact as a H bond if the distance d is smaller than the sum of Van der Waals radii, by taking 1.52 Å for the oxygen atom and 1.0 Å for hydrogen one; that means d ≤ 2.52 Å ^{30, 31}. It is also known that the (link) H bond is strong when its length is short ^{32, 33, 34}. Moreover, studies have shown that the hydrogen bond is stronger when the energy of the reaction is weak ^{35, 36, 37}. The analysis of the intramolecular hydrogen bond across the H bond lengths of the studied adducts shows that the length of the H bond is shorter in the β adduct. It varies from 1.688 to 1.704 Å (Table 1). In α and δ compounds, the values of the hydrogen bond are closer and are comprised between 1.726 Å and 1.733 Å. Concerning the standard free reaction enthalpy, low values are also recorded around the β adduct in water and in acetic acid, which are respectively 21.062 and 22.406. K.cal.mol^-1. All these observations show that there is a strong existence of intramolecular hydrogen bond around the β adduct. This H bond is established between a carboxyl group and the hydrogen atom of the hydroxyl group. Contrary to the β adduct, the hydrogen bond in the other compounds occurs between the two hydroxyl groups (Figure 2). The groups that interact in hydrogen bond formation differ from one adduct to another. This observation assumes the existence of heme conformers based on intramolecular hydrogen bonding between hydroxyl and carbonyl groups. The different structures of the heme conformers are shown in Figure 3. The β-adduct which is the most stable molecule has the smallest distances of hydrogen bond lengths. α and δ compounds with neighboring energy values also have very close hydrogen bonding values. The less stable γ compound does not have hydrogen bond. This observation supposes that there is a correlation between the relative stability of α, β and δ adducts and the hydrogen bond.

Figure 4 shows the different structures of the heme. Calculations of heme’s conformers are carried out in gaseous phase. The energies of enthalpies of formation, boundary orbitals, dipole moments and hydrogen bond lengths are summarized in Table 2.

The results show that the C3 conformer has the lowest value of formation enthalpy (Table 2). It is therefore the most stable. The least stable is the C4 conformer. This assumes that in the gas phase, the heme is in the form of the C3 conformer. The order of relative decreasing stability of these molecules is as follows: C3 > C1 > C2 > C4. The conformers C2 and C4 have the highest dipole moment values. The energy gaps indicate that the conformer C3 is the most stable. The dipole moment of C3 conformer is relatively weak (2.160 Debye) compared with those of the C1, C2 and C4 molecules whose dipole moments are 6.876, 10.081 and 10.408 Debye respectively (Table 2). The relative stability of the conformers increases with their dipole moment. The mentioned values in Table 2 also show that the shorter the hydrogen bond length, the lower the energy of the molecule. These results highlight the probable existence of a correlation between the relative stability of these molecules, their dipole moment and their intramolecular hydrogen bonding.

The energies of the HOMO and LUMO boundary molecular orbitals were calculated at the DFT/B3LYP/genECPs level. The gap ΔE = E_LUMO - E_HOMO is estimated from these energies. The results are shown in Table 3.

HOMO levels of adducts increase and LUMO levels decrease with acidic solvent (Table 3). The molecular orbital boundary HOMO reflects the electro-donor (nucleophilic) character of the molecule. The higher the energy of this OM, the easier the molecule will yield electrons. As for LUMO, it reflects the electro-acceptor (electrophile) character of the molecule. The lower the energy of this molecular orbital, the easier the molecule will accept electrons. The β adduct has the highest energy of HOMO level so this molecule will be susceptible to electrophilic attack. The γ adduct’s LUMO level energy is the lowest one. This molecule can be favorable to a nucleophilic attack. The smallest value of ΔE is obtained with the adduct γ. The β adduct leads to the largest HOMO-LUMO gap. The γ adduct is susceptible to be more reactive than the other adducts and β adduct is the most stable among the studied molecules. The gap values (ΔE) are close to one another (Table 3) and close to that of the adduct γ. This reflects the neigh bourhood of reactivity of these molecules. The continuation of the alkylation reaction which leads to an alkylated heme molecule with a less bulky group (Figure 4) can be said. Thus this last structure could more easily undergo a dimerization. This could explain the resistance of the parasite to artemisinin. These may be able to reduce the steric hindrance (trouble) effect of the artemisinin moiety and thus facilitate the formation of heme dimers.

4. Conclusion

Analysis of formation enthalpies showed that the β-position alkylation of heme gives the most stable adduct. This adduct seems to be the major product of the heme-artemisinin interaction. That in position γ leads to the least stable adduct and the most reactive. This could justify the low level of alkylation of heme observed experimentally in this position. Our results showed a probable correlation between the relative stability of α, β, δ adducts, their dipole moment and intramolecular hydrogen bonding. The energy gap values indicate that the reactivity of these molecules are all closed. In this sense, the possibility of a reaction that would transform these adducts into other molecules is possible. These may be able to reduce the steric hindrance effect of the artemisinin moiety and thus facilitate the formation of heme dimers. Our results also showed that the groups involved in hydrogen bond formation differ from one adduct to another. The results of this work have opened other avenues of scientific investigation such as the probable existence of a transformation of adducts into other molecules. The determination of the alkylation rate in the framework of the study of the mechanism of the reaction of the heme-artemisinin interaction could bring more precision to this first conclusion. All these lines of investigation would contribute to a probable explanation of the resistance of Plasmodium to artemisinin.

Acknowledgments

The calculations were carried out by Dr. Affoué Lucie Bédé during a post-doctoral research internship in the Department of Physical and Theoretical Chemistry of Alexandru Ioan Cuza University of Iasi-Romania, led by Professor Ionel HUMELNICU, as beneficiary of a postdoctoral fellowship under the "Eugen Ionescu" program 2015-2016 of AUF BECO. We thank them for allowing the calculations in this publication.

References