Conceptual DFT-Based Computational Peptidology of Marine Natural Compounds: Discodermins A–H

A methodology based on the concepts that arise from Density Functional Theory named Conceptual Density Functional Theory (CDFT) was chosen for the calculation of some global and local reactivity descriptors of the Discodermins A–H family of marine peptides through the consideration of the KID (Koopmans in DFT) technique that was successfully used in previous studies of this kind of molecular systems. The determination of active sites of the studied molecules for different kinds of reactivities was achieved by resorting to some CDFT-based descriptors like the Fukui functions as well as the Parr functions derived from Molecular Electron Density Theory (MEDT). A few properties identified with their ability to behave as a drug and the bioactivity of the peptides considered in this examination were acquired by depending on a homology model by studying the correlation with the known bioactivity of related molecules in their interaction with various biological receptors. With the further object of analyzing their bioactivity, some parameters of usefulness for future QSAR studies, their predicted biological targets, and the ADME (Absorption, Distribution, Metabolism, and Excretion) parameters related to the Discodermins A–H pharmacokinetics are also reported.


Introduction
Several natural resources in the sea may generate molecules which could be a lead toward the design of new medical drugs. It is for this reason that many recent studies have been undertaken towards new products where the knowledge of marine natural products is used as an important source of information [1]. Among the systems that can be obtained from natural products with marine origin are peptides, which are molecules with a size ranging between that of proteins and amino acids [2]. Peptides are sources of nitrogen and amino acids which are related to numerous potential physiological functions. Bioactive peptides can be protein fragments that acquire functionality when liberated from the parent protein. The first activity assigned to a peptide was neurotoxicity; however, at present, they are associated with other functions, such as cardiotonic, antiviral and antitumor, cardiotoxic and antimicrobial activity. These functions, in addition to their excellent binding properties, low off-target toxicity, and high stability, make peptides promising molecules for the development of new therapeutics. Approximately 60% of described natural products belong to the peptide family. Peptides are present in many marine species, and the extensive research that has been conducted on them has shown that they most often found in sponges, as occurred with polyketides, peptides can be also produced by commensal or symbiotic bacteria or fungi.

Theoretical Background and Computational Details
Kohn-Sham (KS) methodology includes the estimation of the molecular energy and density of a given system, as well as the orbital energies, explicitly connected with the frontier orbitals including the Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO) [20][21][22][23]. This methodology is convenient when thinking of quantitative qualities related to Conceptual DFT descriptors. At present, the use of range-separated (RS) exchange-correlation density functionals in Kohn-Sham DFT is of extraordinary concern [24][25][26][27]. It is essential to the development of these density functionals to think about the partitioning of the exchange and the −1 12 operator into long-and short-ranged parts partnering with a range-separation ω parameter that controls the rate at which long-range behavior is obtained. The estimation of ω can either be fixed or "tuned" by use of a molecule-by-molecule procedure by adhering to some tuning principles. The ideal tuning methodology depends on having the KS HOMO energy related to vertical ionization potential (IP) which is an estimation of the energy difference, E(N-1)-E(N). For the instance of the Generalized KS theory appropriate to an N-electron molecular system, we should have -IP(N) = H (N), which can be considered to be the DFT counterpart of the well-known Koopmans' theorem. In reality, this is valid only for the exact density functional. For the situation where for pragmatic reasons, we have to consider an approximated density functional, there will be possibly some critical distinction between -IP(N) and H (N), Thus, ideal tuning involves setting up ω, system-specific range-separation parameter, by a nonempirical approach and having a RSE useful density functional [28][29][30][31][32][33][34][35]. Indeed, even with the absence of an equal methodology that can be used for the correlation of the electron affinity (EA) combined with the energy of the LUMO, it can be drawn that H ((N + 1) = −EA(N) is conceivable. This makes the acquisition of the optimized ω value easier, provided that the differences between L (N) and H (N + 1) are small. Through this, the forecast of the Conceptual DFT descriptors expectation is upgraded for a given optimized density functional. The concurrent prescription is dubbed as the "KID procedure" (for Koopmans in DFT), this being in reference to the relationship it has with Koopmans' theorem that has been previously quoted [14][15][16][17][18].
Following the methodology considered in our previous studies [14][15][16][17][18], we have done the computational determinations by using the Gaussian 09 series of programs [36] for the implementation of the density functional needed for the development of this work. The Def2SVP basis set [37,38] was chosen for the geometry optimizations and in the verification that the optimized structures corresponded to the minimal ones through the calculation of the associated frequencies. As a larger basis set is usually needed for the calculation and analysis of the electronic properties, the Def2TZVP basis set was chosen [37,38] as a constituent of our standard methodology. According to our previous studies on this kind of molecular systems, water was selected as the solvent through the Solvation Model Density (SMD) parameterization of the Integral Equation Formalism-Polarized Continuum Model (IEF-PCM) [39] for all the DFT calculations. For the determination of the molecular structures and their associated electronic properties of the studied peptides, the MN12SX density functional was chosen because it is already well known that it is capable of giving very good results for several structural and thermodynamic properties [40]. The resulting model chemistry, MN12SX/Def2TZVP/H20, was then considered model chemistry due to the excellent behavior that our previous research has demonstrated owing to the fact that the MN12SX behaves as a Koopmans-complaining density functional which is very useful for obtaining accurate HOMO and LUMO orbital energies, thus avoiding the determination of the energies of the cationic and anionic systems for which convergence is usually hard to obtain for the somewhat large molecules as peptides are [14][15][16][17][18].

