Computational Investigation of Conformational Properties of Short Azapeptides: Insights from DFT Study and NBO Analysis

1. Introduction

Peptides are promising tools with biological and pharmaceutical applications ¹. However, their low metabolic stability because of enzymatic degradation, lack of receptor selectivity, and many other drawbacks made it necessary to search for alternatives with the same function as natural peptides but with better properties ². These alternatives are called peptidomimetics, in which some structural modifications involve peptides’ backbone or side chain ^{2, 3}. Among peptidomimetics, we focus on azapeptides, a change in the backbone of a peptide by substituting the α-carbon with a nitrogen atom ^4-6. This modification confers azapeptides more resistance to physiological degradation and long duration of action ⁷. Many studies revealed the importance of incorporating an azaamino acid in the peptide’s sequence to enhance the stability and bioavailability of drug candidates ^8-10. In addition, synthetic methods ^{4, 9, 11} and biological studies ^{5, 8, 9, 12} of azapeptides have been studied. However, the structural behaviors of azaamino acid residue in peptides remain limited.

The conformational preferences of azaamino acid residues have been studied theoretically ^{13, 14, 15} and experimentally ^{16, 17, 18, 19}. The results showed that the preferred dihedral angles of the backbone of these azaamino acid residues are in the range of ϕ= ±90°±30°, ψ= 0° or ±180°±30°, which is appeared to be the polyproline II (β_P), and β-turn motif (α_R, α_L, δ_R, or δ_L conformation) ^{14, 20}. The nomenclature was denoted by Karplus et al. ²¹. Our groups and others examined the effect of azaamino acid residue on the structures in peptides. For examples, the azaglycine adopted the polyproline II (β_P) in azaGly-Pro-Hyp (hydroxyproline) of the collagen ¹⁷; the δ_L conformer (ϕ= 90°±30°, ψ= 0°±30°) in Ac-Phe-azaGly-NH₂, forming the βII-turn structure; the α_R or α_L conformer (ϕ= ±90°±30°, ψ= 0°) in Ac-Aib-azaGly-NH₂, adopting βI or βI’-turn ²². Meanwhile, incorporating azaamino acids destabilized the β sheet secondary structure compared to the parent peptide [20]. The results showed that the conformational preference of azaamino acid residue in peptides should be sequence dependence. However, it is still largely unknown of these questions: why azaamino acid residue adopts a specific conformation in the context dependence, the influence of the side chain of azaamino acid residue on its backbone structure, the neighboring effect on the conformational behaviors of azapeptides. As a part of the efforts to answer these questions, we previously examined the conformational properties of tripeptides, Ac-Pro-azaXaa-NHMe [Xaa= Asn (1), Asp(2), Ala(3)] using the DFT method, showing that the azaamino acid in this sequence, regardless of the side chain, adopts βⅡ (βⅡ’) turn structure [15]. And we also suggested the importance of the side-chain and main-chain or main-chain and main-chain hydrogen bonds found in these model peptides in the stability of the azapeptides. In this work, we extended our ongoing efforts in understanding conformational properties of short azapeptides containing three azaamino acid residues, Ac-azaXaa-Pro-NHMe [Xaa = Asn (4), Asp(5), and Ala(6)] (Figure 1). While previous theoretical studies focused on the dipeptide models ¹⁴, we wondered the position effect of azaamino acid residue in tripeptides. We used the popular DFT functionals, B3LYP and B3LYP-D3, to compare our previous work ¹⁵ in the gas phase and water. The solvation effect was calculated with the SMD method ²³. We also investigated the influence of the cis-trans isomerization amide bonds on the conformational preferences of model peptides.

2. Computational Methods

All calculations were performed using the Gaussian 09/16 programs package ²⁴. The starting geometries of the three compounds were from the previous work ¹⁵ and were fully optimized in the gas phase at the B3LYP/6-311++G(d, p) level of theory. These geometries combine all possible intramolecular hydrogen bonds formed in the molecule between acceptor groups C=O or N^α and the N-H donor group. The orientation of amide bonds in model azapeptides, which is either cis (ω= 0°) or trans (ω= ±180°) amide bond, was also assessed. All conformations were found to be true local minima, as was indicated by the absence of imaginary frequencies. These conformations were then reoptimized using the B3LYP-D3 ²⁵ in the gas phase using the same basis set. The solvent effect was examined using the solvation model density (SMD) ²³ with and without the dispersion-correction B3LYP functional in water. Intramolecular hydrogen bond interactions found in the molecules in the gas phase (Figure 2) were also inspected using Natural Bond Orbital (NBO) theory at the two levels of theory using NBO 3.1 program implemented in Gaussian 16 package ²⁴. The strength of the hydrogen bond in azapeptides models was estimated using the second-order perturbation analysis (E(2) value) ²⁶. The NBO overlaps were drawn using the Multiwfn program ²⁷.

3. Results and Discussion

The backbone (ϕ and ψ) dihedral angles of the resulting structures for the three azapeptides 4-6 are denoted by the nomenclature given by Karplus ²¹ (Supporting information Tables S1–S12). Besides, the side chain of azaAsn and azaAsp are marked g+(gauche+), g-(gauche-), s+(skew+), and s-(skew-) based on the values of χ₁’ ²⁸ Relative energies of the most stable conformations of the three model compounds 4-5 calculated at the B3LYP/6-311++G(d,p) and B3LYP-D3/6-311++G(d,p) levels of theory in the gas phase and water are listed in Table 1 and Table 2.

3.1. Cis-Trans Isomerization for Minimum Energy Conformations

Ac-azaAsn-Pro-NHME (4): The conformer tttANP-1 (α_Rδ_R(g⁺)), corresponding to the type of βI-turn, and has a trans-trans orientation of the amide bond preceding both residues, was found as the lowest energy minimum for azaAsn model in the gas phase calculated with the B3LYP method. However, the same conformation is the second energy minimum at the B3LYP-D3 functional with a relative energy of 0.94 kcal/mol. Interestingly, the lowest energy minimum conformer of azapeptide 4 calculated at the B3LYP functional is found in the X-ray structure containing azaAsn-Pro sequence, showing the presence of a β-Ⅰ turn with all-trans amide bonds ^{18, 19}. In inclusion of dispersion correction, B3LYP-D3 functional predicts tctANP-2 as the lowest energy minimum containing a cis amide bond between azaAsn-Pro sequence, which is also predicted to be the most stable conformer in water (Table 2). The second lowest energy minimum in solution calculated with B3LYP functional is ttANP-12 (δ_Lβ_P(s⁺)) with a relative energy of 0.94 kcal/mol. This relative energy value was also found for the second energy minimum calculated with B3LYP-D3 functional. Another thing to highlight is the trans amide bond preceding both Pro and azaAsn found for the five lowest energy minima calculated at the B3LYP function in the isolated state, except for the second energy minimum in which the amide bond preceding azaAsn is cis-oriented (Table 1). However, most of the lowest energy minimum conformers calculated with B3LYP-D3 functional in the gas phase are cis-oriented amide bonds preceding either the two residues or at least one. In the solution state, the five lowest energy minima with B3LYP-D3 functional are all trans amide bonds except for the first and the fifth conformers, which have a cis amide bond preceding azaAsn residue. The theoretical results suggest that the azaAsn-Pro sequence may also be involved in cis-trans isomerization in water.

Ac-azaAsp-Pro-NHMe (5): The results of Ac-azaAsp-Pro-NHMe are similar to those of the AzaAsn model. B3LYP functional predict tttADP-1 (α_Rδ_R(g⁺)) with all-trans amide bonds as the lowest energy minimum in the gas phase. The second most stable conformation, tttADP-2 (α_Rγ’(g⁺)) (ΔE= 1.00 kcal/mol), is also oriented to all-trans amide bonds. However, B3LYP-D3 functional predicts the conformer tctADP-3 as the lowest energy minimum and the conformer tctADP-4 as the second energy minimum. The results imply that the dispersion interaction may be essential in stabilizing the cis amide bond in the azaAsp-Pro sequence. Contrary to the isolated state, both SMD/B3LYP and SMD/B3LYP functionals predict the conformer tttADP-9 (δ_Lβ_P(s⁺)) as the lowest energy conformer (Table 2). Similarly, the second most stable conformation in water with the two methods is tttADP-7, with all-trans amide bonds. However, the third lowest energy minimum tctADP-3 (βδ_R(g⁺)) (ΔE=1.07 kcal/mol) has a cis orientation of amide bond preceding Pro residue. The results also imply that the azaAsp-Pro sequence may be involved in cis-trans isomerization in water.

Ac-azaAla-Pro-NHMe (6): The results of azapeptide 6 significantly differ from those of azapeptides 4 and 5. Both B3LYP and B3LYP-D3 functionals predict that the conformer tttAAP-1(α_Rδ_R), related to the β-Ⅰ turn structure, is the lowest energy minimum in the gas and water. In addition, the B3LYP function predicts that the four most stable conformers are trans-trans amide bonds in the gas phase. However, the third and the fourth energy minimum have trans-cis amide bonds with B3LYP-D3 functional. The two most stable conformations calculated by SMD/B3LYP in water have all trans amide bonds, while. In contrast,he most stable conformation calculated by SMD/B3LYP-D3 was found to be all trans amide followed by two trans-cis-trans or trans-cis-trans conformers (Table 2). The results imply that the side chain of azaAsn and azaAsp may play a vital role in cis-trans isomer. At the same time,n while the conformer of azalea residue containing a cis amide bond is found in the higher energy minima.

3.2. Hydrogen Bonds in Ac-azaXaa-Pro-NHMe

