Study of the Reactivity of (100) Felodipine Surface Model Based on DFT Concepts ()
1. Introduction
The chemical reactivity of a compound can be interpreted as the resistance or ease with which it attracts or gives away electrons under the action of an external potential v(r). In this sense, there are global parameters that allow us to characterize this electron transfer from a theoretical point of view, such as the electronic chemical potential (μ) [1] , the molecular hardness (η) [2] and the electrophilicity index (ω) [3] [4]. Since local descriptors of reactivity, such as the Fukui function (FF) [5] [6] [7] and local softness (s), are properties that depend on the position (r), they may explain selectivity in certain parts of a molecule.
Felodipine (FD), (methyl-4-(2,3-dichlorophenyl)-1,4-dihydro-2,6-dimethyl- 3,5-pyridinedicarboxylate), shown in Figure 1, is a calcium channel blocker and an effective and widely prescribed medication in the treatment of hypertension and angina [8] [9].
Four crystalline forms of FD have been reported [8] , all shaping an orthorhombic system. Form I, used in some commercialized products, is the most stable, where its crystalline structure, determined in 1986, has a space group P2(1)/c, and a unit cell of a = 12.086 Å, b = 12.077 Å, c = 13.425 Å and β = 116.13˚ [10]. Many examples have been described in which the compound polymorphs show variation in reactivity [11] [12] [13] [14]. Furthermore, in the same polymorph, the reactivity can be anisotropic with respect to the crystalline faces where remaining chemicals are different, whereby the reaction rates along particular crystallographic directions can vary significantly. Paul and Curtin were pioneers in conducting solid state reaction studies, specifically the acid-base character of gas-solid reactions, such as with ammonia gas and p-chlorobenzoic anhydride crystals [15] [16]. It was observed that, often, certain crystal faces are attacked by gases preferentially and the reactions propagate along specific directions. However, there are still few reported studies on the solid state reactivity of organic crystals at the electronic level. Currently, the use of computational methods as an indispensable tool for the understanding of physical and chemical properties has led to important advances in the molecular sciences, and in the last decades, in particular, there has been an increased interest in the concepts of Density Functional Theory (DFT) [14] [17] - [22]. Important examples are the recent study by Luty et al. based on Fukui’s nuclear functions using DFT to investigate the explosive mechanism of RDX (hexahydro-1,3,5-trinitro-1,3,5-triazine) [17] , and the study conducted by Shaoxin Feng and Tonglei Li, determining that the Fukui nuclear function can be used to characterize the difference in the chemical reactivity of two polymorphs of flufenamic acid [14]. Since a chemical reaction is driven by the change in energy of the system and it is accompanied by electron transfer and atomic displacement, the Fukui nuclear function, a local function for describing the sensitivity
Figure 1. Molecular structure of felodipine (FD), showing the atom-labelling scheme.
of the system to a simultaneous disturbance in the number of electrons N, may be useful to characterize the reactivity of crystals with respect to crystal packing [14].
The compounds of the group of 1,4-DHP drugs undergo hepatic metabolization (by which the oxidation of 1,4-dihydropyridines in pyridines occurs) catalyzed by the CYP3A enzyme of cytochrome P450 [23] [24]. The metabolization reaction takes place in a two-step mechanism: in the first step the hydrogen is abstracted from the dihydropyridine, where as in the second a hydroxyl group is added to the substrate.
There are many computational approaches that attempt to explain the process of metabolism. Some of them use statistical models in molecular descriptors [25] [26] , other works use ligands to predict the preferred sites of hydrogen abstraction by quantum chemistry [27] [28]. Another, simpler, yet computationally expensive, way is to chemically model the complete process of hydroxylation catalyzed by cytochrome [29] [30] [31]. The hydroxylation reaction has been modeled by coupled quantum/molecular mechanics (QM/MM) [31]. These calculations strongly support the “two-state reactivity” (TSR) model, in which, the active site of the cytochrome―where the oxidation takes place―is a ferro-oxyl species called “Compound I” (CpdI) from the current knowledge [32]. Cpd I can be seen as an “electrophilic oxidant” [29]. Therefore, the Fukui function (
) should help to pinpoint sites in a molecule that prefer to be attacked by Cpd I. On the contrary, the Fukui function
, evaluated for Cpd I prefer to be attacked by dihydropyridine [26].
Felodipine shows a particularly rich variety of metabolites [33] [34]. Most involve the oxidation of 1,4-dihydropyridine to pyridine, sometimes together with the loss of one of the esters, which is well reflected in
[26]. For FD and other compounds, Singh et al. [28] estimated the abstraction energies of the H in a semi-empirical way and its correlation with quite successful hydroxylation sites. Several experiments were carried out to study the structure and properties of felodipine. However, few investigations on the crystal structure faces are available to date.
In this report we investigate how crystal morphology can relate to properties of FD. The novelty of this work lies in studying the surfaces (100), (010) and (001), especially the first one, with supercells to model the periodic systems, based on DFT with the Generalized Gradient Approximation (GGA) and a dispersion correction as formulated by Tkatchenko and Scheffler (TS). By using the FF the most reactive sites of the surfaces to electrophilic species are described.
2. Methodology
We will lay the basic definitions of the global descriptors according to DFT. The energy E is expressed in terms of the number of electrons N and an external potential v(r), so that
, where
is the electronic density [35] [36] [37] [38]. The derivatives of
with respect to N and v(r) lead to a set of global and local properties that allows us to quantify the concept of reactivity and site selectivity. As is well known in the literature, global descriptors are properties that explain the stability of a molecule, such as, the electronic chemical potential (µ), molecular hardness (η) and electrophilicity index (ω) [1] - [7]. The variation of energy with respect to the external potential gives local information, i.e., it depends on the position (r), and therefore it is defined as an index of selectivity. Local descriptors such as the Fukui function (FF) [39] [40] [41] and the local softness (s) [42] are properties that explain the selectivity of a region in a molecule. The FF is defined as:
(1)
Equation (1) presents a discontinuity problem in atoms and molecules when combined with the finite difference approximation, resulting in expressions of the FF [43] [44] :
,
(2)
when a molecule accepts electrons, they tend to move around places where
is large because at these locations the molecule is most likely to stabilize additional electrons, and therefore, it is susceptible to nucleophilic attack at such sites. Likewise, a molecule is susceptible to electrophilic attack at sites where
is large, since these are the regions where electron removal least destabilizes the molecule. In chemical density functional theory, the FF are the key of region selectivity indicators for electron-transfer controlled reactions. The function quantification is possible through a system of condensation on an atomic region, where FF can be written by using population analysis techniques [45]. This leads to the following equations:
,
(3)
where
denotes the electronic population of atom k of the reference system, more correctly called
.
The Fukui functions favourably determinate the reactive sites for most chemical systems. However, the values of the FF rely upon the scheme chosen to calculate the charges and the accuracy of the population analysis used. In this study the Hirshfeld method was used [46] , which has been proved to be accurate for organic systems [47].
The evaluation of the surface energies of the crystalline faces can be useful to compare how the resistance of different surfaces affects the kinetics of the reaction. As often observed, a solid state reaction can proceed in a specific crystallographic direction. Therefore, the study of mechanical properties on different faces can give some insight into the solid state reaction. The surface energy can be calculated as [14] :
(4)
where Eslab and Ebulk are total energies of the surface and crystal in bulk, respectively, and n is the thickness (or layers of unit cells) of the surface, and S is the surface area. In the present study, the surface energies, the electronic structures, and Fukui nuclear functions of three surfaces of the crystal form I of FD were calculated.
The crystal structure of FD, form I, was obtained from the Cambridge Structural Database (CSD) (ref code: DONTIJ). FD crystallizes in an orthorhombic lattice with space group P2(1)/c, and cell parameters a = 12.086 Å, b = 12.077 Å, c = 13.425 Å and β = 116.13˚, packed with four molecules [10].
A periodic solid state program was used, with DFT-D (dispersion correction) [48] [49] [50] [51] [52] at the GGA level developed by Perdew-Burke-Ernzerhof (PBE) [52] , with plane waves as a basis set with a cutting energy of 380 eV and a tolerance for the SCF cycle of 10−6 eV/atom. Vanderbilt pseudopotentials [53] were used to model ion-electron interactions with: H: 1s1, C: 2s22p2, N: 2s22p3, O: 2s22p4, Cl: 3s23p5. Surface (100) (a = 12.077 Å, b = 13.425 Å, c = 32.736 Å and α = β = γ = 90˚), Surface (010) (a = 13.425 Å, b = 12.086 Å, c = 43.082 Å and α = β = 90˚, γ = 116.13˚) and Surface (001) (a = 12.086 Å, b = 12.077 Å, c = 42.005 Å and α = β = γ = 90˚) of form I were modeled (Figure 2) and their Fukui nuclear functions FF were analyzed as they may be linked to their chemical reactivity.
3. Results and Discussion
The electronic structure of the three surfaces of Form I of FD was studied. Also, the FF and the surface energy for each were calculated. The bulk crystal structure was optimized with the same methods that were used to calculate the FF. The network parameters were set during the optimization. Surface models of (100), (010) and (001) faces, of two unitary cells of thickness were thus constructed. On the (100) and (001) faces the rings of the pyridine are exposed on the surface unlike the face (010) where the benzyl ring together with the methyl ether group are more exposed.
We focus our attention on the (100) surface model which turned out to be the most stable, in which there are intermolecular hydrogen bonds formed by N-H-O located on the normal plane of the surface (see Figure 3).
Once the optimization of the geometry of the surface model (100) was carried out, it was observed that the hydrogen bonds lying on the surface tend to move downwards, along the z-axis. The above can be seen in Table 1, where the values corresponding to the midpoint of the hydrogen bonds are given as well as the displacement after optimizing such surface of FD form I. The hydrogen bonds corresponding to those in the layer in between have a displacement along the
Figure 2. Surfaces (100), (010) and (001) of form I of the felodipine. The light blue dotted line shows the hydrogen bonds of the surface models.
Figure 3. Intermolecular hydrogen bonds of the surface (100) of FD.
Table 1. Midpoint coordinates of the hydrogen bonds and their displacement after optimizing surface (100) of felodipine form I.
z-axis lower than those on the surface. From the above it can be observed that surface molecules, having no interaction with other molecules in the upper part, tend to flatten the surface. This displacement of the hydrogen bonds causes an activation of the carbon atoms C2 and C4, which favors the reactivity of these atoms as shown by Fukui functions (Table 2), making them more susceptible to nucleophilic substitutions.
Figure 4 shows the boundary orbitals HOMO (Highest Occupied Molecular Orbital) and LUMO (Lowest Unoccupied Molecular Orbital), from which a high density of unoccupied states in the three surfaces can be seen, located mainly in the carbon atoms 1 and 5 of the 1,4-pyridine ring, which gives us an idea of the places most likely to interact with electron donor atoms.
The FF of each atom was obtained from the calculations of the neutral and anionic forms of the surface models. The molecular structure of the anionic form remained the same as its neutral counterpart. The quantitative results of the FF for faces (100), (010) and (001) are listed in Table 2. The electrophilic attack of the crystalline faces is illustrated in Figure 5. Fukui nuclear functions of
Table 2. Calculated Nuclear Fukui Functions of the (100), (010) and (001) slab models of Form I of felodipine with DFT-D.
Figure 4. HOMO and LUMO orbitals of surfaces (100), (010) and (001) of FD form I.
Figure 5. Illustration of Fukui nuclear function (
blue and
red) for (001), (010) and (100) surfaces of felodipine Form I.
the atoms N1, C1 and C4 of the pyridinic ring and the atoms O1 and O3 of the carbonyls groups in the three surfaces are larger than those of other atoms. These results agree with those reported by Michael E. Beck [34].
In Figure 5 it can be seen that the (100) surface has a distribution of molecules favoring electrophilically reactive sites, thus making this surface susceptible to greater reactivity to agents such as Cpd I, which, as already mentioned, is the active site of the cytochrome, where the oxidation of dihydropyridine to pyridine is carried out.
To elucidate the influence of mechanical resistance on chemical reactivity, the surface energies of the surface models of Form I were calculated by the DFT-D with the functional GGA, obtaining the following Esurf: −0.2303 eV/Å2 for (100), 0.0222 eV/Å2 for (010) and 0.0302 eV/Å2 for (001) surfaces. The latter energies remain significantly above the one for (100) surface, thus indicating a closer bond between the molecules on the faces (010) and (001). As the reaction begins on the surfaces, it is expected that their propagation and penetration in the bulk are limited by forces of intermolecular nature. The surface energy characterizes the intermolecular interactions within a crystallographic plane, specifying the mechanical resistance. Therefore, the lower surface energy of (100) may facilitate a faster reaction rate than the other two surfaces
4. Conclusions
The reaction capacity of the three analyzed surfaces of form I of felodipine was investigated by studying their electronic structures, nuclear Fukui functions and their surface energies. The present findings show that Fukui nuclear functions constitute a useful tool in the analysis of surface reactivity for a crystal such as FD. In addition, due to the highly heterogeneous nature of a molecular crystal reaction, the surface intermolecular forces ought to be taken into account to elucidate chemical reactions occurring on this type of crystals.
These results can provide information on experimental work in surface catalysis as based on theoretical knowledge of local reactivity of the compound here analyzed, thereby saving efforts when selecting the best sites a priori. Our findings may also be useful in some pharmaceutical applications of felodipine.
Acknowledgements
One of us (CTC) gratefully acknowledges a PhD fellowship granted by CONACyT (Mexico). JFRS and AFR are grateful to Sistema Nacional de Investigadores (SNI, Mexico) for provision of funds.