The starting molecular structures of the two plant cyclopeptides to be studied, shown in Figure 1, lwere obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov), which is an open chemistry database.

As a first step, the SMILES (Simplified Molecular Input Line Entry Specification) notation of every studied compound, which was obtained by accessing ChemDoodle 11.3.0 software, was fed into the online program Chemicalize, a software developed by ChemAxon (http://www.chemaxon.com), which was used for naming, molecular finger prints, structure generation and the prediction of several properties related to Chemoinformatics (http://chemicalize.com/) (accessed March 2021).

The similarity searches in the chemical space of compounds with molecular structures similar to those that are being studied was accomplished using the online available Molinspiration software from Molinspiration Cheminoformatics (https://www.molinspiration.com/) (accessed, March 2021) which was used for the prediction of the bioactivity scores for different drug targets.

A Webtool named SwissTargetPrediction for efficient prediction of protein targets of small molecules has been considered for the determination of the potential bioactivity of the two terrestrial plant cyclopeptides considered in this study (Daina et al., 2019). The associated website allows the estimation of the most probable macromolecular targets of a small molecule, assumed as bioactive. During the process of development of a new therapeutic drug, it is of the most importance to adcquire a knowledge of the fate of the pharmacokinetics, that is, the fate of a compound in the organism. This is usually performed by through individual indices that are called Absorption, Distribution, Metabolism, Excretion and Toxicity (ADMET) parameters. These parameters are generally obtained using computer models as an alternative to the experimental procedures for their determination. In this research, some ADME parameters were estimated with the aid of Chemicalize and the online available SwissADME software (Daina et al., 2017).Additional information about the Pharmacokinetics parameters and the ADMET properties were obtained by resorting to pkCSM (Pires et al., 2015), a software for the prediction of small-molecule pharmacokinetic properties using SMILES (https://biosig.unimelb.edu.au/pkcsm/) (accessed, March 2021).

The goodness of a given density functional can be determined through a comparison of the results that it renders with the experimental values or with the results that can be obtained by means of high-level calculations. However, the lack of experimental results for the molecular systems under study or the large size of the molecules that made computationally impractical the use of some accurate methodologies. Kohn-Sham (KS) methodology includes the determination of the molecular energy, the electronic density and the orbital energies of a given system, related to the frontier orbitals including the Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO) (Young, 2001; Lewars, 2003; Cramer 2004; Jensen, 2007). This methodology is convenient when thinking of quantitative qualities related with Conceptual DFT descriptors (Parr and Yang, 1989; Chermette, 1999; Geerlings et al., 2003; Toro-Labbé, 2007; Chattaraj, 2009, Geerlings et al., 2020; Chakraborty and Chattaraj, 2021). Range-separated (RS) exchange-correlation density functionals are of extraordinary concern in Kohn-Sham DFT calculations (Iikura et al., 2001; Yanai et al., 2004; Heyd and Scuseria, 2004; Chai and Head-Gordon, 2008; Stein et al., 2009a; Stein et al., 2009b; Stein et al., 2010; Karolewski et al., 2011; Kuritz et al., 2011; Ansbacher et al., 2012; Kronik et al., 2012; Stein et al., 2012). A methodology called KID (Koopmans in DFT) has been established by our research group (Frau et al., 2018; Flores-Holguín et al., 2019a; Frau et al., 2019; Flores-Holguín et al., 2019c; Flores-Holguín et al., 2020; Flores-Holguín et al., 2020a; Flores-Holguín et al., 2020b; Flores-Holguín et al., 2021), t for the validation of a given density functional in terms of its internal coherence. Several descriptors associated with the results of the HOMO and LUMO calculations are related to the results obtained from the estimation of the vertical I and A following the ΔSCF procedure, where SCF refers to the Self-Consistent Field technique. It has been demonstrated that there is a relationship between the KID descriptors and the Koopmans’ theorem or the Ionization Energy theorem, which is its equivalent within the Generalized Kohn-Sham (GKS) version of DFT, by connecting ϵ _H to -I, ϵ _L to -A, and their actions by defining the HOMO – LUMO gap J _I = |ϵ _H + E _gs (N − 1) − E _gs (N)|, J _A = |ϵ _L + E _gs (N) − E _gs (N + 1)|, and J H L = J I 2 + J A 2 . It should be noticed that the J _A descriptor consists of an approximation which is only valid if the HOMO of the radical anion (the SOMO) resembles the LUMO of the neutral system. For this reason, another descriptor ΔSL has been designed by our research group (Frau et al., 2018; Flores-Holguín et al., 2019a; Frau et al., 2019; Flores-Holguín et al., 2019c; Flores-Holguín et al., 2020; Flores-Holguín et al. 2020a; Flores-Holguín et al., 2020b; Flores-Holguín et al., 2021), to help in the verification of the accuracy of the approximation.

Taking into account the KID methodology considered in the previous research being integrated into the finite difference approximation (Frau et al., 2018; Flores-Holguín et al., 2019a; Frau et al., 2019; Flores-Holguín et al., 2019c; Flores-Holguín et al., 2020; Flores-Holguín et al., 2020a; Flores-Holguín et al., 2020b; Flores-Holguín et al., 2021), the following definitions can be used for the global descriptors that help in the understanding of the chemical reactivity of the molecular systems (Parr and Yang, 1989; Chermette, 1999; Geerlings et al., 2003; Gázquez et al., 2007; Chattaraj et al., 2009; Chakraborty and Chattaraj, 2021): Electronegativity as χ = − 1 2 ( I + A ) ≈ 1 2 ( ϵ L + ϵ H ) , Global Hardness as η = (I − A) ≈ (ϵ _L −ϵ _H ), Electrophilicity as ω = μ 2 / 2 η = ( I + A ) 2 / 4 ( I − A ) ≈ ( ϵ L + ϵ H ) 2 / 4 ( ϵ L − ϵ H ) , Electrodonating Power as ω ⁻ = ( 3 I + A ) 2 / 16 ( I − A ) ≈ ( 3 ϵ H + ϵ L ) 2 / 16 η , Electroaccepting Power as ω ⁺ = ( I + 3 A ) 2 / 16 ( I − A ) ≈ ( ϵ H + 3 ϵ L ) 2 / 16 η and Net Electrophilicity as Δω ^± = ω ⁺ − (−ω ⁻) = ω ⁺ + ω ⁻, being ϵ _H and ϵ _L the HOMO and LUMO energies associated with each of the peptides considered in this work. It is worth to mention that for the global indices the chemical power is directly related with the electronic density as well as the corresponding Hohenberg-Kohn functional (Putz, 2011).

As a complement of these global reactivity descriptors that arise from Conceptual DFT (Parr and Yang, 1989; Chermette, 1999; Geerlings et al., 2003; Gázquez et al., 2007; Chattaraj et al., 2009; Chakraborty and Chattaraj, 2021), Domingo and his collaborators (Domingo et al., 2008; Jaramillo et al., 2008; Domingo and Sáez, 2009; Domingo and Perez, 2011; Domingo et al., 2016) have proposed a Nucleophilicity index N through the consideration of the HOMO energy obtained through the KS scheme with an arbitrary shift of the origin taking the molecule of tetracyanoethylene (TCE) as a reference.

The determination of the conformers of the molecules considered in the current study was performed by resorting to MarvinView 17.15 available from ChemAxon (http://www.chemaxon.com) by doing Molecular Mechanics calculations through the overall MMFF94 force field (Halgren, 1996a, Halgren, 1996b, Halgren, 1999; Halgren and Nachbar, 1996; Halgren, 1996c). This was followed by a geometry optimization and frequency calculation by means of the Density Functional Tight Binding (DFTBA) methodology (Frisch et al., 2016). This last step was required for the verification of the absence of imaginary frequencies as a check for the stability of the optimized structures as being a minimum in the energy landscape. The electronic properties and the chemical reactivity descriptors of the studied molecules involved the use of MN12SX/Def2TZVP/H2O model chemistry (Weigend and Ahlrichs, 2005; Weigend, 2006; Peverati and Truhlar, 2012) on the optimized molecular structures due to is ability in the verification of the “Koopmans in DFT” (KID) protocol (Frau and Glossman-Mitnik, 2018a; Frau and Glossman-Mitnik, 2018b; Frau and Glossman-Mitnik, 2018c; Frau and Glossman-Mitnik, 2018d; Frau and Glossman-Mitnik, 2018e; Frau and Glossman-Mitnik, 2018f; Flores-Holguín et al., 2019a; Flores-Holguín et al., 2019b; Flores-Holguín et al., 2019d; Frau et al., 2019; Flores-Holguín et al., 2019c, Flores-Holguín et al., 2020; Flores-Holguín et al., 2020a; Flores-Holguín et al., 2020b; Flores-Holguín et al., 2021) using Gaussian 16 (Frisch et al., 2016) and the SMD model for the simulation of the solvent (Marenich et al., 2009). This model chemistry considers the MN12SX screened-exchange density functional (Peverati and Truhlar, 2012) together with the Def2TZVP basis set (Weigend and Ahlrichs, 2005; Weigend, 2006) and in all cases the charge of the molecules is equal to zero while the radical anion and cation have been considered in the doublet spin state.

The names, identifiers, molecular fingerprints and basic properties of the two Pashinintides A and B plant cyclopeptides are presented in Table 1, while their geometrical and structural properties are displayed in Table 2.

This information could it be of interest for future SAR studies based on these and other peptides as well as for potential derivatives designed for therapeutical purposes using Peptidomimetics.

The effect of the geometrical and structural properties on the bioavailability of the Pashinintides A and B presented in Table 2 can be better visualized considering the Bioavailability Radars displayed in Figure 2.

It can be appreciated that the two more important properties that could prevent the use of the Pashinintides A and B as therapeutic drugs are their molecular size and their polar character. Although these cyclic peptides violate some of the limits traditionally considered to be important for oral bioavailability of drug-like small molecules, it can be expected that the reduced flexibility could ease oral absorption (Nielsen et al., 2017). However, it must be remarked that these ideal values are based on the Lipinski’s Rule of Five (Lipinski et al., 2001), which is not always applicable to peptides (Zhang and Wilkinson, 2007; Doak et al., 2014; Santos et al., 2016; Nielsen et al., 2017; Sable et al., 2017).

The Bioactivity Scores for the Pashinintides A and B are shown in Table 3.

It can be seen from the results on Table 3 that while the bioactivity of Pashinintide B towards the different targets is considered to be low, for Pashinintide A, its interactions as a GPCR Ligand and a Protease Inhibitors could be of importance for its consideration as a potential therapeutic drug. The same conclusion can be extracted by checking visually the predicted biological targets for these plant cyclopeptides shown in Figure 3.

An ADMET study is the assessment of pharmacokinetics of a drug which stands for Absorption, Distribution, Metabolism, Excretion and Toxicity. The prediction of the fate of a drug and the effects caused by a drug inside the body, such as how much drug is absorbed if administered orally and how much is absorbed in the gastrointestinal tract, is an indispensable part of drug discovery. In a similar way, if the absorption is poor, its distribution and metabolism would be affected, which can lead to causing neurotoxicity and nephrotoxicity. Ultimately, the study is to understand the disposition of a drug molecule within an organism. Thus, ADMET study is one of the most essential parts of computational drug design.

The optimized molecular structures of the Pashinintides A and B are displayed in Figure 4.

Although the Koopmans-complaint behavior of the MN12SX density functional has been proved previously for the case of marine peptides (Frau et al., 2018; Flores-Holguín et al., 2019a; Frau et al., 2019; Flores-Holguín et al., 2019c; Flores-Holguín et al., 2020; Flores-Holguín et al., 2020a; Flores-Holguín et al., 2020b; Flores-Holguín et al., 2021), we are now performing a further validation for the plant cyclopeptides considered in the present study. This determination has been done by resorting to the in-house developed CDFT software tool and the resulting values are shown in Table 9.

As can be seen from the values presented in Table 9, the KID descriptors are all very close to zero meaning that the chosen MN12SX density functional displays a Koopmans-complaint behavior. This in agreement with our previous studies on peptides (Frau et al., 2018; Flores-Holguín et al., 2019a; Frau et al., 2019; Flores-Holguín et al., 2019c; Flores-Holguín et al., 2020; Flores-Holguín et al., 2020a; Flores-Holguín et al., 2020b; Flores-Holguín et al., 2021), thus justifying the adequacy of the MN12SX/Def2TZVP/H2O model chemistry for the purpose of this research.

The results for the global reactivity indices were estimated by making use of the mentioned CDFT tool and are presented in Table 10.

As the global hardness η can be regarded as a direct measure of the deformation of the electron density and of the chemical reactivity being related to the HOMO-LUMO gap, it can be seen that Pashininitide A will be slightly more reactive than the other cyclopeptide. The electrodonating ability ω ⁻ is more important that its electroaccepting power ω ⁺ for both cyclopeptides because of their molecular structures. However, after a comparison of the values of ω ⁻ and ω ⁺ for each molecule, it can be concluded that there are not important differences between them. The electrophilicity ω index encompasses the balance between the tendency to acquire an extra amount of electron density by an electrophile and the resistance of a molecule to exchange electron density with the environment Domingo et al. (2016). By studying the electrophilicities of a series of reagents involved in Diels-Alder reactions (Domingo et al., 2002a; Domingo and Sáez, 2009; Pérez et al., 2003), an electrophilicity ω scale for the classification of organic molecules as strong, moderate or marginal electrophiles was proposed being ω > 1.5 eV for the first case, 0.8 < ω < 1.5 eV for the second case and ω < 0.8 eV for the last case (Domingo et al., 2002a; Domingo and Sáez, 2009; Pérez et al., 2003). By inspection of Table 10, it can be said that both peptides considered in this study may be regarded as moderate electrophiles. Notwithstanding, the overall chemical reactivity is about the same for both cyclopeptides. This information could be of interest for future studies on the potential therapeutic ability of these compounds.

Besides global reactivity descriptors, their local counterparts have been developed to get an idea of the differences in chemical reactivity between the atoms within the molecule. Among these local reactivity descriptors are the Fukui functions (Parr and Yang, 1989; Chermette, 1999; Geerlings et al., 2003) and the Dual Descriptor (Toro-Labbé, 2007; Morell et al., 2005, Morell et al., 2006; Martínez-Araya, 2012a; Martínez-Araya, 2012b; Martínez-Araya, 2015), which have been defined as: Nucleophilic Fukui Function (NFF) = f ⁺(r) = ρ _N+1(r) − ρ _N (r), Electrophilic Fukui Function (EFF) = f ⁻(r) = ρ _N (r) − ρ _N−1(r), and Dual Descriptor (DD) = Δf(r) = ∂ f ( r ) / ∂ N υ ( r ) , relating the electronic densities of the neutral, positive and negative species.

The NFF, f ⁺(r), is associated with the sites within a molecular system which are prone to nucleophilic attacks while the EFF, f ⁺(r), describes those sites that are more susceptible to electrophilic attacks. Although the NFF and the EFF have been used successfully for the identification of reactive sites, the Dual Descriptor Δf(r) or DD, can describe unambiguously nucleophilic and electrophilic sites within a molecule (Martínez-Araya (2015)). A graphical representation of the DD for the Pashinintides A and B cyclopeptides is displayed in Figure 5 showing the zones where DD > 0 and DD < 0.

Although these graphical representations allowed to distinguish the regions within the molecules where the Dual Descriptor will be greater or smaller than zero, it can appreciated that there is some overlap between them. Thus, for a better estimation of these reactivity areas it is worth to determine the values of the Condensed Dual Descriptor (Δf _k ) (Morell et al., 2008; Frau and Glossman-Mitnik, 2018d) over the all the atoms (excluding H) in comparison with the condensed versions of the Electrophilicity, that is, the Condensed Electrophilicity (ω _k ) (Domingo et al., 2002b), and of the Nucleophilicity, being the Condensed Nucleophilicity N _k (Pérez et al., 2009). The resulting values are displayed in Tables 11 and 12 for the Pashinintides A and B, respectively.

Even if every atom within the peptides cannot be graphically individualized due to the large size of the molecules, it is clear from the results in Tables 11 and 12 about which are the sites prone to electrophilic and nucleophilic attacks based on the agreement between the values for the Condensed Electrophilicity ω _k with the positive values of the Condensed Dual Descriptor Δf _k , for one side, and with the values of the Condensed Nucleophilicity N _k and the negative results for the Condensed Dual Descriptor Δf _k .

From Table 11 it can be seen that for Pashinintide A the maximum values for the Condensed Electrophilicity ω _k (shown in bold) are located over the C (27) and O (4) which correlate well with the maximum positive results for the Condensed Dual Descriptor Δf _k over those atoms. The same situation is found for the case of the Condensed Nucleophilicity N _k , whose maximum values over the C (29) and C(35) correlate with the maximum negative values Condensed Dual Descriptor Δf _k localized on those atoms. The closeness between C(27) and C (29) explains the overlap between the two regions within the graphical representation of the Dual Descriptor.

For the case of Pashinintide B, it can be appreciated from Table 12 that the maximum values are located over the O (1), C (24), O (5) and C (37) (in that order) correlating in agreement with the greatest positive results for the Condensed Dual Descriptor Δf _k , while for the Condensed Nucleophilicity N _k the order of reactivity will be C (27) > C (36) > C (40) > C (45) > N (13) being the same as for those derived from the Condensed Dual Descriptor Δf _k . As for the case of the other peptide, the partial overlapping between the different reactive areas could be attributed to nearness between C (24) and C (27).