Results and Discussion
The starting molecular structures of the Discodermins A-H peptides to be studied were obtained from ChemSpider (http://www.chemspider.com/), which is a website that acts as an online free resource of chemical information related to physical and biological properties, interactive spectra and literature references. In order to get a glimpse of the potential therapeutic properties of the considered peptides, Simplified Molecular Input Line Entry Specification (SMILES) notations for the molecular systems under consideration were fed into the online Molinspiration software from Molinspiration Cheminformatics (Slovensky Grob, Slovak Republic) allowing for the estimation of several molecular properties that are known to be related to druggability and that could be useful in QSAR studies [41,42] that are commonplace in the process of drug design and development. The results of this determination are presented in Table 1: A further and complementary step can be performed by resorting to SwissADME [43] which is an online free tool that allows the evaluation of drug-likeness through a graphical representation of the properties of interest called Bioavailability Radar which is obtained for every peptide by with the aid of its SMILES representation and where the pink area exhibits the zone with the optimal range for a particular property, as shown in Figure 2: Another information that can be obtained from the precedent study is that related to the pharmacokinetic properties of the potential therapeutic peptides, i.e., how the living organism will interact with the drugs since they are delivered to the body up to when they or their metabolites are finally excreted. This information is collectively known as ADME properties, and this important information in the process of drug discovery is presented in Table 2 for the Discodermins A-H marine peptides.
Indeed, the efficacy of a given therapeutic drug will be highly dependent on the way it interacts with a given receptor in the cells, which in turn with the diseases that the medicines will be designed to fight. This can be predicted beforehand by resorting to a homology procedure using databases of known molecules of similar and related structures to the ones under study. This is an important tool in the drug discovery process and has been used in this work to analyze the Discodermins A-H family of peptides by resorting to the free online SwissTargetPrediction software [44] with the results for the predicted biological targets presented in Figure 3. All the presented results for the potential therapeutic activity of the Discodermins A-H family of peptides are an indication that the chemical reactivity of these molecules is an area worth to be explored. As we did in the past [14][15][16][17][18]. A methodology called Conceptual DFT-based Computational Peptidology developed in our research group will be considered now for the calculation and analysis of the chemical reactivity properties of these interesting molecules. The starting point consists of a search for the most stable conformers of each peptide on the basis of the molecular structures taken from the ChemSpider webpage by resorting to the MarvinView 17.15 program (ChemAxon, Budapest, Hungary) relying on the overall MMFF94 force field [45][46][47][48][49].
The most stable conformer for each peptide obtained through the described procedure were then subjected to a geometry optimization in the gas phase by considering the Density Functional Tight-Binding Approximation (DFTBA) model that is accessible in Gaussian 09 [36]. As it was mentioned in the Theoretical Background and Computational Details section, the resultant structures were subsequently reoptimized by considering the MN12SX/Def2SVP/H2O model chemistry that have proven to be adequate for this purpose [14][15][16][17][18]. Upon verification that every optimized structure corresponded to a minimum in the energy potential curve by employing the frequency-calculation analysis technique, the electronic properties were calculated by resorting to a similar model chemistry but considering the Def2TZVP basis set instead of the Def2SVP one because it has been demonstrated [14][15][16][17][18] that this is a better choice for the prediction and analysis of the HOMO and LUMO orbitals and of the chemical reactivity properties derived from them.
A graphical display of the tridimensional optimized molecular structures of the Discodermins A-H peptides obtained with the aid of the UCSF Chimera Visualization System [50] are presented in Figure 4.
It is usually assumed that the goodness of a given density functional can be estimated by comparing the results that it gives with the experimental values that are trying to be reproduced or with the results that can be obtained through post Hartree-Fock calculations like MP2, MP4 or CCSD. However, this is not always possible due to the lack of experimental results for the molecular systems that are being studied or the large size of the molecules that keep some accurate methodologies to be computationally practical. For this reason, we have developed a protocol named KID (Koopmans in DFT) [14][15][16][17][18], which is an attempt to validate a given density functional in terms of its internal coherence. 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, where SCF refers to the Self-Consistent Field technique. It has been shown that there is a connection between the key descriptors and the simplest conformity to the theorem of Koopmans 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 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 [14][15][16][17][18], to help in the verification of the accuracy of the approximation. Although the Koopmans'-complaining behavior of the MN12SX density functional has been proved previously for the case of peptides [14][15][16][17][18], we think that it is worth performing a further validation for the case of the molecules considered in the present study. This determination has been achieved by making use of the in-house developed CDFT software tool and the results of this analysis are shown in Table 3.  As can be seen from the results in Table 3, the values for the KID descriptors are all very close to zero which is an indication that the chosen MN12SX density functional behaves as a Koopmans-complaining one and that for that reason the MN12SX/Def2TZVP/H2O is a model chemistry that has been further demonstrated very adequate for the purpose of this research.
Taking into account the KID methodology considered in the previous research being integrated into the finite difference approximation [14][15][16][17][18], the following definitions can be used for the global descriptors that help in the understanding of the chemical reactivity of the molecular systems [5][6][7]51,52]: Electroaccepting Power Net Electrophilicity being H and L the HOMO and LUMO energies associated with each of the peptides considered in this work. It is worth mentioning that for the global indices the chemical power is directly related with the electronic density as well as the corresponding Hohenberg-Kohn functional [53]. As a complement of these global reactivity descriptors that arise from Conceptual DFT [5-7,51,52], Domingo and his collaborators [54][55][56][57][58] 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.
By making use of the mentioned CDFT software tool applied to the results of the calculation of the electronic properties of the Discodermins A-H peptides, the values of the defined global reactivity descriptors (including the Nucleophilicity N) could be obtained and they are displayed in Table 4: Table 4. Global reactivity descriptors for the Discodermin family of marine cyclopeptides: Electronegativity (χ), Hardness (η), Electrophilicity (ω) (all in eV), Softness S (in eV −1 ), Nucleophilicity N, Electrodonating Power (ω − ), Electroaccepting Power (ω + ) and Net Electrophilicity (∆ω ± ) (also in eV). The Global Hardness η can be seen as a measure of the resistance of the electronic density to be deformed and thus as an indication low reactivity of a given molecular system. From the results of Table 4, it can be concluded that Discodermin E will be the more reactive peptide in this family, while Discodermin G and Discodermin H will be the less reactive ones, being the chemical reactivity of the other peptides approximately the same. An analog behavior is observed for the Electrophilicity ω descriptor which encompasses the balance between the tendency of an electrophile to acquire an extra number of electrons and the resistance of a molecule to exchange electrons with the environment. As expected from the molecular structure of these species, their electron donating ability is more important than their electron accepting character, again the Discodermin E peptide possessing an electron donating power larger than the others. On the basis of the previous definition and the scale established by these authors [55], it can be concluded that the Discodermins A-H peptides can be regarded as strong nucleophiles because the values for the Nucleophilicity N are greater than 3 eV in all the cases. As a summary of these results the following scales can be given for each of the considered descriptors: for the Electronegativity (χ), Discodermin E > Discodermin A > Discodermin D > Discodermin C > Discodermin B > Discodermin F > Discodermin H > Discodermin G; for the Hardness (η), Discodermin H > Discodermin G > Discodermin F > Discodermin B > Discodermin D > Discodermin C > Discodermin A > Discodermin E; for the Electrophilicity (ω), Discodermin E > Discodermin A > Discodermin D > Discodermin C > Discodermin B > Discodermin H > Discodermin F > Discodermin G; for the softness S it will be Discodermin E > Discodermin A > Discodermin B ≈ Discodermin C ≈ Discodermin E > Discodermin F > Discodermin G ≈ Discodermin H, for the Nucleophilicity N, Discodermin G > Discodermin C ≈ Discodermin F > Discodermin A ≈ Discodermin B > Discodermin D > Discodermin H > Discodermin E, while for the Electrodonating Power (ω − ), Electroaccepting Power (ω + ) and Net Electrophilicity (∆ω ± ) will be the same scale, The presented global descriptors are a representation of the chemical reactivity of a molecule as a whole. However, local reactivity descriptors have been developed that can give an idea of the differences between the reactivity of each of the atoms that form the molecule. One of the most important groups of such reactivity descriptors are the Fukui functions [5][6][7] and the Dual Descriptor [59][60][61][62][63][64], which has been defined as: which are relationships between the electronic densities of the neutral, positive and negative species as well as between the Nucleophilic and Electrophilic Fukui functions. The Nucleophilic Fukui function, f + (r), reveals the sites on a molecular system that are susceptible to nucleophilic attacks and the Electrophilic Fukui function, f + (r), describes those sites that are more susceptible to electrophilic attacks. These local reactivity descriptors are very useful and have been used successfully for the identification of reactive sites. However, sometimes there is an overlap between the results of both descriptors and no conclusions can be accurately obtained. Instead, the Dual Descriptor ∆ f (r) or DD, can describe unambiguously nucleophilic and electrophilic sites within a molecule [64]. Thus, a graphical representation of the DD for the Discodermins A-H marine peptides is presented in Figure 5 showing clearly the areas within the molecules where DD > 0 and DD < 0 for a better understanding of the local chemical reactive of these molecules: Another local reactivity descriptors are the Parr functions [65,66] that can be considered to be an alternative to the Fukui functions for describing sites or areas within the molecules where nucleophilic or electrophilic attacks will be favored. The Parr functions can be expressed as [65,66]: where ρ rc s (r) and ρ ra s (r) are related to the atomic spin density of the radical cation or anion of the considered system, respectively [58].
To perform a comparison between the results that can be obtained from both formulations, the predictions for the specific reactions sites coming from the Electrophilic and Nucleophilic Fukui functions analysis for each peptide have been compiled in a series of tables that are presented as Supplementary Materials in the form of an Appendix to this work together with the values that are the outcome of the Parr functions analysis. The Radical Fukui function f 0 (r), which can be considered to be an average of f + (r) and f − (r), and denotes the favorable sites for a radical attack, has also been included for the sake of completeness. In a complementary way to the tables, a graphical representation of these descriptors have been included as a series of figures such as the comparison between the results coming for both kinds of studies can be done accurately.  Tables A1-A8 and Figures A1-A8 in Appendix A, it can be concluded that there is a nice agreement between the values coming from both formulations. In checking the numerical part, it can be seen that in all cases the reactive sites predicted by the Fukui functions or the Parr functions are the same. Indeed, the numerical values for the Parr functions are greater than for the Fukui functions with the meaning that the first ones are well defined. This observation is reinforced by checking at the graphical representations of the Fukui functions and the Parr functions, where although are extended over the same group of atoms or regions within the molecules, it can be seen that the Parr functions are more compact than the Fukui functions. Notwithstanding, this nice agreement could lead to the conclusion that if the Dual Descriptor DD can be built as the difference between both Fukui functions, a similar new reactivity descriptor could be defined in terms of the difference between both Parr functions giving the same information and being the starting point for future research on the field of chemical reactivity.

