A computer-aided method for controlling chemical resistance of drugs using RRKM theory in the liquid phase

The chemical resistance of drugs against any change in their composition and studying the rate of multiwell-multichannel reactions in the liquid phase, respectively, are the important challenges of pharmacology and chemistry. In this article, we investigate two challenges together through studying drug stability against its unimolecular reactions in the liquid phase. Accordingly, multiwell-multichannel reactions based on 1,4-H shifts are designed for simplified drugs such as 3-hydroxyl-1H-pyrrol-2(5H)-one, 3-hydroxyfuran-2(5H)-one, and 3-hydroxythiophen-2(5H)-one. After that, the reverse and forward rate constants are calculated by using the Rice Ramsperger Kassel Marcus theory (RRKM) and Eckart tunneling correction over the 298–360 K temperature range. Eventually, using the obtained rate constants, we can judge drug resistance versus structural changes. To attain the goals, the potential energy surfaces of all reactions are computed by the complete basis set-quadratic Becke3 composite method, CBS-QB3, and the high-performance meta hybrid density functional method, M06-2X, along with the universal Solvation Model based on solute electron Density, SMD, due to providing more precise and efficient results for the barrier heights and thermodynamic studies. To find the main reaction pathway of the intramolecular 1,4-H shifts in the target molecules, all possible reaction pathways are considered mechanistically in the liquid phase. Also, the direct dynamics calculations that carry out by RRKM theory on the modeled pathways are used to distinguish the main reaction pathway. As the main finding of this research, the results of quantum chemical calculations accompanied by the RRKM/Eckart rate constants are used to predict the stability of drugs. This study proposes a new way to examine drug stability by the computer-aided reaction design of target drugs. Our results show that 3-hydroxyfuran-2(5H)-one based drugs are the most stable and 3-hydroxythiophen-2(5H)-one based drugs are more stable than 3-hydroxy-1H-pyrrol-2 (5H)-one based drugs in water solution.

www.nature.com/scientificreports/ calculation was carried out at the M06-2X/6-311+g(2df,2p) and M06-2X/Jun-cc-pVTZ (and M06-2X/Jun-cc-pV(T+d)Z for S atom) levels of theory for all species of the liquid phase reactions. We used the calculated frequencies to verify the nature of stationary points as true minima, such as reactant (R), intermediate (IN), and product (P) and true first-order saddle points (TSs), and also to gain zero-point energies (ZPEs), partition functions, and thermodynamic parameters (ΔE°, ΔH°, and ΔG°).
To ensure the connectivity between saddle points and corresponding minima, the intrinsic reaction coordinate calculations (IRC) 32,33 were performed at the M06-2X/6-311+g(2df,2p) level in the mentioned phases. Intending to achieve the accurate and more reliable energetic parameters, i.e. relative energies, we performed higher-level calculations by the CBS-QB3 method for all stationary points. All of the electronic structure calculations in this study were done with the Gaussian 09 program package 34 . The images of all molecular structures were drawn by the Chemcraft program (Version 1.7) 35 .

Rate constant calculations.
Here, we discuss the main approach and the used theory for the calculation of the rate of multiwell-multichannel reactions. For multi-step reactions, we can write where l stands for liquid, k 1 , k 2 , k 3 , and k 4 are forward rate constants, and k -1 , k -2 , k -3 , and k -4 are reverse rate constants. Using the steady-state approximation, the overall rate constant for generation of B(l), C(l), D(l), and E(l) are as follows: 36 In Eqs. (2)-(5), the forward and reverse rate constants are calculated by employing the Master Equation/ Rice-Ramsperger-Kassel-Marcus 37 (ME/RRKM) theory. In this research the same as gas-phase reactions 38,39 for calculating the rate constants of multiwell-multichannel reactions, we used the Ssumes program 40 to solve ME/ RRKM for all forward and backward unimolecular reactions. In the framework of RRKM theory, the microcanonical rate constant k(E) for unimolecular reactions with internal excitation energy E is computed by means of very simple expression [41][42][43] where α is the reaction path degeneracy. G*(E − E0) is the integrated density for activated complex without the degree of freedom related to the passage via transition state. E 0 is the critical energy of a reaction. h is Planck's constant. N(E) is the density of the rovibrational states of the reactant.
The RRKM calculation input requires the Lenard Jones parameters (collision diameter, σ (Å), and energy parameter, ε/k B (K)) of reactant and third body. The third body is solvent molecules. As far as we know, the Lenard Jones parameters among colliding molecules studying here are not available experimentally but water. The collision diameter, σ (Å), and the energy parameter, ε/k B (K), for water molecule are 2.74 Å and 506 K 44 . The Lennard Jones parameters for 3-hydroxyfuran-2(5H)-one, 3-hydroxyl-1H-pyrrol-2(5H)-one, and 3-hydroxythiophen-2(5H)-one is calculated approximately. First, an approximation method based on the Sanghvi et al. 45 and Varshni 46 studies is used for derivation of the critical parameters of the investigating molecules. Finally, we used Tee 47 relations for extracting the Lennard-Jones parameters of the target molecules using the following equations where T c (in K) and P c (in atm) are critical temperature and pressure, respectively. The predicted Lennard-Jones parameters for 3-hydroxyfuran-2(5H)-one, 3-hydroxyl-1H-pyrrol-2(5H)-one, and 3-hydroxythiophen-2(5H)one is 5.34 Å and 519.88 K, 5.44 Å and 580.79 K, and 5.34 Å and 553.67 K, respectively.