Optimized structures of Ac-azaXaa-Pro-NHMe are presented with the denotation of hydrogen bonds (Supporting Information Tables S1–S12 and Figures S1–S12. NBO orbital interactions corresponding to hydrogen bonds are depicted in the Supporting information Tables S13–S17, calculated with the B3LYP method. We used the E(2) value to describe the strength of the hydrogen bond ²⁶.

The most stable conformation in the gas phase at the B3LYP method is tttANP-1, which is stabilized with two intramolecular hydrogen bonds (HBs) closing C10 and C6 pseudo cycles between carbonyl and N-H groups C²=O¹….H²-N⁴ and C⁸=O⁴….H¹-N¹, respectively (Figure 3). The C₁₀ pseudo cycle is a HB interaction related to the β-Ⅰ turn which is found in X-ray structure of Z-azaAsn(Me)-Pro-NHiPr ¹⁹. The HB is made up of three interactions: π_CO → σ^*_NH of 1.43 kcal/mol, n_O→σ^*_NH of 1.98 kcal/mol, and n’O → σ^*_NH of 0.81 kcal/mol (Figure 3). The sum ΣE_HB of the individual E(2) stabilization energies of the three H.B. interactions is 4.22 kcal/mol. It is not the strongest interaction, yet tttANP-1 is the most stable conformation. Moreover, the strongest interaction is that related to the C8 pseudo cycle in tctANP-7 conformation in which the sum of the individual E(2) stabilization energies, associated with π_CO, n_O, and n’_O / σ^*_NH overlaps involved in the C.O….HN HB is 9.34 kcal/mol. The same result was found with the B3LYP-D3 method (Supporting information Table S16) since the tctANP-7 conformation has the highest value of the sum of E(2) (9.98 kcal/mol) for the C8 HB pseudo cycle. However, the most stable conformation, tttANP-1, has a sum value of 5.39 kcal/mol related to the C10 pseudo cycle. We notice a slight increase in the values of the sum of the individual E(2) energies calculated with B3LYP-D3 compared to those found with B3LYP.

The opposite conclusions are drawn in water; the SMD/B3LYP method predicted the tttANP-12 (δ_Lβ_P) conformer to be the lowest energy minimum in water. This conformer is stabilized with one H.B. closing a C8 pseudo cycle (Supporting information Figure S7) between the C=O of the main chain and the N-H of the azaAsn side chain. The length of the hydrogen bond is calculated to be 1.92Å. Unlike the gas phase, results in water showed that the second most stable conformation is the conformer tttANP-1, stabilized with an H.B. interaction closing a C10 pseudo cycle, and the H.B. length is equal to 2.30Å. We observed the absence of the C6 pseudo cycle found in the gas phase for the tttANP-1 conformer. Also, the HBs length becomes longer than that in the gas phase.

Unlike the B3LYP functional, the B3LYP-D3 functional predicted the conformer tctANP-2 to be the lowest energy minimum in gas and solution states. This conformer is stabilized with two intramolecular HBs in the gas phase. One is a backbone-side chain HB closing a C6 pseudo cycle, and a backbone HB closing a C5 ring, with the length of HBs estimated to be 2.26Å for both interactions. The C6 ring was not observed in the solution for this same conformer. However, the length of the C5 HB was slightly larger than that found in the gas phase (2.33Å). Our result indicates that the side chain of azaAsn in solution may be involved in an intermolecular interaction with water.

The results of compound (5) showed few similarities with those of compound (4). The most stable conformer in the gas phase with the B3LYP functional is tttADP-1, with two H.B.s interactions closing C₁₀ and C₆ rings. The C₁₀ HB is related to the βⅠ turn motif for which the sum of the individual E(2) stabilization energies (Supporting information Table S14), associated with π_CO, n_O and n’_O / σ^*_NH overlaps (Figure 4) involved in the CO….HN HB is 3.71 kcal/mol. This turn motif was also found in the X-ray structure of Z-azaAsp(OEt)-Pro-NhiPr¹⁹. Comparing these results with those observed with B3LYP-D3 functional, we see clearly that the dispersion-corrected B3LYP-D3 functional might not predict the lowest energy minimum; since it is tctADP-3 conformer with one H.B. closing a C5 pseudo cycle (Supporting information Figure S5) with a sum of E(2) is 1.22 kcal/mol related to the n’N→ σ^*_NH NBO overlap of the N.H…..N backbone H.B. ²⁹(Figure 4).

The tttADP-9 conformer was found as the most stable conformation in water with SMD/B3LYP and SMD/B3LYP-D3 methods. This conformer is stabilized with one H.B. interaction, a C=O…H-O backbone-side chain HB. We also notice that the length of this type of H.B. is slightly larger with the dispersion-corrected functional (1.75Å) than that found with the traditional B3LYP functional (1.64Å).

The most stable conformation of compound (3) in the gas phase with and without dispersion functional is ttAAP-1, which is stabilized with an H.B. closing a C10 pseudo cycle related to βⅠ turn (Supporting information Figures S3 and S6). The length of such bond is predicted to be 2.16 and 2.08 Å with B3LYP and B3LYP-D₃ functional, respectively. The HB length has slightly decreased from B3LYP to B3LYP-D3_, contrary to azaAsn and azaAsp models. The NBO overlaps of the C10 HB are shown in Figure 5. The sum of E(2) associated with π_CO, n_O, and n’_O / σ^*_NH overlaps calculated with B3LYP functional is 3.62 kcal/mol (Supporting Information Table S15). While that calculated with B3LYP-D3 functional is 4.15 kcal/mol (Supporting information Table S18), which is only made of two π_CO and n_O / σ^*_NH overlaps.

Results of the azaAla model in the water show significant similarities with those of the gas phase. The ttAAP-1 conformer is the lowest energy minimum and stabilized with a C10 backbone C=O…H-N HB. The length of this H.B. is found to be 2.32Å with B3LYP functional and 2.14Å with B3LYP-D3 functional, which is the same finding as in the gas phase.

3.3. Asx Turn Vs βⅠ-tun in Ac-azaAsx-Pro-NHMe (x= n or p)

Hydrogen-bonded Asx turns are similar to β turns since both have a C10 pseudo cycle. However, in Asx turns, it is the side chain’s carbonyl group of Asn or Asp at the i+1 position, which is hydrogen bonded to NH of the backbone of residue at the i+3 position. Thus, the nomenclature of Asx-turns is the same as β-turns. Four categories were found: type-Ⅰ, type-Ⅱ, type-Ⅰ′, and type-Ⅱ′ Asx turn (Table 3). This nomenclature is based on the dihedral angles ϕ and ψ of residues i+1 and i+2. Since Asx turns are mimicry of β turns, the dihedral angles ϕ_e and ψ_e of Asx turn (Figure 6) should resemble those of i+1 in β turns, while ϕ and ψ of residue i+1 of Asx turns should resemble those of residue i+2 in β turns ^{30, 31}.

Furthermore, the azaAsx-Pro chain contains another C10 pseudo-cycle similar to an Asx turn, which is related to the H.B. involving the N^δH or O^δ-H of the side chain of Asn or Asp, respectively, and carbonyl group two residues ahead which will be noted (C10SC). Results of the azaAsn model compound in the gas phase calculated with the B3LYP functional showed that all three conformations adopting the βⅠ-turn motif are more stable than those adopting the Asx turn. However, the most stable Asx turn conformation tttANP-9 is only 2.89 kcal/mol higher in energy, which means there is a competition between Asx and βⅠ turn secondary structures. The C10 C=O….H-N backbone side chain HB with the sum of E(2) being 5.76 kcal/mol which is associated with n_O and n’_O / σ^*_NH overlaps (Supporting information Table S13). On the other hand, the most stable Asx turn conformation with B3LYP-D3 functional is cctANP-17, which has a sum of E(2) value of 1.88 kcal/mol related to n_O and n’_O / σ^*_NH overlaps (Supporting information Table S16). Moreover, the majority of βⅠ-turn conformations are more stable than the Asx turn, and the cctANP-16 conformation related to the C10SC interaction is more stable than the Asx turn. The sum of the E(2) value (3.05 kcal/mol) related to the C10 HB of this conformer is higher than that of cccANP-17, which means a stronger HB bond interaction. The Asx turn was not reported in the azaAsn-Pro sequence yet. Instead, the N^δ-H of the side chain of azaAsn is bounded to N^α of the same residue ^{18, 19}. Moreover, based on the dihedral angles of a typical Asx turn (Table 3), the most stable Asx turn conformation tttANP-9 (ΔE =2.89 kcal/mol for B3LYP; ΔE=4.52 kcal/mol for B3LYP-D3) is of type-Ⅰ Asx turn. This Asx turn was also found in Asn-Pro studied with NMR and IR spectroscopies ³⁰. On the contrary, a study of a database of 500 proteins found that the frequency of occurrence of Asx turn is type-Ⅱ′, while the most common type is βI turn ³⁰. These findiA recent study of Asx also supported these findingsn-Ala and Asn-Gly in crystallized protein structures in the PDB ³¹, revealing that the most common Asx turn in these two sequences are type-Ⅱ′ and type-Ⅱ, respectively. In sum, our results suggest that azaAsn and azaAsp residues in models 4 and 5 could have different folding patterns than natural peptides and the,y could also enhance the formation of a type-Ⅰ Asx turn instead of type-Ⅱ′ mostly found in peptides. Even though the Asx turn is more stable than the C10S,C which is similar to the Asx turn, most conthe formers with C10SC are more stable than those adopting the Asx turn, especially with B3LYP-D3 functional.

In solution, results showed few similarities regarding the most stable Asx turn conformation with B3LYP functional, which was found to be the conformer tttANP-9. However, the tctANP-14 conformer related to the C10SC turn also has the same relative energy as tttANP-9 (∆E=2.95 kcal/mol). The C10SC turn was found to be more stable than Asx turn in solution with the B3LYP-D3 functional since the tctANP-14 conformer adopting the C10SC turn has a relative energy of 2.07 kcal/mol and the tttANP-9 conformer related to the Asx turn is 4.20 kcal/mol energy. The result indicates that the side chain of azaAsn in the azaAsn-Pro sequence acts as a proton donor in water.

Similarly, the βⅠ turn structure in the azaAsp-Pro sequence with B3LYP functional is more stable than Asx turn in the isolated state. This is observed when comparing tttADP-1 and tttADP-8 conformers. The latter conformer related to the type-Ⅰ Asx turn is 3.70 kcal/mol higher in energy than the first one, which adopts the βⅠ turn. It is also noticed that the Asx turn in this sequence is less stable than that found for the azaAsn-Pro sequence. In addition, conformers cctADP-15 and tctADP-16 adopt an undefined type of Asx turn stabilized with the C10 cycle. Like the azaAsn model compound, the Asx turn in the azaAsp model compound is more stable than the C10SC motif as we compare the most stable Asx turn conformer tttADP-8 and the most stable C10SC conformer cctADP-19 (ΔE=7.40 Kcal/mol). Also, most Asx turn conformations are energetically more favorable than C10SC conformations. Like the B3LYP functional, the B3LYP-D3 functional predicted that the Asx turn is more stable than the C10SC turn in the isolated state. The B3LYP-D3 functional showed that the tctADP-16 conformer is the most stable Asx turn with a relative energy of 3.83 kcal/mol, and the sum of E(2) value is equal to 2.47 kcal/mol and associated with n_O and n’_O / σ^*_NH overlaps.

Similar results were found in the solution state since the βⅠ turn is more stable than C10SC and Asx turn conformations. However, the C10SC is more stable than Asx turn with B3LYP functional, conformation tcADP-22 is the most stable C10SC turn conformation with a relative energy value of 5.15 kcal/mol while conformation cctADP-15 adopting the Asx turn has a relative energy value of 9.29 kcal/mol. The B3LYP-D3 functional reveals that the C10SC turn is more stable than the Asx turn, with the conformer tctADP-22 being the most stable C10SC conformer in water with a relative energy of 3.45 kcal/mol.

4. Conclusions

In summary, conformational preferences of three azapeptides models Ac-azaXaa-Pro-NHMe [Xaa= Asn, Asp, Ala] were studied using the B3LYP/6-311++G(d,p) and B3LYP-D3/6-311++G(d,p) levels of theory in isolated and water to examine the effect of side chain’s nature and position effect on the folding patterns of azapeptides. Our result demonstrates that the position of azaamino acid residue largely affects the folding patterns of these compounds. Another finding is the orientation of amide bonds preceding azaamino acid residue and Pro. The most stable conformations of all three compounds in gas and solution states calculated with B3LYP functional are all trans amide bonds. Moreover, the stability order of the three secondary structures stabilized with a 10-membered HB in azaAsn and azaAsp models is βⅠ ˃ C10S ˃Asx for azaAsn-Pro and azaAsp-Pro sequences in water calculated by SMD/B3LYP and SMD/B3LYP-D3, and for athe zaAsp-Pro model in the gas phase also calculated by two methods. In contrast, the stability order for the azaAsn-Pro in the gas phase calculated with the two functionals is βⅠ ˃ Asx ˃ C10SC. Note the hat azaXaa-Pro sequence prefers βI-turn whereas Pro-azaXaa sequence prefers βII-turn structure ¹⁵. Our computational investigation of the conformational properties of short azapeptides demonstrates that the nature of the side chain of azapeptides could largely influence the folding patterns of these components.