Conclusions
The Discodermins A-H family of cyclodepsipeptides of marine origin has been studied by resorting to some techniques of common use in the process of drug discovery and development showing that these kind molecules can be regarded as potential therapeutic drugs.
With this knowledge in mind, the chemical reactivity of the studied peptides has been exhaustively analyzed through the optimization of their structures using an MN12SX/Def2SVP/H2O model chemistry and the determination of their electronic properties by means of an improved model chemistry, namely MN12SX/Def2TZVP/H20, already used in previous works for the study of peptides, demonstrating their usefulness for this kind of calculations.
The chemical reactivity of the considered molecular systems was subject to an analysis based on a particular methodology developed by our research group named Conceptual DFT-based Computational Peptidology being the use of the MN12SX density functional validated once again by resorting to the KID procedure and the in-house software tool CDFT.
The analysis of the global and local reactivity descriptors arising from Conceptual DFT together with some proposals like the Nucleophilicity N and the Parr functions allowed for a complete understanding of the chemical reactivity of the studied peptides, by distinguishing the different chemical reactivities through the analysis of the global descriptors and the further identification of the reaction sites or regions within the molecules by resorting to Fukui and Parr functions as well as the Dual Descriptor.
Finally, a nice agreement was found between the outcome of the numerical and graphical analysis of the Fukui and Parr functions on each peptide leading to the conclusion that it could be possible to define a new chemical reactivity descriptor analog to the Dual Descriptor on the basis of the difference between the Parr functions which could act as an additional tool for the study of the chemical reactivity of molecular systems.

Abbreviations
The following abbreviations are used in this manuscript: