Monitoring Organic Synthesis via Density Functional Theory

Abstract

A preliminary molecular structure for a system, which may or may not be known, is the first step in a typical investigation using ab initio techniques. A stable system is generated by a geometry search using an energy minimization method (usually a local minimum or transition state). Subsequently, it is easy to obtain any energetic properties (such as atomization energies, formation temperatures, binding energies) or expectation values or quantifiable quantities from the wave function of the molecular system and its fragments. The stability of such a system can be determined by considering the second derivative of the energy with respect to the spatial coordinates (also known as the Hessian matrix). It could be a goal to find out how the system interacts with other systems and eventually to decipher the synthesis pathways. Therefore, this chapter presents a recent application of approaches based on density functional theory (DFT) to study chemical processes at the catalytic sites of enzymes. The focus is on the interaction of small organic molecules with the ability to inhibit a catalytic cysteine of the malaria parasite, in the area of drug design.

1. Introduction

Two quantum-chemical theories have been established since Lewis first proposed the idea of a chemical bond at the beginning of the twentieth century [1]: the valence bond (VB) theory and the molecular orbital (MO) theory, both of which are based on Schrödinger’s equation [2]. The density functional theory (DFT), which states that the ground state energy of a non-degenerate N-electron system is a special functional of the density (r), was developed in the 1960s of the previous centuries on the basis of the Hohenberg and Kohn theorems.

$$\mathrm{E}\left[\mathrm{\rho }\left(\mathrm{r}\right)\right]=\int \mathrm{\rho }\left(\mathrm{r}\right)\mathrm{v}\left(\mathrm{r}\right)\mathrm{dr}+\mathrm{F}\left[\mathrm{\rho }\left(\mathrm{r}\right)\right]$$E1

The “external one-electron potential” or the electron-nucleus Coulomb interaction is represented by v(r), and F[(r)] is the Hohenberg-Kohn universal functional obtained by adding the kinetic energy functional T[(r)] and the energy functional of the electron-electron interaction V_ee[(r)]. The rigorous theoretical basis of DFT is this theorem. If the number of electrons is kept constant within the DFT framework, the electron density can be represented as a functional derivative of the energy with respect to the external potential:

$$\mathbf{\rho }\left(\mathit{r}\right)={\left(\frac{\mathbf{\delta }E}{\mathbf{\delta }v\left(\mathit{r}\right)}\right)}_{\mathrm{N}}$$E2

Density functional theory (DFT) calculations therefore require the construction of an expression for the electron density. Similar to the quantum chemical theory based on the Schrödinger equation, it is computationally impossible to resolve the electron density functional ρ(r) for a complicated system. The definition of the individual terms in the functional F[ρ] is the crux of the mathematical problem. The Hartree-Fock equation served as a rough analogy to introduce the Kohn-Sham formalism [3]. In recent years, numerous empirical DFT functions have been developed, including B3LYP [4, 5], MPWB1K [6], and more recently M06 and related functions [7], which provide precise energies and allow the study of organic reactions with a computational cost comparable to MO calculations.

A thorough quantum chemical study of molecular electron density in terms of nonbinding and binding molecular areas was made possible by the invention of topological analysis of the electron localization function in the late twentieth century [8]. This investigation made it possible to build a molecular model that was connected to the Lewis binding pattern. The molecular mechanism of the majority of organic reactions has been identified through topological electron localization function analysis of binding alterations throughout a reaction pathway. The construction of a reactivity model in which these bonds are produced by the C-to-C coupling of two pseudo radical centers [9] formed along the reaction pathway [10] has been made possible by several studies of organic reactions in which C-C bonds are formed. This pattern is interestingly present in C-C double bond reactions that are nonpolar, polar, and ionic. High activation energies in nonpolar reactions are a result of the energy needed to break C-C double bonds in order to produce pseudo radical structures. It is noteworthy that when the reaction’s polarity rises, these high activation energies diminish. Domingo has shown that the global electron density transfer that takes place in polar processes favors the variations in electron density necessary for the creation of C-C single bonds.

