Comparative Study of the Chemical Reactivity of Helical Peptide Models for Protein Glycation

Non-enzymatic glycation of proteins has been implicated as an important cause of the complications associated with diabetes and Alzheimer disease. It is well known that glycation involves the reactivity of, primarily, the ε -amino group of the lysines present in the protein. The immediate chemical environment of an amino group modulates the glycation reaction. In this work, several model helical peptides for protein glycation has been studied by resorting to QM:MM calculations through the ONIOM methodology. Some Conceptual DFT descriptors have been calculated that allowed the comparison of the chemical reactivity between the different model peptides in terms of the position of the Lys group and other spatially proximate amino acid residues.


Introduction
The nonenzymatic reaction between reducing carbohydrates and the amino groups of amino acids, peptide and proteins is called glycation and through a series of chain reaction leads to the formation of Advanced Glycation Endproducts (AGEs).These molecules have been related to some diseases like diabetes, Alzheimer and Parkinson.Nonenzymatic glycation of proteins has been implicated as an important cause of the complications associated with diabetes and Alzheimer disease.It is well known that glycation involves the reactivity of, primarily, the ε-amino group of the lysines present in the protein.
The immediate chemical environment of an amino group modulates the glycation reaction.Some years ago, Venkatraman et al [1] studied several helical peptide models in order to elucidate the proximity effect in the catalysis of the Amadori rearrangement.Their conclusions were similar to those of Howard et al [2] and Povey et al [3] who found that the glycation reaction affected the helicity, charge and other properties of the resulting molecules.
Therefore, we consider that it will be of practical interest to develop a theoretical and computational tool that could be an aid for the prediction of the preferred glycation sites for the different possible conformational structures of model peptides because they could be an aid in the development of AGE inhibitors.
Thus, the objective of this work is to predict the preferred glycation sites of the model helical peptides proposed by Venkatraman et al [1]

Theoretical Background
As this work is part of an ongoing project, the theoretical background is analog to that presented in previous research [7]- [13] and will be shown here for the sake of completeness.As it has been described in detail before, [7]- [13] within the conceptual framework of DFT, [5] [14] the chemical potential µ is defined as:

(
) ( ) where χ is the electronegativity, while the global chemical hardness η is: ω + = + − .It follows that a larger ( ω + ) value corresponds to a better capability of accepting charge, whereas a smaller value of ( ω − ) makes it a better electron donor.In order to compare ( ω + ) with ( ω − ), the following defini- tion of net electrophilicity has been proposed: that is, the electroaccepting power relative to the electrodonating power.

