For the consideration of therapeutic peptides from the point of view of medicine, it is necessary to know their molecular properties and their bioactivity. It is our belief that the bioactivity of these peptides is intimately related to their chemical reactivity from a molecular perspective. For this reason, we consider it essential to study the chemical reactivity of natural products that have the potential to become medicines through the tools provided by Computational Chemistry and Molecular Modeling. Probably the most powerful tool currently available to study the chemical reactivity of molecular systems from the point of view of Computational Chemistry and Molecular Modeling is the Conceptual DFT   , also called Chemical Reactivity Theory, which using a series of global and local descriptors allow to predict the interactions between molecules and understand the way in that chemical reactions proceed.
When one is working with Computational Chemistry and Molecular Modeling, there is not a universal methodology that could be applied to the entire spectra of known or unknown molecular systems. As it, also our own devised methodology cannot be considered useful for all the systems. The workaround that researchers in this field have found relies in studying with a particular methodology different but related families of molecules in order to see if the proposed methodology can be applied with the same degree of success to the different group of molecules and that the results obtained with a particular family of molecules are good not by chance or some kind of serendipity. Indeed, the larger the number of studies with different groups of molecules increases the validity of the used methodologies.
Considering that the knowledge of the chemical reactivity is essential for the development of new medicines, we have decided to study in this work the Leucine-Enkephalin which is an endogenous opioid peptide neurotransmitter that is found naturally in the brains of many animals and humans   and that could be the foundation for the design of new therapeutic peptides by relying on a protocol that we have devised recently as a way to validate it with the study of the chemical properties of this important peptide.
Then, rather than employing the usual methodology in Medicinal Chemistry research based on Molecular Docking and QSAR and QSPR approximations, we prefer to rely on that mentioned protocol for our research on peptides. Thus, the objective of this work is to study the chemical reactivity of the Leu-Enkephalin opioid peptide neurotransmitter using the techniques of the Conceptual DFT, determining its global properties (of the molecule as a whole) as well as the local properties that allow to understand and predict active reaction sites, both electrophilic and nucleophilic by considering the recent methodology proposed by us  -  . Likewise, the pKa value of the peptide will be predicted based on a methodology previously developed by us  , the ability of this potentially therapeutic peptide to act as inhibitor of the formation of Advanced Glycation Endproducts (AGEs) will be established according to our previous ideas  , and the descriptors of bioavailability and bioactivity (Bioactivity Scores) will be calculated through different procedures described in the literature. Thus, The knowledge of the values of the global and local descriptors of the molecular reactivity of the Leu-Enkephalin peptide studied could be useful in the development of new drugs based on this compound or some analogs.
2. Computational Methodology
Consistent with our previous work  -  , the computational studies were performed with the Gaussian 09  series of programs that implement density functional methods. The basis set Def2SVP was used in this work for geometry optimization and frequency determination, while the Def2TZVP basis set was used for calculating electronic properties   . All calculations were performed in the presence of water as solvent under the Solvation Model Density (SMD) parameterization of the Integral Equation Formalism-Polarized Continuum Model (IEF-PCM)  .
To calculate the molecular structure and properties of the studied system, we have chosen eight density functionals known to consistently provide satisfactory results for several structural and thermodynamic properties:
CAM-B3LYP Long-range-corrected B3LYP by the CAM method 
LC-ωPBE Long-range-corrected ωPBE density functional 
M11 Range-separated hybrid meta-GGA 
MN12SX Range-separated hybrid nonseparable meta-NGA 
N12SX Range-separated hybrid NGA 
ωB97 Long-range corrected density functional 
ωB97X Long-range corrected density functional 
ωB97XD ωB97X version including empirical dispersion 
In these functionals, GGA denotes the generalized gradient approximation, in which the density functional depends on the up/down spin densities and their reduced gradient, and NGA denotes the nonseparable gradient approximation, which is similar to GGA but also adopts a nonseparable form.
The SMILES notation of the studied compound was fed in the online Molinspiration software from Molinspiration Cheminformatics (www.molins.piration.com) for the calculation of the molecular properties (Log P, Total polar surface area, number of hydrogen bond donors and acceptors, molecular weight, number of atoms, number of rotatable bonds, etc.) and for the prediction of the bioactivity score for different drug targets (GPCR ligands, kinase inhibitors, ion channel modulators, enzymes and nuclear receptors). The bioactivity scores were compared with those obtained through the use of other software like MolSoft from Molsoft L.L.C. (http://molsoft.com/mprop/) and ChemDoodle Version 9.02 from iChemLabs L.L.C. (www.chemdoodle.com).
3. Results and Discussion
3.1. Geometry Optimization and Global Reactivity Descriptors Calculation
The molecular structure of Leu-Enkephalin, which graphical sketch is shown in Figure 1, was preoptimized in the gas phase by considering the DFTBA model
Figure 1. Graphical sketch of the Leu-Enkephalin molecule.
available in Gaussian 09 and then reoptimized using the eight density functionals mentioned in the previous section together with the Def2SVP basis set and the SMD solvent model using water as the solvent. After verifying that each of the structures corresponded to the minimum energy configurations through a frequency calculation analysis, the electronic properties were determined by using the same model chemistry but with the Def2TZVP basis set instead of that used for the geometry optimization.
The analysis of the results obtained in the study aimed at verifying that the KID procedure was fulfilled. On doing it previously, several descriptors associated with the results that the HOMO and LUMO calculations obtained are related with results obtained using the vertical I and A following the ΔSCF procedure. A link exists between the three main descriptors and the simplest conformity to the Koopmans’ theorem by linking with-I, with-A, and their behavior in
describing the HOMO-LUMO gap as ,
, and . Notably, the descriptor consists of an approximation that remains valid only when the HOMO that a radical anion has (the SOMO) shares similarity with the LUMO of the neutral system. Consequently, we decided to design another descriptor ΔSL, to guide in verifying the accuracy of the approximation  -  . The results of this analysis are presented in Table 1.
As can be seen from Table 1, the results for the descriptors show values that are consistent with our previous findings for the case of the melanoidins  -  , that is, only the MN12SX and N12SX density functionals are capable of giving HOMO and LUMO energies that allow verifying the agreement with the approximate Koopmans’ theorem. This is not only true because the values are almost zero, but due to the fact that the ΔSL descriptor, which relates to the difference between the LUMO of the neutral and the HOMO of the anion, is also close to zero. Indeed, these values cannot be exactly equal to zero, but the small differences mean that errors in the prediction of the global reactivity descriptors will be negligible. Moreover, it can be seen from Table 1 that the MN12SX and N12SX density functionals are the only ones that predict negative values for the LUMO energies which will represent positive values of the electron affinity A.
Table 1. Electronic energies of the neutral, positive, and negative molecular systems (in au) of Leu-Enkephalin; the HOMO, LUMO, and SOMO orbital energies (in eV); and the , , and ΔSL descriptors calculated with the eight density functionals and the Def2TZVP basis set using water as the solvent simulated with the SMD parametrization of the IEF-PCM model.
An opposite and incorrect (unphysical) behavior is observed from Table 1 for the other density functionals considered in this work.
By taking into account the KID procedure presented in our previous works together with the finite difference approximation, the global reactivity descriptors can be expressed as:
Electronegativity  
Global Hardness  
Electrodonating Power 
Electroaccepting Power 
Net Electrophilicity 
where and are the energies of the HOMO and LUMO, respectively.
According to our previous discussion, the results for the calculated global reactivity descriptors based on the values of the HOMO and LUMO energies according to the previous definitions will be significative only for the MN12SX and N12SX density functionals. Thus, these results are presented in Table 2.
As expected from the molecular structure of this species, its electrodonating ability is more important that its electroaccepting character. There is no significative differences between the values obtained by using either of the density functionals for the calculation of the global reactivity descriptors. Notwithstanding, after an inspection of Table 1 it can be said that the MN12SX density functional is a somewhat better that the N12SX density functional in verifying the approximate Koopmans behavior. Thus, only the MN12SX density functional will be considered for the remaining of this work.
3.2. Local Reactivity Descriptors Calculation
Applying the same ideas as before, the definitions for the local reactivity descriptors will be:
Nucleophilic Fukui Function  
Electrophilic Fukui Function  
Dual Descriptor  - 
Nucleophilic Parr Function  
Electrophilic Parr Function  
where , , and are the electronic densities at point for a system with , N, and electrons, respectively, and and are related to the atomic spin density (ASD) at the atom of the radical cation or anion of a given molecule, respectively  .
As it has been stated by Martínez-Araya in some recent works    , while the Fukui function is a nice descriptor to understand the local reactivity of the molecules, it can be demonstrated that the Dual Descriptor in its condensed form will perform better for the prediction of the preferred sites for the electrophilic and nucleophilic attacks. For this reason, we have decided to
Table 2. Global reactivity descriptors for the Leu-Enkephalin molecule calculated with the MN12SX and N12SX density functionals with the Def2TZVP basis set and the SMD solvation model using water as the solvent.
Figure 2. (a) Electrophilic Fukui function and (b) Nucleophilic Fukui function for the Leu-Enkephalin molecule.
present the results for the Condensed Dual Descriptor as calculated from either Mulliken Population Analysis (M) or Natural Population Analysis (N) in comparison with the Nucleophilic Parr Function and Electrophilic Parr Function proposed by Domingo et al.   considering atomic spin densities coming from the mentioned Mulliken Population Analysis (M) or from Hirshfeld Population Analysis (H).
The results for the calculation of these local reactivity descriptors for the Leu-Enkephalin molecule are presented in Table 3. It must be noted that we are presenting only the results for those atomic sites where the (which is itself multiplied by 100) are greater than 1. Also, the H atoms are not shown.
As can be seen from Table 3, the local reactivity descriptors calculated from the different formulations are able to recognize the nucleophilic and electrophilic sites for chemical reactivity with great accuracy. Moreover, there is an impressive agreement between the results coming from the Condensed Dual Descriptor and the Nucleophilic and Electrophilic Parr Functions and which means that their use in this and future works related to the study of therapeutic peptides will be a warranty of success.
3.3. Determination of pKa Value of the Peptide
We have recently presented a study of the computational prediction of the pKas of small peptides through Conceptual DFT descriptors  . In that work, we concluded that the relationship pKa = 16.3088 - 0.8268η could be a valuable starting point for the prediction of the pKa of larger peptides of interest for the development of AGE inhibitors.
Table 3. Local reactivity descriptors for the Leu-Enkephalin molecule calculated with the MN12SX density functional with the Def2TZVP basis set and the SMD solvation model using water as the solvent: Condensed Dual Descriptor , Nucleophilic Parr Function and electrophilic parr function ; M stands for mulliken population analysis, N corresponds to natural population analysis and h means hirshfeld population analysis.
Thus, we have now applied the mentioned relationship to the calculation of the pKa of the Leu-Enkephalin molecule giving a result of 11.90. This result could be of interest when designing pharmaceutical drugs starting from these peptide allowing to explain the mechanisms of action and the drug delivery procedures.
3.4. Quantification of the AGEs Inhibitor Ability
The Maillard reaction between a reducing carbonyl and the amino group of a peptide or protein leads to the formation of a Schiff base which through a series of steps renders different molecules known as Advanced Glycation Endproducts or AGEs. It is believed that the presence of these AGEs is one of the main reasons for the developing of some diseases like Diabetes, Alzheimer and Parkinson.
Among several strategies that have been considered for the prevention of the formation of AGEs, it is worth to mention the use of compounds presenting amino groups in their structure capable of interacting with the reducing carbonyl of carbohydrates and being competitive with the amino acids, peptides and proteins present in our body. Many compounds have been devised as drugs to achieve this goal and to name a few, we can include Pyridoxamine, Aminoguanidine, Carnosine, Metformin, Pioglitazone and Tenilsetam.
It can be proposed that peptides having amino and amido groups could be thought as potential therapeutic drugs for preventing the formation of AGEs. In a previous work, we have studied the ability of a group of proposed molecules to act as inhibitors of the formation of AGEs by quantifying their behavior in terms of Conceptual DFT reactivity descriptors  . It was concluded that the key factor in the study of the chemical reactivity of the potential AGEs inhibitors was on their nucleophilic character and although there are several definitions of nucleophilicity  , our results suggested that the inverse of the net electrophilicity could be a good definition for the nucleophilicity N. On the basis of the mentioned analysis, we were able to find some qualitative trends for the studied molecular systems.
In this work, we will extend this correlation to the Leu-Enkephalin peptide in order to see if it can be considered as a precursor of therapeutic drugs for the inhibition of the formation of AGEs. As the model chemistry employed in both works is the same, the comparison is straightforward:
Aminoguanidine > Metformin > Carnosine > Leu-Enkephalin
> Tenilsetam > Pyridoxamine > Pioglitazone
This qualitative trend is representative of the known pharmacological properties of the studied AGEs inhibitors   and it can be seen that Leu-Enkephalin possesses better AGEs Inhibitor Ability that Tenilsetam and Pyridoxamine, but lower than Aminoguanidine, Metformin and Carnosine.
3.5. Bioactivity Scores
When considering a given molecular system as a potential therapeutic drug, it is customary to check if the considered species follows the Lipinsky Rule of Five which is used to predict whether a compound has or not has a drug-like character  . The molecular properties related to the drug-like character were calculated with the aid of the MolSoft and Molinspiration software and are presented in Table 4 where miLogP represents the octanol/water partition coefficient, TPSA is the molecular polar surface area, natoms is the number of atom of the molecule, nON and nOHNH are the number of hydrogen bond acceptors and hydrogen bond donors respectively, nviol is the number of violations of the
Table 4. Molecular properties of the Leu-Enkephalin peptide calculated to verify the Lipinsky Rule of Five.
Lipinsky Rule of Five, nrotb is the number of rotatable bonds, volume is the molecular volume, and MW is the molecular weight of the studied system.
However, what the Lipinsky Rule of Five really measures is the oral bioavailability of a potential drug because this is desired property for a molecule having drug-like character. Indeed, this criteria cannot be applied to peptides, even when they are small, as we can see from Table 4, due to the inherent molecular weight and number of hydrogen bonds.
In a more recent work, Martin  have developed what she called “A Bioavailability Score” (ABS) for avoiding these problems. The rule for the ABS established that the Bioavailability Score for neutral organic molecules must be 0.55 if they pass the Lipinsky Rule of Five and 0.170 if they fail. The ABS value for all the Leu-Enkephalin peptide considered in this work have been calculated by using the ChemDoodle software and the result was equal to 0.170.
Then, a different approach was followed by considering similarity searches in the chemical space of compounds with structures that can be compared to those that are being studied and with known pharmacological properties.
As has been mentioned in the Settings and Computational Methods section, this task can be accomplished using the online Molinspiration software for the prediction of the bioactivity score for different drug targets (GPCR ligands, kinase inhibitors, ion channel modulators, enzymes and nuclear receptors). The results are named Bioactivity Scores and the values for the Leu-Enkephalin are presented in Table 5.
These bioactivity scores for organic molecules can be interpreted as active (when the bioactivity score > 0), moderately active (when the bioactivity score lies between −5.0 and 0.0) and inactive (when the bioactivity score < −5.0) as it is summarized in Table 6. The Leu-Enkephalin peptide was found to be bioactive towards the Protease Inhibitor and the GPCR Ligand considered in the study.
Table 5. Bioactivity scores of the Leu-Enkephalin molecule calculated on the basis of GPCR ligand, ion channel modulator, nuclear receptor ligand, kinase inhibitor, protease inhibitor and enzyme inhibitor intteractions.
Table 6. Factors for evaluation of bioactivity.
In this paper, we have presented the results of a study of the chemical reactivity of the Leucine-Enkephalin opioid neurotransmitter peptide based on the Conceptual DFT as a tool to explain the molecular interactions.
The knowledge of the values of the global and local descriptors of the molecular reactivity of the Leu-Enkephalin peptide studied could be useful in the development of new drugs based on this compound or some analogs.
In a similar manner, the pKa value for the potentially therapeutic peptide has been predicted by resorting to the value of the chemical hardness following a previously proposed methodology and the information that resulted would be helpful in understanding not only the chemical reactivity but other important properties like the water solubility.
A point of special interest has been the quantification of the ability of the peptide to act as an inhibitor in the formation of AGEs and this could be of importance for the design of medicines for fighting diseases like Diabetes, Alzheimer or Parkinson.
Finally, the molecular properties related to bioavailability have been predicted using different methodologies already described in the literature, and the descriptors used for the quantification of the bioactivity allowed characterizing the studied peptide as being bioactive towards the Protease Inhibitor and the GPCR Ligand considered in the study. In this way, our work is on the same line as previous studies on the bioactivity of organic molecules   with the novelty of being applied to a peptide of great therapeutic importance.
This work has been partially supported by CIMAV, SC and Consejo Nacional de Ciencia y Tecnologa (CONACYT, Mexico) through Grant 219566-2014 for Basic Science Research. Daniel Glossman-Mitnik conducted this work while a Visiting Lecturer at the University of the Balearic Islands from which support is gratefully acknowledged. Norma Flores-Holguín and Daniel Glossman-Mitnik are researchers of CIMAV and CONACYT. This work was cofunded by the Ministerio de Economía y Competitividad (MINECO) and the European Fund for Regional Development (FEDER) (CTQ2014-55835-R).
 Frau, J. and Glossman-Mitnik, D. (2018) Conceptual DFT Study of the Local Chemical Reactivity of the Dilysyldipyrrolones A and B Intermediate Melanoidins. Theoretical Chemistry Accounts, 137, 1210.
 Frau, J. and Glossman-Mitnik, D. (2018) Computational Study of the Chemical Reactivity of the Blue-M1 Intermediate Melanoidin. Computational and Theoretical Chemistry, 1134, 22-29.
 Frau, J. and Glossman-Mitnik, D. (2018) Chemical Reactivity Theory Applied to the Calculation of the Local Reactivity Descriptors of a Colored Maillard Reaction Product. Chemical Science International Journal, 22, 1-14.
 Frau, J., Flores-Holgun, N. and Glossman-Mitnik, D. (2018) Chemical Reactivity Properties, pKa Values, AGEs Inhibitor Abilities and Bioactivity Scores of the Mirabamides A-H Peptides of Marine Origin Studied by Means of Conceptual DFT. Marine Drugs, 16, 302-319.
 Frau, J., Hernández-Haro, N. and Glossman-Mitnik, D. (2017) Computational Prediction of the pKas of Small Peptides through Conceptual DFT Descriptors. Chemical Physics Letters, 671, 138-141.
 Frisch, M.J., Trucks, G.W., Schlegel, H.B., Scuseria, G.E., Robb, M.A., Cheeseman, J.R., Scalmani, G., Barone, V., Mennucci, B., Petersson, G.A., Nakatsuji, H., Caricato, M., Li, X., Hratchian, H.P., Izmaylov, A.F., Bloino, J., Zheng, G., Sonnenberg, J.L., Hada, M., Ehara, M., Toyota, K., Fukuda, R., Hasegawa, J., Ishida, M., Nakajima, T., Honda, Y., Kitao, O., Nakai, H., Vreven, T., Montgomery, J.A., Peralta, J.E., Ogliaro, F., Bearpark, M., Heyd, J.J., Brothers, E., Kudin, K.N., Staroverov, V.N., Kobayashi, R., Normand, J., Raghavachari, K., Rendell, A., Burant, J.C., Iyengar, S.S., Tomasi, J., Cossi, M., Rega, N., Millam, J.M., Klene, M., Knox, J.E., Cross, J.B., Bakken, V., Adamo, C., Jaramillo, J., Gomperts, R., Stratmann, R.E., Yazyev, O., Austin, A.J., Cammi, R., Pomelli, C., Ochterski, J.W., Martin, R.L., Morokuma, K., Zakrzewski, V.G., Voth, G.A., Salvador, P., Dannenberg, J.J., Dapprich, S., Daniels, A.D., Farkas, O., Foresman, J.B., Ortiz, J.V., Cioslowski, J. and Fox, D.J. (2016) Gaussian 09 Revision E.01. Gaussian Inc., Wallingford.
 Weigend, F. and Ahlrichs, R. (2005) Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Physical Chemistry Chemical Physics, 7, 3297-3305.
 Marenich, A., Cramer, C. and Truhlar, D. (2009) Universal Solvation Model Based on Solute Electron Density and a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. Journal of Physical Chemistry B, 113, 6378-6396.
 Yanai, T., Tew, D.P. and Handy, N.C. (2004) A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chemical Physics Letters, 393, 51-57.
 Henderson, T.M., Izmaylov, A.F., Scalmani, G. and Scuseria, G.E. (2009) Can Short-Range Hybrids Describe Long-Range-Dependent Properties? The Journal of Chemical Physics, 131, Article ID: 044108.
 Peverati, R. and Truhlar, D.G. (2011) Improving the Accuracy of Hybrid Meta-GGA Density Functionals by Range Separation. The Journal of Physical Chemistry Letters, 2, 2810-2817.
 Peverati, R. and Truhlar, D.G. (2012) Screened-Exchange Density Functionals with Broad Accuracy for Chemistry and Solid-State Physics. Physical Chemistry Chemical Physics, 14, 16187-16191.
 Chai, J. and Head-Gordon, M. (2008) Systematic Optimization of Long-Range Corrected Hybrid Density Functionals. Journal of Chemical Physics, 128, Article ID: 084106.
 Chai, J. and Head-Gordon, M. (2008) Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections. Physical Chemistry Chemical Physics, 10, 6615-6620.
 Martnez-Araya, J.I. (2012) Revisiting Caffeate’s Capabilities as a Complexation Agent to Silver Cation in Mining Processes by Means of the Dual Descriptor—A Conceptual DFT Approach. Journal of Molecular Modeling, 18, 4299-4307.
 Martnez-Araya, J.I. (2012) Explaining Reaction Mechanisms Using the Dual Descriptor: A Complementary Tool to the Molecular Electrostatic Potential. Journal of Molecular Modeling, 19, 2715-2722.
 Martínez-Araya, J.I. (2015) Why Is the Dual Descriptor a More Accurate Local Reactivity Descriptor than Fukui Functions? Journal of Mathematical Chemistry, 53, 451-465.
 Domingo, L.R., Pérez, P. and Sáez, J. (2013) Understanding the Local Reactivity in Polar Organic Reactions through Electrophilic and Nucleophilic Parr Functions. RSC Advances, 3, 1486-1494.
 Chamorro, E., Pérez, P. and Domingo, L.R. (2013) On the Nature of Parr Functions to Predict the Most Reactive Sites along Organic Polar Reactions. Chemical Physics Letters, 582, 141-143.
 Domingo, L.R., Ríos-Gutiérrez, M. and Pérez, P. (2016) Applications of the Conceptual Density Functional Theory Indices to Organic Chemistry Reactivity. Molecules, 21, 748.
 Rajasekaran, S., Rao, G. and Zonunsiami (2017) Molecular Properties and Bio-Activity Score of 2[2-(4-chlorophenyl)-4-oxoquinazolin-3(4H)-yl]amino-N (substitutedphenyl) Acetamides. Journal of Pharmaceutical Research, 16, 95.
 Erylmaz, S. (2018) The Theoretical Investigation of Global Reactivity Descriptors, NLO Behaviours and Bioactivity Scores of Some Norbornadiene Derivatives. Sakarya University Journal of Science, 22, 1.