An essential and typical first step in most quantum chemical studies is structural optimization. It is a crucial element of any computational chemistry study that deals with the structure and reactivity of a molecule. There are numerous methods for structural optimization. These techniques are used to find transition structures (TS), find minimum energy paths (MEP) that correlate with reaction pathways, and optimize equilibrium geometries. In this chapter, a concise protocol is presented to understand how to easily obtain energetic properties (such as atomization energies, formation temperatures, binding energies) by considering the first and second derivatives of the energy with respect to the Hessian matrix.

1.1 Local optimization

The last 50 years have seen rapid progress in optimization techniques for ab initio molecular orbital simulations. The introduction of energy gradient techniques, the development of algorithms, and the increase in processing capacity have played an important role [11, 12, 13]. Optimization algorithms based solely on energy are orders of magnitude slower than analytical gradient-based optimization strategies. Searching for transition structures is now feasible, and optimization of equilibrium geometries has become routine, even for quite large systems.

1.1.1 Potential energy surface

The relative positions of the individual atoms within a molecule define its molecular structure. The molecule has a clear energy for a particular location and electronic state. The potential energy surface (PES) describes this energy, which varies depending on the electronic state and atomic coordinates. The potential energy surface of two geometric variables is roughly shown in Figure 1. The Born-Oppenheimer approximation [15], in which the motions of nuclei and electrons are studied independently, leads to the ideas of potential energy surfaces. Nuclei move slower than electrons because they are much heavier. This makes it possible to separate the nuclear motions from the electronic motions.

Figure 1.

The minima, saddle points, and inflection points are displayed as an intriguing feature of the PES model [14].

The minima, maxima, and saddle points—the stationary points that form the surface of the potential energy—define it [16]. Each point, which can be identified from the first and second derivatives of the molecule, indicates a different state of the molecule. The first energy derivatives of each atom with respect to its coordinates combine to give a vector called a gradient. The combined second derivatives are used to form a matrix called a Hessian. A vanishing gradient is a property of any stationary point. Also, the Hessian matrix is positively defined at a minimum (all eigenvalues of the Hessian matrix are positive) and the Hessian matrix has only one negative eigenvalue at a first-order saddle point. In chemistry, minima indicate stable structures, while first-order saddle points can be associated with transition states (TS). For structural optimization, higher-order saddle points on the potential energy surface are not relevant (but they are relevant for electronic structure calculations). When using mass-weighted coordinates, the minimum energy path (MEP) or intrinsic reaction coordinate (IRC) consists of the steepest descent reaction paths (SDP) from the transition state to the minima on either side of the saddle point. Taylor expansion allows to represent the potential energy surface, E(x), as an infinite sum in neighborhood of a point x₀, using the step x and gradient g^T vectors, and Hessian matrix, H,

$$\mathit{\text{E}}\mathit{\text{(x)}}=\mathit{\text{E}}{\mathit{\text{(x}}}_{\mathit{\text{0}}})+{\mathit{\text{g}}}^{\mathit{\text{T}}}\mathit{\text{x}}+1/2{\mathit{\text{x}}}^{\mathit{\text{T}}}{\mathit{\text{H}}}_{\mathit{\text{x}}}+\dots$$E3