Results and discussion
The intramolecular 1,4-H-shifts in five-member heterocyclic compounds (that are used frequently in drug designing) are modeled with the CBS-QB3, M06-2X/Jun-cc-pVTZ, M06-2X/Jun-cc-pV(T+d)Z, and M06-2X/6-311+g(2df,2p) levels. All total energies, ZPE corrections, thermal energies, enthalpy energies, and Gibbs free  Tables S1-S6 (see Supplementary Information). Also, the computed reverse and overall rate constants (in s -1 ), calculated frequencies (in cm -1 ), and cartesian coordinates are collected in Tables S7-S16 (see Supplementary Information). We simplify the structure drawn in scheme 1by substituting hydrogen atoms in lieu of R i groups. Thus, we use the following compounds.
The mechanism of 1,4-hydrogen shifts leading to yield the aromatic products is the conventional hydrogen atom transfer. However, all steps of the considered multiwell-multichannel reactions are the conventional hydrogen atom transfer. The most common reaction pathways for the abovementioned compounds are summarized as follows: where X refers to O, S, and N atoms). In Scheme 2, the molecular structures and names of reactants, intermediates, and products are brought. The prefix Y-N, Y-O, and Y-S denote reactant, intermediates, and product in the reaction pathways of 3-hydroxyl-1H-pyrrol-2(5H)-one, 3-hydroxyfuran-2(5H)-one, and 3-hydroxythiophen-2(5H)-one compounds, respectively.
As far as we are aware, neither experimental nor theoretical rate data are available for the abovementioned 1,4-H shifts. This is the first report for evaluating the mechanism and rate constants of intramolecular 1,4-H shifts of 3-hydroxyl-1H-pyrrol-2(5H)-one, 3-hydroxyfuran-2(5H)-one, and 3-hydroxythiophen-2(5H)-one compounds in the liquid phase. Throughout the text, three consecutive energetics and rate constants belong to the CBS-QB3, M06-2X/Jun-cc-pVTZ (or M06-2X/Jun-cc-pV(T+d)Z for S atom), and M06-2X/6-311+g(2df,2p) levels, respectively. Intramolecular H-shifts in 3-hydroxyl-1H-pyrrol-2(5H)-one. Figure 1 shows all bond lengths of the stationary point's structures in 1,4-H shift of 3-hydroxyl-1H-pyrrol-2(5H)-one. The schematic potential energy www.nature.com/scientificreports/ surface for the multiwell-multichannel unimolecular reaction of 3-hydroxyl-1H-pyrrol-2(5H)-one is sketched in Fig. 2. Graph of overall rate constants for the production of all minimum stationary points is displayed in Fig. 3. Thermodynamic parameters such as thermal energies, enthalpies, and Gibbs free energies, and also relative energies corrected by ZPE for all involved species are listed in Table 1. Tables 2 and 3 contain the computed unimolecular rate constants and half-lifes of the intermediates, respectively. All of the mentioned energetic parameters and rate constants are calculated at the CBS-QB3, M06-2X/Jun-cc-pVTZ, and M06-2X/6-311+g(2df,2p) levels. According to Fig. 2 Fig. 3 illustrates that the production rate of IN2-N is lower than IN1-N as expected for multi-step reactions if the next step barrier is higher than the previous. This statement remains true over the 298.15-360 K temperature range. As depicted in Fig. 1 Table 3 Table 1. Relative energies corrected by ZPE, thermal energies (ΔE°), enthalpy energies (ΔH°), and Gibbs free energies calculated for all stationary points of 3-hydroxyl-1H-pyrrol-2(5H)-one unimolecular reaction at different levels. A, B, and C are CBS-QB3, M06-2X/Jun-cc-pVTZ, and M06-2X/6-311+g(2df,2p) levels, respectively.  It is worth noting that the relative energies of all species, but IN4-N, are not sensitive to the applied methods in water. For IN4-N, the computed energetic parameters have a large difference ~ 8.3-9.7 kcal/mol between the CBS-QB3 and DFT levels. This difference causes to obtain different rates and so half-lifes in the computed methods.
After all, the total rates of all minimum stationary points show that all species produced in the course of reaction have the same behavior from the kinetic point of view (see Fig. 3). The generation rates increase with temperature. Moreover, by emphasizing this point that any change in the drug structure leads to arise a change in its properties. And also, Fig. 3 reveals that 1,4 H-shift of 3-hydroxy-1H-pyrrol-2(5H)-one compound occurs hardly. So, it can be concluded that the drugs based on the considered compound are stable against 1,4-H shifts over the considered temperature range.  Figure 5 shows the schematic potential energy surface. The calculated overall rate constants are shown in Fig. 6. Tables 4 contain relative energies (corrected by ZPEs), thermodynamic parameters such as thermal energies, enthalpies, Gibbs free energies. Tables 5 and 6 Table 2. Unimolecular rate constants (s −1 ) calculated by RRKM theory at the CBS-QB3, M06-2X/Jun-cc-pVTZ, and M06-2X/6-311 + g(2df,2p) levels for all channels of 1,4-H shift of 3-hydroxyl-1H-pyrrol-2(5H)-one.    Tables 3, 6 In the last stage of the reaction IN3-O isomerizes to the final product P-O through hydrogen migration from C1 atom to O1 by TS4-O (see Fig. 5). Activation enthalpy for this process is 54.13, 52.70, and 53.40 kcal mol −1 . The same as 1,4-H shift of 3-hydroxyl-1H-pyrrol-2(5H)-one this step is the rate-determining step of the overall Table 3. Half-lifes (in s) of the predicted intermediates in 1,4-H shift of 3-hydroxyl-1H-pyrrol-2(5H)-one at the CBS-QB3, M06-2X/Jun-cc-pVTZ, and M06-2X/6-311 + g(2df,2p) levels.  www.nature.com/scientificreports/ 2(5H)-one compound is displayed in Fig. 8. Figure 9 represents the calculated overall rate constants. Table 7 lists the relative energies and thermodynamic parameters such as thermal energies enthalpies, Gibbs free energies for all maximum and minimum points. Unimolecular rate constants for all intermediate and product formation are tabulated in Table 8, and half-lifes of the reactant the intermediates IN1-S, IN2S, and IN3-S are listed in Table 9.
About the less important pathway (the second pathway), our results show that the initial step of this path is more suitable in comparison with the 3-hydroxyl-1H-pyrrol-2(5H)-one and 3-hydroxyfuran-2(5H)-one unimolecular reactions. The activation free energy for TS5-S is 4.6, 2.35, and 2.26 kcal mol −1 and 16.47, 14.16, 14.01 kcal mol −1 lower than that of TS5-N and TS5-O, respectively. Thus, the rate of IN4-S production is 2.91E+01, 5.92E+00, and 3.06E+00 times faster than IN4-N and is 7.41E+09, 2.05E+08, and 5.17E+07 faster than IN4-O at 298.15 K. Also, the second step of this path has the activation free energy of 12.08, 14.20, and 14.06 kcal mol −1 that is 16.64, 6.79, and 7.42 kcal mol −1 and 8.95, 7.73 and 8.23 kcal mol −1 lower than the activation free energy of TS4-N and TS5-N, respectively. Therefore, it will have a faster rate constant for conversion to the third intermediate. Table 5. Unimolecular rate constants (s −1 ) calculated by RRKM theory at the CBS-QB3, M06-2X/Jun-cc-pVTZ, and M06-2X/6-311+g(2df,2p) levels for all channels of 1,4-H shift of 3-hydroxyfuran-2(5H)-one. and the product P-X are the most stable compounds after reactant R-X in the 1,4-H shifts of 3-hydroxyl-1H-pyrrol-2(5H)-one, 3-hydroxyfuran-2(5H)-one, and 3-hydroxythiophen-2(5H)-one compounds. It is important to note that the intermediates IN2-X cannot be assumed as a product since it has a small half-life and also the forward and reverse free activation energies of these intermediates are lower than the reverse free activation energies of the products P-X (a reaction with a high saddle point). In fact, the product P-X is the major adduct of all investigated reactions because of their attributed heats of formation and Gibbs free energies of reaction, and having reverse high free energy barriers.
As the last point for expressing more clearly the purpose of this study, we mention some concepts here. Our mean by chemical resistance is that when chemical compounds such as drugs do not react from a specified center in defined conditions, we can say that the center has chemical resistance under specified conditions. Also, it is important to note that reactivity and stability are two contrary challenges. Stable compounds are less reactive. But, the condition in which a compound can be assumed stable must be considered. So, for examining the reactivity of a center, we must simulate the most probable reactions of the center by applying reliable conditions. In the sequel, after constituting the potential energy surface, computing the rate of studying reactions should preferably be done to have complete pieces of information. After that, it can be easily judged whether the reaction occurrs or not (in the designed condition).
In summary of the three last sections, it is preferable to say that the calculated PESs and overall rate constants have useful information for the final judge of the drugs' stability. As seen in Figs. 3, 6 and 9, conversion of the reactant to the final product encounters with small rate constant that is in line with the high the energy barriers. This creates a chemical resistance against to drugs' intramolecular changes. Also, the computed forward rate constants uncloak this fact that up to the IN3-X formation may occur at high temperatures, but the computed half-lifes and reverse rate constants indicate that the reaction back to the reactant with high percentage and a very small amount of IN3-X may transform to P-X, which is negligible. Because all of the computed rate constants Table 6. Calculated half-life (in s) for all intermediates in 1,4-H shift of 3-hydroxyfuran-2(5H)-one at the CBS-QB3, M06-2X/Jun-cc-pVTZ, and M06-2X/6-311+g(2df,2p) levels.    Tables 2,5,8), in a condition of the body (37 °C), any of these reactions do not occur. Therefore, all investigated target drugs are stable versus the above-discussed reactions in the human body. But, it should be noted that instead of different functional groups which are available in real drugs we substituted hydrogen atoms for simplification. Also, it has been widely proven that the attachment of substituents to chemical compounds facilitates their reactions. And as aforementioned, this study is an example of drug stability evaluation against its possible unimolecular reactions. In general, we can say that if 3-hydroxyfuran-2(5H)-one, 3-hydroxythiophen-2(5H)-one, 3-hydroxy-1H-pyrrol-2 (5H)-one based drugs have the same functional groups, the first is the most stable and the second is more stable than the third.

