A Computational Study of Microhydrated N-Acetyl-Phenylalaninylamide (NAPA): Kinetics and Thermodynamics ()
1. Introduction
Water plays a significant role in changing the structural, kinetic, and thermodynamic properties of peptides, proteins, nucleic acids, lipids, carbohydrates, and natural products [1] [2] [3] [4] . As well as regulates the biological activities of such biomolecules [5] [6] . However, it is still unknown how water molecules control the structures and biological functions. Therefore, it becomes a challenging issue in biophysical science from the experimental and theoretical points of view and also becomes a promising research era in which the aim is to understand the physicochemical properties and biological/pharmaceutical activities at the molecular level. Experimentally it is very difficult to understand the working mechanism due to the extremely rapid phenomena occurring in the bulk water as well as in the biomolecule-water interface. Nevertheless, some experimental techniques have been developed to study the biomolecule-water interface and reported that the molecules are more tightly connected in the interface than in the bulk and hence, fluctuations occur slowly in the interface [7] [8] .
In contrast, computational approaches to biomolecule-water systems have been of great importance to the experimental society because of complementary and independent information to accurately interpret the nature of biomolecule-water interactions. Several experimental [9] [10] [11] [12] and theoretical [13] [14] [15] researches have been performed to interpret the structural, thermodynamic, and kinetic behavior of biomolecule-water interactions using explicit and implicit solvent models. Recently, Lanza et al. reported by comparing and analyzing the computed thermodynamics data of Ala3H+·22H2O and Val3H+·22H2O clusters that the β-strand structure of Val3H+·22H2O clusters is more stable than Ala3H+·22H2O clusters due to the hydration enthalpy and entropy [16] . Arsiccio et al. reported that they used a novel computational approach to study protein unfolding using pressure in an implicit solvent model and found that the volume reduction of the bound water is the key energetic factor for the denaturation of proteins under the effect of pressure [17] . Lanza et al. reported by computationally studying AcAlaNH2·nH2O (n = 1 - 13) complexes with the peptide in the fully extended (FE) conformation that the interaction energy increases upon increasing the number of water molecules [2] . Fischer et al. recently reported on the reaction of protonated triglycine (Gly3H+) with a single water molecule in the gas phase both experimentally and theoretically [18] . The highest barrier height of the transition state (TS) for cis isomer is computed to lie at 1.8 kJ/mol below the reactants. On the contrary, the highest barrier height of the TS for trans-isomer is found to lie 7.0 kJ/mol below the reactants.
The studies of kinetics and thermodynamics of peptide-water interactions are focused on intermediate states (transition states), intermolecular interactions (cooperative hydrogen bonding), and free energy barriers [19] . The addition of water molecules in proteins, protein-ligand interfaces, protein-drug interfaces, or biopolymers modulates thermodynamic properties such as enthalpy, entropy, and Gibb’s free energy depending on the number of hydrogen bond presence [20] [21] . As Gibb’s free energy of activation depends on enthalpy and entropy [22] , therefore, the addition of water in peptides or proteins changes the free energy of activation and finally changes the rate of reaction. In this manuscript, we studied the kinetics and thermodynamic properties of N-acetyl-phenylalaninylamide (NAPA) which is a model dipeptide, and its microhydration using the explicit solvent model in the gas phase.
2. Computational Details
The conformational analysis of NAPA was explored with the aid of AMBER force field [23] implemented in the HyperChem professional 7.51 package [24] . The structure of reactants (R), transition states (TS), and products (P) were computed employing density functional theory (DFT) calculations. The DFT computations were implemented at dispersion correction functional (wB97XD) in conjunction with the cc-pVTZ basis set in the gas phase. Vibrational frequency analyses were performed for the minimum energy structures to evaluate the character of stationary points and also to obtain zero-point vibration energy (ZPVE). The characteristics of local minima for reactants and products were proved by ascertaining that the second derivatives energy of matrices (Hessian) have no imaginary frequency. Transition state geometries were found by QST3 (quadratic synchronous transit-guided quasi-Newton approach) method and confirmed that Hessian matrices have an imaginary frequency. Thermal corrections were done to get thermochemical energies such as enthalpies and free energies at 298.15 K under atmospheric pressure. The IRC calculations were also performed using the same level of computational approach for all the transition states (TS) to prove whether these TS are associated between the reactants and products to the accurate minima or not. Total calculations were done with the aid of Gaussian16 computational program [25] and GaussView 6.0 was applied to visualize the geometry of the compounds. The rate of the reaction was calculated by renowned transition-state theory (TST) [26] via the below-mentioned equation:
where k(T) is the rate of reaction at temperature T, kB is the Boltzmann constant (=1.380662 × 10−23 J/K), h is the Planck’s constant (=6.627176 × 10−34 J∙s), R is the gas constant (=8.31441 J/mol∙K), c0 is the standard concentration (=1 mol/L), and Δ‡G0 is Gibb’s free energy of activation.
3. Results and Discussion
3.1. Geometry Optimization and Thermodynamic Stability
By applying MP2/6-31 + G* and DFT/wB97XD/cc-pVTZ computational levels, N-acetyl-phenylalaninylamide (NAPA) molecules can assume four different conformations and they are NAPA-A, -B, -C, and -D. Among them, the minimum energy conformer is NAPA-A with the structural motif of βL(a) and that is proved both experimentally as well as theoretically [27] [28] [29] [30] [31] . To avoid complexation, only NAPA-A conformer has been taken for microhydration. Several conformers are calculated for mono-, di-, tri-, and tetra-hydrated NAPA-A complexes using DFT/wB97XD/cc-pVTZ computational approach. The most stable conformers for all the hydrated NAPA complexes are taken into account for kinetics and thermochemical analysis in the gas phase. Figure 1 shows the most stable structures of NAPA-A and its hydrated cluster. The Gaussian 16.0 program package is utilized to carry out thermochemical analysis at the default temperature (298.15 K) and pressure (1 atmosphere). All the thermodynamical parameters of hydrated NAPA-A complexes are computed and presented in Table 1.
![]()
Figure 1. DFT-optimized the most stable structures of [NAPA-A(H2O)n (n = 1, 2, 3, 4)] complexes.
![]()
Table 1. Calculated thermodynamic parameters of the most stable structures of NAPA-A, water and [NAPA-A(H2O)n (n = 1, 2, 3, 4)] complexes with wB97XD/cc-pVTZ computational approach.
H = Enthalpy (kcal/mol), G = Gibb’s free energy (kcal/mol), C = Heat Capacity (cal/mol-Kelvin) and S = Entropy (cal/mol-Kelvin).
The total energy for a molecular system is assumed by Etotal = E0 + Evibrational + Erotational + Etranslational, where E0 is the sum of Eelectronic and ZPE. This ZPE to the electronic energy considers for the effect of vibrations occurring even at 0 K in the molecule. Therefore, zero point corrected vibrational energy E, enthalpy H (H = E + RT), and Gibb’s free energy G (G = H − TS) are accounted for the thermal correction. The values of H, G, Cv, and S of NAPA-A are 160.952, 124.195, 54.414, and 125.263 in the respective units mentioned in Table 1. Successive addition of water molecules in NAPA-A increases the value of H, G, Cv, and S at an almost constant rate. The total number of vibrational modes increases with the addition of water molecules. Ultimately the total electronic energies of [NAPA-A(H2O)n (n = 1 - 4)] complexes are increased and therefore enthalpy, Gibb’s free energy, heat capacity, and entropy increase.
3.2. Formation of Hydrated NAPA-A Complex and Thermochemical Analysis
The energy of a reaction is a very important parameter to understand the nature of the reaction. The reaction for the formation of [NAPA-A(H2O)n] complexes are given below:
NAPA-A + (H2O)n (n = 1, 2, 3, 4) → [NAPA-A(H2O)n (n = 1, 2, 3, 4)] (1)
The energy of the reaction can be calculated by subtracting the total electronic energies of the bare NAPA-A and isolated water from the total electronic energy of [NAPA-A(H2O)n (n = 1, 2, 3, 4)] complex. The equations for the calculation of the reaction energy (ΔEr), the change of enthalpy (ΔHr), the change of Gibb’s free energy (ΔGr), the change of heat capacity (ΔCr) and the change of entropy (ΔSr) are shown below:
ΔEr = E[NAPA-A(H2O)n (n = 1, 2, 3, 4)] − E[bare NAPA-A] − E[(H2O)n (n = 1, 2, 3, 4)]
ΔHr = H[NAPA-A(H2O)n (n = 1, 2, 3, 4)] − H[bare NAPA-A] − H[(H2O)n (n = 1, 2, 3, 4)]
ΔGr = G[NAPA-A(H2O)n (n = 1, 2, 3, 4)] − G[bare NAPA-A] − G[(H2O)n (n = 1, 2, 3, 4)]
ΔCr = C[NAPA-A(H2O)n (n = 1, 2, 3, 4)] − C[bare NAPA-A] − C[(H2O)n (n = 1, 2, 3, 4)]
ΔSr = S[NAPA-A(H2O)n (n = 1, 2, 3, 4)] − S[bare NAPA-A] − S[(H2O)n (n = 1, 2, 3, 4)]
To calculate ΔEr, ΔHr, ΔGr, ΔCr, and ΔSr the structures of reactants and products are very important. The structures of reactants, TS, and products are optimized using DFT/wB97XD/cc-pVTZ computational approach and are presented in Figure 2. The TS structures were confirmed by checking the Hessian matrices that have an imaginary frequency. The DFT-calculated ΔEr, ΔHr, ΔGr, ΔCr, and ΔSr are summarized in Table 2. The values of ΔEr, ΔHr, and ΔSr are negative (−) which means the formation of [NAPA-A(H2O)n (n = 1, 2, 3, 4)] complexes from NAPA-A, [NAPA-A(H2O)n (n = 1, 2, 3)] and H2O are exothermic. In the formation
![]()
Figure 2. DFT-optimized the most stable structure of NAPA-A, transition state structures of [NAPA-A(H2O)n (n = 1, 2, 3, 4)] complexes, and the most stable structure of [NAPA (H2O)n (n = 1, 2, 3, 4)] complexes.
![]()
Table 2. DFT-calculated energy, enthalpy, Gibb’s free energy, change of heat capacity, and entropy change for the reactions of NAPA-A and [NAPA-A(H2O)n (n = 1, 2, 3)] complexes with water.
ΔEr = Reaction energy, ΔHr = Enthalpy of reaction, ΔGr = Gibb’s free energy of reaction, ΔCr = Change of heat capacity and ΔSr = Change of Entropy.
of [NAPA-A(H2O)2] complex, the values of ΔEr, ΔHr, and ΔSr are more negative (−) compared to the other complexes’ formation. Because the most stable conformer of [NAPA-A(H2O)2] complex has a structure that blocks both C5 and C7 backbone interactions by two water molecules which makes the complex more compact than others. The values of ΔEr, ΔHr, ΔGr, ΔCr, and ΔSr for n = 1, 3, and 4 complexes are increased (more negative to more positive) as the total number of vibrational modes increases. The change ΔHr, ΔGr, ΔCr, ΔSr and the reaction energy (ΔEr) for the successive addition of water molecule is depicted in Figure 3.
3.3. The Barrierless Reactions and Kinetics
Calculated energy of transition states, free energy of activation, and rate of the reaction for the reactions of NAPA-A and [NAPA-A(H2O)n (n = 1, 2, 3)] complexes with water are shown in Table 3. A sketch of the calculated potential energy surface (PES) with DFT/wB97XD/cc-pVTZ computational approach for the NAPA-A + (H2O)n (n = 1, 2, 3, 4) → [NAPA-A(H2O)n (n = 1, 2, 3, 4)] are shown in Figure 4 and the structures of reactants, TS and products are presented in Figure 2. Transition state geometries were obtained by QST3 method and were confirmed that Hessian matrices have an imaginary frequency. Transition state complexes have longer hydrogen bond distances than the stable [NAPA-A(H2O)n (n = 1, 2, 3, 4)] complexes (products) (Figure 2).
It is very interesting that the barrier height for the transition states (TS) of [NAPA-A(H2O)n (n = 1, 2, 3, 4)] complexes are computed to lie at 4.41, 4.05, 3.72, and 2.26 kcal/mol below the reactants.
![]()
Figure 3. Graphs representing reaction energy, the change of enthalpy, the change of Gibb’s free energy, the change of heat capacity and the change of entropy for the successive addition of water.
![]()
Table 3. Calculated energy of transition states, free energy of activation, and rate of the reaction for the reactions of NAPA-A and [NAPA-A(H2O)n (n = 1, 2, 3)] complexes with water.
![]()
Figure 4. Potential energy surface (PES) for the reactions of NAPA-A and [NAPA-A(H2O)n (n = 1, 2, 3)] complexes with water having no energy barrier. Relative energies are given in kcal/mol.
Consequently, the barrier height for the production of the stable complex of [NAPA-A(H2O)1] is found to lie 9.29 kcal/mol below the reactants, and for [NAPA-A(H2O)2] complex is computed to lie 0.03 kcal/mol above the reactants. These calculated results confirmed that the products formed global minima [32] . The NAPA-A conformer has both C5 and C7 backbone interaction positions. For the 1st reaction NAPA-A + (H2O)1 → [NAPA-A(H2O)1] complex, where water molecule goes C5 interactions position and formed very stable complex and for the 2nd reaction [NAPA-A(H2O)1] + (H2O)1 → [NAPA-A(H2O)2] complex, where water molecule goes C7 interactions position and formed very stable complex.
As water molecules are strongly bonded via two hydrogen bonds, the value of entropy is more negative compared to the higher cluster and therefore these two reactions are barrierless [33] .
On the other hand, the 3rd reaction [NAPA-A(H2O)2] + (H2O)1 → [NAPA-A(H2O)3] complex and the 4th reaction [NAPA-A(H2O)3] + (H2O)1 → [NAPA-A(H2O)4] complex, the 3rd and 4th water molecule are not directly connected to the backbone of NAPA-A. That means these two water molecules are loosely bound to NAPA-A and hence entropy value is less negative compared to 1st and 2nd complexes. Due to the increase of entropy, the 3rd and 4th reactions are endothermic, and the barrier heights are predicted to lie at 6.30 and 10.54 kcal/mol above the reactants. To understand the kinetic behaviors, Gibb’s free energy of activation (Δ‡G0) as well as the rate of reaction (k) have been calculated for all the reactions. The calculated values of Δ‡G0 are 4.43, 4.28, 3.83, and 5.11 kcal/mol for 1st, 2nd, 3rd, and 4th reactions respectively. Finally, the rate of reactions is 3.490 × 109, 4.514 × 109, 9.688 × 109, and 1.108 × 109 per second for 1st, 2nd, 3rd and 4th reactions respectively. That means all the reactions are very fast but the formation of trihydrated NAPA-A complex from dihydrated NAPA-A and water is faster than other hydrated complex formation reactions.
4. Conclusion
In this work, we studied the thermodynamics and kinetics of N-acetyl-phenylalaninylamide (NAPA) and its microhydration using the explicit solvent model in the gas phase. The most stable conformers of [NAPA-A(H2O)n (n = 0 - 4)] complexes are taken into account for kinetics and thermochemical analysis. The values of H, G, Cv, and S of NAPA-A and its hydrated complexes increase with the successive addition of water molecules due to the increase in the total number of vibrational modes. On the other hand, the values of ΔEr, ΔHr, and ΔGr become more negative to less negative which indicates the sequential addition of water in NAPA-A makes exothermic to endothermic reactions. Transition state geometries were calculated by the QST3 method and the barrier height for the transition states (TS) of [NAPA-A(H2O)n (n = 1, 2, 3, 4)] complexes are predicted to lie 4.41, 4.05, 3.72 and 2.26 kcal/mol below the reactants. According to our calculations, the first two reactions means the formation of [NAPA-A(H2O)1] and [NAPA-A(H2O)2] complexes are barrierless reactions because both water molecules are strongly bonded via two hydrogen bonding in the backbone of NAPA-A peptide. Whereas [NAPA-A(H2O)3] and [NAPA-A(H2O)4] complex formation reactions are endothermic and the barrier heights are calculated to lie at 6.30 and 10.54 kcal/mol respectively above the reactants. To see the kinetic effect, we calculated the rate of reactions for the formation of [NAPA-A(H2O)1], [NAPA-A(H2O)2], [NAPA-A(H2O)3], and [NAPA-A(H2O)4] complexes are 3.490 × 109 s−1, 4.514 × 109 s−1, 9.688 × 109 s−1 and 1.108 × 109 s−1 respectively. This theoretical observation messages us that the microhydration of NAPA-A is very fast and spontaneous.
Data Availability
All generated data is available from the authors under request.