Most optimization methods are constructed based on this representation. Since two coordinate systems connected by a linear transformation are equivalent in all gradient optimization methods, the significance of the coordinates for the optimization process may not be immediately apparent. This is because the gradient vector and the Hessian matrix must be properly transformed from one system to the other. On this basis, it has been argued that the Cartesian coordinates are comparable to the valence-type internal coordinates [17, 18, 19] or even better [20]. This claim [17, 18, 19] ignores the potential importance of cubic and higher valence couplings. However, it has been shown that optimization in Cartesian coordinates can be as effective as in Z-matrix coordinates [20] or natural internal coordinates for medium-sized systems with a trustworthy initial Hessian matrix and suitable initial geometries [21]. However, optimization in Cartesian coordinates is incredibly wasteful, especially if the system is large enough, since it lacks initial curvature information (i.e., without a Hessian matrix) [20, 21]. Even if the initial geometries are suboptimal for optimization, Cartesian coordinates are less effective than internal coordinates, even if the exact Hessian matrix is always known [22]. A good and affordable initial Hessian matrix for ab initio calculations is usually a molecular mechanics force constant matrix.

1.2 Monitoring transition states

Numerous methods have been proposed to find suitable starting geometries because the calculations of the transition structures are so sensitive to the starting geometry. One very helpful method is to start with reactant and product structures, which are easier to obtain than transition structures. The simplest way to determine the shape of a transition structure is to assume that each atom is exactly halfway between its initial and final positions. The term “linear synchronous transit” refers to this nearly linear motion (LST). Although this is a good first approximation, it is not perfect. Consider the motion of an atom whose bond angle varies relative to the other molecules in the system. The bond length in the middle of the line connecting its starting and ending points will be shorter than expected and therefore have a larger energy (perhaps much larger). The quadratic synchronous transit (QST) method is the logical development of this technique. These techniques assume that a parabola connecting the geometries of reactant and product is formed by the positions of the atoms in the transition structure. Even if it is only a very small improvement, QST is usually better than LST. A weighting factor can often be entered by the user (e.g., to specify a structure containing 70% products and 30% reactants). This allows the application of the Hammond postulate, which states that the transition structure resembles the reactants in an exothermic reaction or the products in an endothermic reaction [23]. These methods have their limitations, but have proven to be quite useful for simple reactions. The drawback is that each of these methods, however good, is based on the premise that the reaction proceeds in a single step with coordinated motion of all atoms [24]. For multistep reactions, each of these methods can be used on its own. For a reaction with a single transition structure but uncoordinated motion (e.g., breaking one bond and building another), the use of hand-drawn initial geometries or eigenvalue tracking may be preferable.

1.2.1 Reaction coordinate

The maximum on the reaction pathway is by definition a transition structure. There is a single, well-defined reaction pathway, namely the intrinsic or low energy (IRC) pathway. Numerous groups have developed methods to derive the IRC pathway from quasi-Newton optimization. These methods are based on the observation that quasi-Newton methods oscillate around the IRC path from one iteration to the next.

By following the reaction path from the equilibrium geometry to the transition structure, one can also obtain a transition structure. The user chooses which mode of oscillation to trigger a reaction given sufficient kinetic energy, which is why this method is also known as eigenvalue tracking [25]. This is not the most effective approach to obtaining an IRC, nor is it the fastest or most reliable way to find a transition structure. The advantage is that there are no assumptions about the transition structure or the coordinated motions of the atoms.

Using a pseudo-reaction coordinate is another method. This method can be a lot of work for the user and takes more time than most other methods. However, it has the advantage of being extremely reliable, so it will work even when all other methods have failed. By first selecting a geometric parameter that is directly related to the reaction, a pseudo-reaction coordinate is calculated (e.g., the bond length for a bond that is formed or broken). Then, a series of calculations is performed while all other geometric parameters are optimized and this parameter is fixed at various values ranging from those in the reactants to those in the products.

The result is an approximation of the reaction coordinate rather than the actual reaction coordinate, which is close to the actual reaction coordinate only for equilibrium geometries and transition structures. Normally, the geometry for a quasi-Newton optimization starts with the calculation from this set with the highest energy. Quasi-Newton optimization may still fail in some rare cases where extremely flat potential surfaces are present. In this case, the transition structure can be estimated with arbitrary accuracy (within the theoretical model) by determining the maximum energy by adjusting the selected geometric parameter in smaller and smaller steps.

1.2.1.1 Solvent effects

The choice of solvent can influence the reaction rate. Interactions with the solvent often change the geometry of the transition structure only slightly, but have a large effect on the energy of the structure. If the solvent effects are taken into account in the calculation, all transition structure discovery methods can be applied. The transition structure identification method is not affected by the presence of solvent interactions.

1.2.2 Evaluating the Hessian

Partitioned-rational function optimization (P-RFO), a useful technique, is based on partitioning the eigenvalues of the Hessian into modes with negative curvatures along which the energy is maximized and all other modes along which the energy is minimized. To guarantee convergence to the nearest transition states, P-RFO requires an initial Hessian with a single negative eigenvalue along the reaction coordinate. If all eigenvalues of the hessian are positive, the P-RFO search may not find the desired transition states because the smallest positive eigenvalue is selected as the mode for energy maximization, which may not match the reaction coordinate. The same problem occurs when the Hessian has numerous negative eigenvalues of comparable magnitude. In this case, the most negative eigenvalue is tracked upward, even though it may not be the reaction coordinate. Even if the transition states estimate is accurate, a P-RFO calculation may fail if the Hessian is not calculated with high accuracy, since inaccuracies in the sign of the eigenvalues create the same ambiguity in determining which mode to follow upward [26]. Monitoring a reaction rate using the activation energy in the Arrhenius equation is the simplest technique for determining a reaction rate. Experimental results or a simple theoretical approach, such as the kinetic theory of gasses, can be used to determine the pre-exponential factor. The activation energy can be roughly calculated by subtracting the energy of the reactant and transition structures. The addition of the zero-point vibrational energy leads to an easily determined additional correction to these energies.

The simple use of activation energy implies that the intrinsic reaction coordinate is the only direction in which a reaction can proceed. It would be more accurate to consider the possibility that reactions can also occur that pass through a geometry that closely approximates the transition structure. The variable transition state calculations take this into account. The vibrational frequencies of the transition structure, the full reaction coordinate, or the total potential energy surface may need to be used in these calculations. Tunneling at the reaction barrier may also be considered in these calculations. These calculations can provide accurate answers, but they are particularly susceptible to minute parameters such as the choice of a mass-weighted coordinate system for the geometry.

Dynamic analyses can be used to study how the direction and velocity of the entering reactants affect the reaction rate. These studies often begin with ab initio calculated potential energy surface. Although obtaining the potential energy surface was not an easy task in itself, the effort required to study a reaction using these methods can be significantly greater.

2. Structural model based on computational method

One of the conventional methods for modeling matter based on quantum mechanics, including atoms, molecules, solids, nuclei, and both quantum and classical fluids, is density functional theory (DFT). The first category includes techniques often referred to as ab initio techniques, such as Hartree-Fock (HF), Moller-Plesset perturbation theory (MP), configuration interaction (CI), and linked-cluster methods (CC) [27]. Given that the wave function for a system of N electrons depends on 3N spatial variables, it should be noted that as the number of electrons increases, so does the complexity of the wave function. DFT approaches [28], on the other hand, use functionals of the electron density which are a function of only 3N spatial variables and independent of the system size.

The first H-K theorem shows how the Hamilton operator, and consequently all the properties of the system, is uniquely determined by the electron density, while the second theorem asserts that the functional related to the ground-state energy of the system yields the lowest energy if and only if the input density is the true ground-state density (i.e., nothing other than the variational principle). Currently, the correctness and effectiveness of DFT-based approaches depend on the basis chosen for the expansion of the Kohn-Sham orbitals, but in particular on the caliber of the exchange correlation (XC) functional used [3].