Conclusion
A combination of quantum chemical calculations accompanied by the RRKM/Eckart rate constants was used to predict the stability of some drugs. In this respect, the stability of drugs based on a change in their structures was considered mechanistically and kinetically. In the mechanistic part, multiwell-multichannel reactions based on 1,4-H shifts in water were designed for pyrrolidinone-based drugs as an example. So, three target drugs simplified as 3-hydroxy-1H-pyrrol-2 (5H)-one, 3-hydroxyfuran-2(5H)-one, and 3-hydroxythiophen-2(5H)-one were  composite method, and the meta hybrid M06-2X method in conjunction with two well-behaved basis sets, Juncc-pVTZ (or Jun-cc-pV(T+d)Z) and 6-311+g(2df,2p), in water. Also, for more clarification of the mechanism, we reported the thermodynamic parameters of the mentioned compounds and the associated stationary points in the ground-state potential energy surface using high-level theoretical calculations. In the kinetic part, for the formation and consumption of the intermediates INT1-X-IN4-X, and the product P-X, the steady-state approximation assumption accompanied by the computed thermal rate constants were used to calculate the www.nature.com/scientificreports/ overall rate constants for all intermediates and the main product over 298.15-360 K temperature range. Therefore, the RRKM theory along with Eckart tunneling correction was utilized to calculate forward and reverse rate constants. Also, the overall rate constants were calculated by using Eqs. (2)-(5). Our calculated potential energy surface accompanied by the predicted rate constants indicated that the main reaction pathway for the formation of P-X was the first pathway and the rate of the second pathway due to having a large energy barrier has a very small contribution to the generation of the product P-X. In addition, the large energy barrier in the reverse direction of the last step, P-X → IN3-X, indicates that the reaction has a small likelihood to return if P-X is formed (see Tables 1, 4, 7). It is merit pointing out that our results demonstrated that the hetero atom plays a key role in a hydrogen atom jumping on five-membered heterocyclic compounds due to make different energetic results and rate constants in the same carbon skeleton. Therefore, it affects directly the stability of drugs. For example, our results showed that 3-hydroxyfuran-2(5H)-one based drugs are the most stable, and 3-hydroxythiophen-2(5H)-one based drugs are more stable than 3-hydroxy-1H-pyrrol-2 (5H)-one based drugs in water. Table 8. Unimolecular rate constants (s -1 ) calculated by RRKM theory at the CBS-QB3, M06-2X/Jun-cc-pV(T + d)Z, and M06-2X/6-311 + g(2df,2p) levels for all channels of 3-hydroxythiophen-2(5H)-one.   Table 9. Temperature dependent half-lifes (in s) of the predicted intermediates in 1,4-H shift of 3-hydroxythiophen-2(5H)-one at the CBS-QB3, M06-2X/Jun-cc-pV(T + d)Z, and M06-2X/6-311 + g(2df,2p) levels.