Settings and Computational Methods
Following the lines of our previous work [7]- [13], the molecular structures of the studied compounds were constructed by starting with the readily available MOL structures for the amino acids (ChemSpider: www.chemspider.com,PubChem: https://pubchem.ncbi.nlm.nih.gov/).The most stable conformers were obtained by means of the Avogadro 1.2.0 program [18] [19] through a random sampling with Molecular Mechanics techniques and a consideration of all the torsional angles through the general AMBER force field [20].
The energies of the neutral, positive and negative peptides were obtained within the framework of QM:MM calculations performed through the ONIOM method [21] in the presence of water as a solvent, by doing IEF-PCM computations according to the SMD solvation model [22].For the QM region, the MN12SX density functional [23] in connection with the Def2TZVP basis set was chosen [24] [25].The MN12SX density functional is a range-separated hybrid nonseparable meta-NGA [23].The Amber force field [26] was considered for the MM part.All the ONIOM calculations were performed as implemented in the Gaussian 09 [27] series of programs.

Results and Discussion
The following peptides proposed by Venkatraman et al [1] were studied as mentioned before: From the cluster analysis, some members were selected for further analysis: 3 for KD4, 4 for KD2, 3 for KH4, 3 for RKD4 and 1 for K14.
Within Conceptual DFT, the Fukui function is defined in terms of the derivative of ( ) r ρ with respect to N [5]: reflects the ability of a molecular site to accept or donate electrons.High values of ( ) f r are related to a high reactivity at point r [5].
By applying a finite difference approximation to the previous expression, two definitions of Fukui functions depending on total electronic densities are obtained: ( ) ( ) , are the electronic densities at point r for the system with Morell et al. [28]- [34] have proposed a local reactivity descriptor (LRD) which is called the dual descriptor (DD) ( ) ( ) ( ) The dual descriptor can be condensed over the atomic sites: when 0 k f ∆ > the process is driven by a Table 1.The ionization potentials I and electron affinities A (in eV), and global electronegativity χ , otal chemical hardness η , global electrophilicity ω , electrodonating power, ( ω − ), electroaccepting power ( ω + ), and net electrophilicity ω ± ∆ of the optimized conformers of the KD4, KD2, KH4, RKD4 and K14 model peptides calculated with the ONIOM method (MN12SX/Def2TZVP:Amber) using water as solvent simulated with the SMD parametrization of the IEF-PCM model.
The first thing that can be observed from Table 2 is that the results for the descriptors through calculations based on the Mulliken Population Analysis (MPA) or Hirshfeld Population Analysis (HPA) are roughly the same and that there are not qualitatively differences between them.
However, it can be observed from Table 1 that the chemical reactivity will be  different for the different conformational structures of the model peptides.The amino groups of the Lys residues in the model helical peptides make these systems to behave as electrodonating molecules in the glycation reaction.Thus, KD4-3 and KD4-2 will be more reactive than KD4-1.The same trend is observed for the chemical hardness η , the global electrophilicity ω and the net electrophilicity ω ± ∆ .If we now turn to Table 2, it can be seen that the values of the condensed dual descriptor from Table 1 and Table 2 for the other model helical peptides.
Indeed, the sites for glycation will be also the preferred sites for protonation.
Thus, the pKa's of the Lys residues will be also dependent on the conformational structures of the model peptides.Moreover, the trend in the pKa of the different conformers may be predicted in terms of the local hypersoftness (LHS), which is a local reactivity descriptor that has been defined so that it permits to measure local reactivitiesaccording to the molecular size [30] [32].The working equation is expressed as follows: ( ) ( ) ⋅ , where S stands for the global softness [4] [38] [39].As the local hypersoftness can be condensed over the atomic sites, the condensed local hypersoftness is simply computed as . The procedure is explained as follows: ( ) f r ∆ is expressed in atomic units, meanwhile S is measured in mili eV raised to the power of −1, however before performing the multiplication, the mili factor is turned back into 10 −3 and then S is raised to the power of 2; the resulting value uses the unit mili eV raised to the power of −2, meaning m(eV −2 ); the parenthesis is put in order to make clear that the prefix mili is not raised to the power of -2.By considering the most reactive conformer of each model peptide, the following trend for the pKa's of the Lys residues can be found: KD4 ≈ KD2 > KH4 > RKD4 > K14.That is, although the exact value of the pKa cannot be predicted, a qualitative information about their relative ordering can be established by considering the global and local Conceptual DFT descriptors.

Conclusion
The preferred glycation sites of several model helical peptides have been by resorting to the calculation of some Conceptual DFT descriptors like the Fukui function indexes, the condensed dual descriptor ( ) f r ∆ and the electrophilic and nucleophilic Parr functions.As the preferred glycation sites will be the same as those for the protonation reaction, the results of the calculations will allow us to qualitatively predict the pKa's of the different lysine residues on the light of the obtained values for the Conceptual DFT descriptors [4] [5] [6].
, where I and A are the vertical ionization potential and elec- tron affinity, respectively.The electrophilicity index ω has been defined as: electrodonating ( ω − ) and electroaccepting ( ω + ) powers have been defined as:

1 N + , N and 1 N
− electrons, respectively.The first one, associated to reactivity for a nucleophilic attack so that it measures the intramolecular reactivity at the site r towards a nucleophilic reagent.The second one, associated to reactivity for an electrophilic attack so that this function measures the intramolecular reactivity at the site r toward an electrophilic reagent[14].
attacks) which are related to the atomic spin density (ASD) at the r atom of the radical cation or anion of a given molecule, respectively.The ASD over each atom of the radical cation and radical anion of the molecule gives the local nucleophilic k P − and electrophilic k P + Parr functions of the neutral molecule [37].The condensed Fukui functions k f − , condensed dual descriptors k f ∆ and electrophilic Parr functions k P − over the N atoms of the amino groups of the Lys residues of the optimized conformers of the KD4, KD2, KH4, RKD4 and K14 model peptides calculated with the ONIOM method (MN12SX/Def2TZVP: Amber) using water as solvent simulated with the SMD parametrization of the IEF-PCM model are presented in Table and the electrophilic Parr functions k P − over the N atoms of the amino groups of the Lys residues are larger for the KD4-3 and KD4-2 model peptides than for the KD4-1 system.Therefore, the global and local Conceptual DFT descriptors predict the same behavior and can discriminate between the different conformational structures.The same conclusions may be extracted from the values of the Conceptual DFT descriptors established by resorting to the calculation of some Conceptual DFT descriptors like the Fukui function indexes, the condensed dual descriptor Parr functions.The results were obtained within the framework of QM:MM calculations performed through the ONIOM method in the presence of water as a solvent.The pKa's of the different Lys residues could be qualitatively predicted on the light of the obtained values for the Conceptual DFT descriptors.Therefore, the techniques presented in this work can be regarded as a practical computational tool for the prediction of the preferred glycation sites of peptides and proteins and for a qualitative description of the pKa's associated to the ionizable side-chain groups.

Table 2 .
The condensed Fukui functions k f − , condensed dual descriptors − over the N atoms of the amino groups of the Lys residues of the optimized conformers of the KD4, KD2, KH4, RKD4 and K14 model peptides calculated with the ONIOM method (MN12SX/Def2TZVP:Amber) using water as solvent simulated with the SMD parametrization of the IEF-PCM model.MPA: Mulliken Population Analysis, HPA: Hirshfeld Population Analysis.