Local density approximation (LDA) was the first method originally proposed for modeling XC functionals [29]. A simple generalization of the local spin density approximation (LDA) to include electron spin in the functional was proposed by Slater [30]. The well-known SVWN functional, whose exchange component was created by Slater [31] and whose correlation component was created by Vosko, Wilk, and Nusair [31] is indeed one of the early LSDA-XC functionals. However, experience has shown that LSDA leads to inflated binding energies and underestimation of barrier heights, although it usually provides the bond lengths of molecules and solids with an amazing accuracy (about 2%). Thus, the average accuracy of LSDA is insufficient for most applications in monitoring organic synthesis in chemistry.

Compared to LSDA functions, GGA (Generalized Gradient Approximation) functions have been shown to provide better predictions for total energies, atomization energies, and structural properties. The most commonly used XC functions of this type are the Perdew 86 (P86), Becke (B86, B88), Perdew-Wang 91 (PW91), Laming-Termath-Handy (CAM), Perdew-Burke-Ernzerhof (PBE), revised Perdew-Burke-Ernzerhof (RPBE), Perdew-Burke-Ernzerhof revised for solids (PBEsol), Becke exchange and Lee-Yang-Parr correlation (BLYP), and Armiento-Mattsson 2005 (AM05). However, they still provide too low barrier heights and usually fail in describing van der Waals interactions.

An alternative strategy has been developed, called the SCC-DFTB (self-consistent-charge density functional tight binding, later abbreviated as DFTB) method. The DFTB approach, which is an alternative to the traditional semi-empirical methods in quantum chemistry, including the well-known MNDO, AM1, and PM3 schemes derived from RF theory, is based on a second-order expansion of the DFT total energy expression. Since its parameterization technique is based entirely on DFT calculations and does not require fitting empirical data, DFTB is not, strictly speaking, a semi-empirical method.

3. Case study: monitoring enzymatic inhibition

Typically, a two-step mechanism is used to describe how the enzyme and substrate function. First, favorable placement of the molecule in the binding pocket ensures the formation of a non-covalent enzyme-inhibitor complex. When the thiolate group of the catalytic CYS attacks one of the carbon atoms of the ring, ring opening occurs (Figure 2), leading to irreversible alkylation of this amino acid.

Figure 2.

Schematic representation of the ring opening mechanism of the alkylation of a methyl thiolate by an N-substituted aziridine ring. X = O, N-R. Adapted from [27].

The effect of the environment at the carbon atoms in the rings on the kinetics and thermodynamics was evaluated by Helten and colleagues using quantum chemical simulation [32, 33]. According to these estimates, acidic media have a greater effect on the opening of the aziridine ring than for epoxide. They discovered that the protonation of the nitrogen center of aziridine-based inhibitors occurs at the beginning of the reaction course, more precisely before the transition state for the ring opening. It was also discovered, and this is important, that aziridine is inactive without the prior N-protonation step. Therefore, it is suggested that protonation significantly accelerates the reaction rate by stabilizing both the transition state and the ring opening product. Hence, substituents that deprive the nitrogen atom of electrons should lead to lower reaction barriers and, consequently, a greater inhibition effect [32]. Since the products of this reaction are highly stabilized, only the thermodynamics of this reaction are favored by O-protonation, while the kinetics remain unchanged, in contrast to the behavior of aziridines [32]. This is due to the fact that, in the case of epoxide-based inhibitors, the protonation of the oxygen center occurs only after the appearance of the transition state. The attacking cysteine was represented by a methyl thiolate (H₃C-S-), while the inhibitors were represented by tripartite ring systems (H₂C)₂X with X = O, N-R. The authors characterized this simplified model as typical for the inhibition mechanism of the enzyme. The heteroatom of the three-membered rings and the methyl thiolate were placed near solvent molecules with increasing proton donor capacity to capture the effects of decreasing pH on the reaction profile. Water molecules were used to model settings with weak proton donor capabilities (pKa = 15.74). In contrast, NH₄⁺ (pKa 9.3) and HCO₂H (pKa 3.8) molecules were used to mimic environments with stronger proton donor capabilities [33]. While the energies were derived using B3LYP [21] single point calculations, the geometry optimizations of the relevant stationary points, verified by frequency calculations, were calculated using the BLYP [34, 35] functional. A TZVP basis set and both functionals were merged [36]. The accuracy of the theoretical approach was evaluated and the researchers found that BLYP significantly underestimated the barrier heights, while the response profile derived with B3LYP showed excellent agreement with CCSD (T) data [32].

A schematic representation of the mechanism of action at the catalytic site is formed by a cysteine and a histidine, whose side chains form a thiolate/imidazolium ion pair, and an asparagine, which plays a crucial role in the appropriate alignment of the ion pair [37], is outlined in Figure 3. One of the key events in the catalytic hydrolysis of hemoglobin is the nucleophilic attack of the thiolate anion on the corresponding electron-deficient carbonyl group of the substrate, resulting in a negatively charged tetrahedral intermediate stabilized by the “oxyanion hole” formed by the side chains of GLN36 and TRP206 in FP2 and by GLN38 and TRP208 in FP3 [38]. The kinetics, thermodynamics, and regioselectivity of the ring-opening reaction of epoxide- and aziridine-based compounds can be evaluated by using standard quantum mechanics (QM) calculations [32, 33].

Figure 3.

Schematic reaction mechanism of the cysteine protease-mediated cleavage of a peptide bond.

The enzyme, the warhead of the inhibitor, and water are the components of the model system used for the quantum mechanics/molecular mechanics (QM/MM) calculations. The part subjected to quantum mechanics corresponded to the zwitterionic side chains of the catalytic residues CYS and HIS and the electrophilic warhead of the inhibitors. The calculations in the QM domain both with and without the mediating water molecule can identify if water could mediate proton transfer. One-point calculations at the B3LYP/TZVP level of theory are used to more accurately determine the activation and reaction energies, while QM calculations can be performed during geometry optimization and potential energy surface scanning. The MM part can be implemented using the CHARMM force field [39]. The inhibition reaction, which included the ring opening of the inhibitor and the formation of the S(CYS)-C bond, can be captured by the calculated potential energy surface. The results rule out a direct proton shift from the HIS to the organic ligand and show that even a single water molecule is sufficient to create a highly effective relay system that allows protons to move from the HIS to the inhibitor. This result strongly suggests that the effectiveness of the inhibitors can be reduced by substituents that can prevent proton transfer.

The speed-up of the calculations can be even greater when a parameterized quantum mechanics technique such as the density functional tight binding (DFTB) approach is considered or when it is combined with molecular mechanics. This makes QM/MM MD simulations in the ns range easily feasible, which is an essential requirement for reliable determination of free energies for biomolecular reactions.

4. Conclusions

The purpose of this chapter is to present a brief protocol of conceptual DFT reactivity in enzymatic reactions. The calculations can be used to easily determine energetic properties (such as atomization energies, formation temperatures, binding energies) by considering the first and second derivatives of the energy with respect to the Hessian matrix. Chemical intuition comes into play during the procedure, such as understanding the shape of transition states in a system related to the system under study. It is likely that the transition states can be found if the conjecture is strong enough to indicate that it lies in the quadratic basin of transition states, or in other words, has a negative Hessian eigenvalue. The main advantage of this homolog-based method is that it is computationally free and effective when dealing with simple or very closely related systems with known transition states. On the other hand, with the help of realistic models and DFT techniques, enzymatic catalysis is becoming more understandable at the electronic level. This is particularly helpful for the development of irreversible enzyme inhibitors that can covalently bind the catalytic amino acids of the enzyme.

Acknowledgments

N.N acknowledges funding from the Ministry of Higher Education Malaysia for FRGS grant, KPT (FP118-2020)/1/2020/STG04/UM/01/1, and RU grant (ST018-2021).

Conflict of interest

The authors declare no conflict of interest.