QuBiLS-MAS, open source multi-platform software for atom- and bond-based topological (2D) and chiral (2.5D) algebraic molecular descriptors computations

Background In previous reports, Marrero-Ponce et al. proposed algebraic formalisms for characterizing topological (2D) and chiral (2.5D) molecular features through atom- and bond-based ToMoCoMD-CARDD (acronym for Topological Molecular Computational Design-Computer Aided Rational Drug Design) molecular descriptors. These MDs codify molecular information based on the bilinear, quadratic and linear algebraic forms and the graph-theoretical electronic-density and edge-adjacency matrices in order to consider atom- and bond-based relations, respectively. These MDs have been successfully applied in the screening of chemical compounds of different therapeutic applications ranging from antimalarials, antibacterials, tyrosinase inhibitors and so on. To compute these MDs, a computational program with the same name was initially developed. However, this in house software barely offered the functionalities required in contemporary molecular modeling tasks, in addition to the inherent limitations that made its usability impractical. Therefore, the present manuscript introduces the QuBiLS-MAS (acronym for Quadratic, Bilinear and N-Linear mapS based on graph-theoretic electronic-density Matrices and Atomic weightingS) software designed to compute topological (0–2.5D) molecular descriptors based on bilinear, quadratic and linear algebraic forms for atom- and bond-based relations. Results The QuBiLS-MAS module was designed as standalone software, in which extensions and generalizations of the former ToMoCoMD-CARDD 2D-algebraic indices are implemented, considering the following aspects: (a) two new matrix normalization approaches based on double-stochastic and mutual probability formalisms; (b) topological constraints (cut-offs) to take into account particular inter-atomic relations; (c) six additional atomic properties to be used as weighting schemes in the calculation of the molecular vectors; (d) four new local-fragments to consider molecular regions of interest; (e) number of lone-pair electrons in chemical structure defined by diagonal coefficients in matrix representations; and (f) several aggregation operators (invariants) applied over atom/bond-level descriptors in order to compute global indices. This software permits the parallel computation of the indices, contains a batch processing module and data curation functionalities. This program was developed in Java v1.7 using the Chemistry Development Kit library (version 1.4.19). The QuBiLS-MAS software consists of two components: a desktop interface (GUI) and an API library allowing for the easy integration of the latter in chemoinformatics applications. The relevance of the novel extensions and generalizations implemented in this software is demonstrated through three studies. Firstly, a comparative Shannon’s entropy based variability study for the proposed QuBiLS-MAS and the DRAGON indices demonstrates superior performance for the former. A principal component analysis reveals that the QuBiLS-MAS approach captures chemical information orthogonal to that codified by the DRAGON descriptors. Lastly, a QSAR study for the binding affinity to the corticosteroid-binding globulin using Cramer’s steroid dataset is carried out. Conclusions From these analyses, it is revealed that the QuBiLS-MAS approach for atom-pair relations yields similar-to-superior performance with regard to other QSAR methodologies reported in the literature. Therefore, the QuBiLS-MAS approach constitutes a useful tool for the diversity analysis of chemical compound datasets and high-throughput screening of structure–activity data.Graphical abstract . Electronic supplementary material The online version of this article (doi:10.1186/s13321-017-0211-5) contains supplementary material, which is available to authorized users.


If I have seen further it is by standing on the shoulders of giants.
Isaac Newton, 1676.

Background
The codification of chemical information using mathematical-computational methods to accelerate small-molecule drug discovery constitutes one of the fundamental tasks of mathematical chemistry [1,2]. In recent years, the number and diversity of molecular features, also known as molecular descriptors (MDs), has significantly increased and corresponding educational and commercial computational implementations developed [3][4][5][6][7][8][9]. The absence of an ultimate universal chemical descriptor emphasizes the need of defining alternative methods to codify relevant and orthogonal chemical information.
In previous reports, Marrero-Ponce et al. proposed algebraic formalisms for characterizing topological (2D) and chiral (2.5D) molecular features through atom-and bond-based ToMoCoMD-CARDD (acronym for Topological Molecular Computational Design-Computer Aided Rational Drug Design) molecular descriptors [10][11][12][13]. These MDs codify molecular information based on the bilinear, quadratic and linear algebraic forms and the graph-theoretical electronic-density and edge-adjacency matrices in order to consider atom-and bond-based relations, respectively. The ToMoCOMD-CARDD MDs have been successfully applied in the screening of chemical compounds of different therapeutic applications ranging from antimalarials [14], trichomonacidals [15,16], antitrypanosomals [17], paramphistomicides [18], antibacterials [19], tyrosinase inhibitors [20,21] and others [22,23]. To compute these descriptors, a computational program with the same name was developed. However, this software barely offered the functionalities required in contemporary molecular modeling tasks, in addition to the inherent limitations that made its usability impractical, for instance: (a) it did not support standard input formats (i.e. MDL MOL/SDF files) and the only input method for the chemical structures entailed the sketching of molecular pseudographs using a built-in manual drawing mode; (b) parameter configurations could not be exported or saved for posterior experiments; (c) no option for batch processing of descriptors was offered; (d) lacked the distributed computing functionality which permits the correct utilization of current multi-core architectures; (e) could not be used as a standalone library thus preventing the its integration in other applications; and (f) presented ambiguities in the labeling of the descriptors' names in the output file.
In addition, in several mathematical procedures employed to compute MDs (e.g. GT-STAF [24,25], DIVATI [26] and QuBiLS-MIDAS [27][28][29][30]), the molecules are not analyzed as a whole, that is, these are partitioned in order to univocally characterize each atom independently. In this way, several mathematical operators (also known as aggregation operators) may be applied over the atom-level indices to compute different global/ local MDs. The use of several aggregation operators is based on the idea that the most suitable global definition of a system may not necessarily be additive. In fact, it is reported in the literature that operators other than the sum could yield better correlations with determined chemical properties [24][25][26][27][28]. In this sense, in the present report strategies are defined to generalize the procedure of obtaining global or local QuBiLS-MAS (acronym for Quadratic, Bilinear and N-Linear mapS based on graphtheoretic electronic-density Matrices and Atomic weight-ingS) indices using the so-called aggregation operators. Moreover, several new atom-based properties, chemical local-fragments (e.g. terminal methyl groups), distancebased cut-offs (for the analysis of the most important non-covalent or covalent interactions) and probabilistic transformations of the matrix representations are introduced. Furthermore, initiatives to deal with the computational and practical limitations inherent to the original ToMoCoMD-CARDD program were carried out, with the ultimate goal of improving its applicability in presentday cheminformatics tasks. extensions and generalizations implemented in this software is demonstrated through three studies. Firstly, a comparative Shannon's entropy based variability study for the proposed QuBiLS-MAS and the DRAGON indices demonstrates superior performance for the former. A principal component analysis reveals that the QuBiLS-MAS approach captures chemical information orthogonal to that codified by the DRAGON descriptors. Lastly, a QSAR study for the binding affinity to the corticosteroid-binding globulin using Cramer's steroid dataset is carried out.
� ∃ e ij : e ij ∈ E 0 otherwise is built with values computed from the properties corresponding to the atoms that each bond connects [10,13,20,40]: where, w ij constitutes the weighting scheme computed for the edge e ij , w i and w j are the atomic weights (e.g. electronegativity) for atoms i and j forming the considered bond (e ij ), δ i and δ j are the corresponding vertex degrees which also account for bond multiplicity. Moreover, in order to codify information on the 3D structure of the molecule, a trigonometric 3D-chirality correction factor is applied to the molecular vectors aforementioned, which has been comprehensively explained in several reports [40][41][42].
From the previous molecular vectors and matrix formalisms, the algebraic calculation of the NS and SS total (whole-molecule) bilinear indices may be condensed in the following equations, for atom-(see Eq. 3) and bondbased (see Eq. 4) indices, respectively: where, n (or m) is the number of atoms (or bonds) in the molecule, k = 1, 2, …15 is the matrix power, m k ij (or e k ij ) represents the elements of the M k ns,ss (or E k ns,ss ) nonstochastic (ns) and simple stochastic (ss) matrices, and x i and y j are the elements of the x and ȳ atom-based (or bond-based) property vectors. On one hand, when the vectors x and ȳ encode the same atomic property (i.e. x =ȳ), the Eqs. 3 and 4 define the NS and SS total atombased and bond-based quadratic indices, respectively. On the other hand, if x is a vector with all entries equal to 1 and ȳ an atom/bond-based property vector, then the Eqs. 3 and 4 define the NS and SS total atom-based and bond-based linear indices, respectively.
In addition, local-fragment (group or atom-type) quadratic, bilinear and linear atom/bond-based indices can be defined to characterize a predetermined molecular fragment (F) instead of the whole molecule (total indices). These are computed using the kth local-fragment matrix F M k ( F E k ), which is computed from the corresponding kth total matrix M k (E k ) considering only those vertices (or edges) belonging to the selected molecular fragment. These fragments F may be heteroatoms (X), halogens (G) It is important to remark that for each partitioning of a molecule into Z molecular exclusive fragments, there will be Z local-fragment matrices. In this case, if a molecule is partitioned into Z molecular fragments, then the original kth power of matrix M k ns,ss (or E k ns,ss ) is exactly the sum of the kth power of the local-fragment matrices. Consequently, the total algebraic form-based indices are the sum of the exclusive contributions of the respective local-fragment algebraic form-based indices, as long as there is not overlap among the fragments. Therefore, taking into consideration the previous elements, the next sections address in detail the improvements related with the mathematical definition corresponding to the 2D algebraic indices introduced by Marrero-Ponce et al. [10,31,32,43,44].

The QuBiLS-MAS MDs: new definitions, generalization and extension of algebraic indices
As previously explained, up to date, the 2D atom/bondbased algebraic indices have been computed as wholemolecule (total) indices or from specific chemical groups (local indices), where the simplest fragment could be the atom itself, known as a LOcal Vertex Invariant (LOVI) and in case of a bond as LOcal Edge Invariant (LOEI). In this manuscript the LOVEIs term is adopted to refer both LOVIs and LOEIs of a molecule, and is denoted as  L. Therefore, if a molecule is comprised of n atoms or m bonds then the kth total bilinear, quadratic and linear indices for each atom "a" (known as total atom-level index) or each bond "e" (known as total bond-level index) may be computed as two-linear algebraic forms (maps) in R n , in a canonical basis set, and whose values are components (entries) of the vector L denoted as L a and L e for atom-and bond-level indices, respectively. In this way, the kth total atom-level and bond-level bilinear indices are mathematically defined as follows, respectively: where x 1 , …, x n(m) and y 1 , …, y n(m) are the coordinates or components of the molecular vectors x and ȳ [45]. To compute these molecular vectors the following atomic properties have been selected: (1) atomic mass, (2) the Van der Waals volume, (3) the atomic polarizability, (4) atomic electronegativity according to Pauling scale, (5) atomic Ghose-Crippen LogP, (6) atomic Gasteiger-Marsili charge, (7) atomic polar surface area, (8) atomic refractivity, (9) atomic hardness and (10) atomic softness. These properties are calculated using the CDK library [9]. Note that when x =ȳ atom-and bond-level quadratic indices are obtained [i.e. q L a = q a,k (x,x) and q L e = q e,k (x,x)], while if all coefficients of x are equal to 1 then linear indices for atoms (or bonds) may be obtained [i.e. f L a = f a,k ū,ȳ and f L e = f e,k ū,ȳ ]. The coefficients m ij a,k (see Eq. 5) are the elements corresponding to the kth NS (or SS) total atom-level pseudograph-theoretic electronic-density matrix [NS(SS)-GEDM] M a,k for atom "a", while the entries e e,k ij (see Eq. 6) belonging to kth NS (or SS) total bond-level edge-adjacency matrix [NS(SS)-EAM] E e,k for bond "e". These atom/bondlevel coefficients are obtained from the entries m ij k of the M k total matrix and e k ij of the E k total matrix, respectively, using the described procedure to compute local-fragment matrices but considering the fragment F as an atom "a" or bond "e" of the molecule. Moreover, the diagonal coefficients m ii 1 could have two distinct values in order to achieve greater discrimination of molecular structures: (1) aromatic ring sensibility for setting up aromatic atoms hooked on full aromatic rings instead of mapping individual atom loops as shown in the molecular pseudograph of the Table 1, and/or (2) the number of lone-pairs for each atom. The e 1 ii entries are always zero.
It is important to highlight that as an extension of the former ToMoCoMD 2D-MDs several local-fragments have been aggregated: H-bond acceptors (A), carbon atoms in aliphatic chains (C), H-bond donors (D), halogens (G), terminal methyl groups (M), carbon atoms in an aromatic portion (P) and heteroatoms (X). Thus, from these local-fragments the kth NS (or SS) local-fragment atom-level pseudograph-theoretic electronic-density matrices F M a,k for atom "a" and the kth NS (or SS) local-fragment bond-level edge-adjacency matrices F E e,k for bond "e", may be computed. Consequently, local-fragment atom-and bond-level bilinear, quadratic and linear indices are determined from the Eqs. 5 and 6 using F M a,k and F E a,k as matrix forms, respectively. Note that the coefficients F m a,k ij ∈ F M a,k and F e e,k ij ∈ F E e,k are calculated from the elements F m k ij ∈ F M k and F e k ij ∈ F E k , respectively. In addition, two normalization procedures are introduced as novel extensions. The atom-based simple stochastic scheme defined in the original ToMoCoMD 2D-MDs [18,39,43] describes changes in the electron distribution over time throughout the molecular backbone. This SS matrix is not symmetrical and the probability for atom i to interact with atom j is different from the probability for the atom j to interact with the atom i. Therefore, with the aim of balancing the probabilities in both senses a double-stochastic (DS) matrix is employed, that is, a matrix with real non-negatives entries whose column and row sums are equal to one. In this way, the kth total (or local-fragment) DS graph-theoretical electronicdensity (DS-GEDM, (F ) M k ds ) and edge-adjacency (DS-EAM, (F ) E k ds ) matrix approaches can be calculated from the corresponding M k ns and E k ns matrices, respectively, using the Sinkhorn-Knopp algorithm [46]. Additionally, the kth total (or local-fragment) mutual probability (MP) graph-theoretical electronic-density matrix (MP-GEDM, (F ) M k mp ) and edge-adjacency matrix (MP-EAM, (F ) E k mp ) are introduced. The mutual probability matrices are obtained dividing each entry between the total sum of their elements, in this way, symmetrical matrices where the total sum is equal to 1 are obtained. The Scheme 1 shows the steps followed in the computation of the NS-, SS-, DS-and MP-GEDMs, while Tables 2 and 3 illustrate the calculation of these matrices with and without taking in consideration the lone-pair electrons.
Lastly, in order to obtain the global kth total (or localfragment) bilinear, quadratic and linear indices from the corresponding atom-level (L a ) or bond-level (L e ) definitions, the summation operator is used. The global indices obtained using this operator over components of vector L coincide with those indices calculated through the original procedure vector-matrix-vector detailed in Eqs. 3 and 4. Note that the summation operator is equivalent to the Manhattan norm applied to elements of the vector L relative to the origin, which is in turn a specific case of Minkowski norm when p = 1. Motivated by this understanding, a generalization in which different p values are used, i.e. p = 2 and 3, where the former (p = 2) is the Euclidean norm (see Additional file 1: Figure SI1 for geometrical interpretation) was introduced. Additionally, other operators (see Additional file 1: Table SI2) applicable to the vector of LOVEIs were applied with the aim of generalizing the use of the linear combination to obtain global indices. It has been demonstrated in several reports [24][25][26][27][28] that better correlations for bioactivities may be attained when operators other than the sum are employed.

Neighborhood topological constraints in the graph-theoretical electronic-density and edge-adjacency matrix
The (F ) M k and (F ) E k matrices contain information on the connectivity for all atoms and bonds that constitute a molecule, respectively. However, some biological properties do not depend on the chemical structure as a whole but rather on interactions at particular topological distances, for example, short-, middle-and large-range contacts. Thus, with the aim of considering interactions that satisfy specific topological criteria, three graph-theoretical constraints (cut-offs) are introduced: (1) keeping only the diagonal elements of the matrix, denoted as "Self-Returning Walks" (SRW), (2) keeping only the off-diagonal elements of the matrix, denoted as "Non-Self-Returning Walks" (NSRW), and (3) keeping only the elements within a given interval, based on the topological distance for a path cut-off, denoted as Lag p.
The application of these cut-offs over the matrices where, p ij is a user-defined topological distance threshold, and min and max are the minimum and maximum cut-off values (rank). Table 4 shows an illustrative example where three topological constraints are calculated for an atom-level matrix. A custom cut-off allows to distinguish the interaction types, for example, when a topological graph-theoretical cut-off is applied, then atomic Table 2 The

molecular structure considering lone-pair electrons (n) for the first and second powers of the molecular pseudograph's atom adjacency mutual probability (mp), non-(ns), double (ds)-and stochastic (ss) matrices for Isonicotinic Acid
Isonicotinic Acid displaying Lone-pairs Electrons indices could be calculated for atoms separated by 1 step (covalent interactions) or for those atoms separated by more than 1 step (p ≥ 2). The present approach could be viewed as a threshold that generalizes the use of lag p in 2D-Moreau-Broto autocorrelations [1]. Likewise, these matrices based on cut-offs may be employed to determine the corresponding atom-level and bond-level representations to be used in the calculation of QuBiLS-MAS 2D-MDs. In Scheme 2, a complete workflow to compute the QuBiLS-MAS indices is represented.

The QuBiLS-MAS module
The QuBiLS-MAS module was designed as standalone software, with the extensions and generalizations discussed in "The QuBiLS-MAS MDs: new definitions, generalization and extension of algebraic indices" section. This software was developed in Java v1.8 and the Chemistry Development Kit (CDK) library (version 1.4.19) [9] was used in the manipulation of the chemical structures, as well as in determining the atom-and fragment-based chemical properties involved in the calculation process. The QuBiLS-MAS software is comprised of a front-end and back-end. The front-end is composed of a desktop and command-line user interface, while the back-end is developed as an Abstract Programming Interface (API) to enable its use as an independent Java library in the development of other cheminformatics applications or in the implementation of other user-friendly interfaces either graphical or command-line based. With these two components, independence between the software presentation layer and the processing logic implemented in the back-end is achieved and thus, any modification in the latter does not provoke changes in the front-end (GUI), and vice versa.

Back-end: the QuBiLS-MAS molecular descriptors library-computational complexity of algorithms
All the requests performed by the users through the GUI are processed by the QuBiLS-MAS library. This component is structured in packages according to the goals of Table 4 First, second and third order NS-matrices for Isonicotinic Acid, obtained by applying three types of topological constraints (cut-off): Self-Returning Walks (SRW), Non-Self-Returning Walks (NSRW) and a topological path cut-off distance from 2 to 5 (LAG the functionalities (see Additional file 1: Figure SI3 for UML diagram). The main package is tomocomd.cardd. qubils which contains the packages descriptors, matrices, metrics and workers that encapsulate the main concepts utilized in the definition of the QuBiLS-MAS MDs.
The descriptors package includes the classes related to the calculation of the total and local-fragment bilinear, quadratic and linear algebraic maps. The matrices package contains the objects responsible for building the pseudograph-theoretic electronic-density matrix and the edge-adjacency matrix corresponding to atom-and bond-based representations, respectively. Additionally, the simple-stochastic, double-stochastic and mutual probability normalization strategies, as well as the topological constraints (cut-offs) are defined in this package. The tools package includes classes for the identification of the local-fragments, as well as the considered aggregation operators. Lastly, the workers package comprises the classes for the configuration and control of the algebraic MDs calculation process. The algorithms responsible for performing the multiplication based on bilinear, quadratic and linear algebraic forms constitute the principal procedures to compute the QuBiLS-MAS indices. This procedure consists of a loop that iterates for each atom of the molecule to determine the corresponding atom-or bond-level matrix. Next the atom/bond-level matrices are multiplied by the corresponding property vectors in order to obtain the atom/ bond-level indices. The corresponding sequential implementations have a computational complexity of O(n 3 ). Nonetheless, when the atom/bond-level matrices are computed according to the mentioned procedure, it is noted that the only entries with values different from zero correspond to the atom with respect to which the atom/bond-level matrix is built. Therefore, instead of iterating for each atom in order to build the atom/bondlevel matrix used posteriorly to determine the corresponding index, it is more suitable to compute the atom/ bond-level indices at the same time as the original matrix is analyzed. Taking this into account, the algorithms have been optimized to an inferior polynomial order, achieving a complexity of O(n 2 ) in the computation of the atom/ bond-based contributions for the QuBiLS-MAS indices.

Graphic user interface of the QuBiLS-MAS software
To facilitate the calculation of the QuBiLS-MAS MDs, a friendly Desktop GUI was developed in order to provide a simple and intuitive way to configure the Schema 2 Workflow followed in the computation of the ToMoCoMD-CARDD QuBiLS-MAS MDs different parameters used, such as: algebraic forms, matrix approaches, atomic properties, topological cutoffs and so on. Figure 1 shows the main GUI and the dialog windows designed to configure some of these parameters. These configuration sections allow the users to personalize the bilinear, quadratic and linear indices according to their necessities and thus predefined MDs are not calculated.
In the "Algebraic Form" panel, the specific algebraic maps to be used in the computation of the MDs are chosen according to the selected option in the "Constraints" panel, which could be atom-based or bond-based. Also, chirality detection may be configured in the "Constraints" panel. The matrix normalization formalisms (MP, NS, SS, and DS) used in the algebraic forms are configured in the "Matrix Form" panel, as well as the maximum order (k value) to which the coefficients of the matrices are raised. In the "Cut-Off " panel the option to "keep all" (KA) atomic interactions is selected by default, but other options [i.e. "Self-Returning Walks" (SRW), "Non-Self-Returning Walks" (NSRW) and/or the value-rank(s) of threshold p] may be considered to take into account only the non-covalent interactions according to the established criterion. The "Local-Fragments" panel contains the options to configure the seven chemical groups (or atom-types) that may be employed to compute either the total or local-fragment indices. Likewise, in the "Properties" panel the atomic properties used to setup different weighting schemes are chosen. Finally, the mathematical operators used to compute the global total or local indices from the atomic contributions are selected in the "Invariants" panel.
It is important to highlight that the selected options to compute the descriptors can be exported into an XML configuration file, called the project file, which can be used to calculate the same QuBiLS-MAS indices for other datasets when the software is run again. Another important feature is that the software can be executed on computer clusters using a command-line interface, which uses the project files to obtain the configuration of the indices to be computed. Also, the QuBiLS-MAS software has incorporated the "On/Off H-Atoms" option to consider (or not) the H-atoms during the calculation, the "On/Off Lone-Pair Electron" option to consider (or not) Fig. 1 Main graphic user interface for QuBiLS-MAS software (a) and dialog windows to configure the following parameters: invariants or aggregation operators (b), atom properties (c) and local-fragment chemical groups (d) the number of lone-pairs for heteroatoms and the "Show Debug Report" option to track the algebraic processes that take place during the calculation (see Additional file 1: SI4).
The supported input file format for the chemical structures to be analyzed is the MDL MOL/SDF format and these are sequentially read in order to employ suitable memory allocation according to the size of the molecule. Moreover, the path of the output file may be specified where the values of the computed MDs are saved. To this end, the QuBiLS-MAS software supports the following output file formats: CSV, ARFF, and TXT (space-and tabseparated ASCII format) which are easily interpretable in popular statistical and/or machine learning software.
The calculation procedure is monitored in real time through the main interface and controlled with the interactive mode of the GUI. Indeed, more than one project file can be calculated over different datasets. This is a feature implemented in the QuBiLS-MAS software encapsulated into a batch processing module, which is useful for carrying out high-throughput and routine MD calculations. This module is designed to manage the configuration of up to eight independent tasks (see Additional file 1: SI5), where each task consists of one or several datasets for which one or several projects files previously saved with the QuBiLS-MAS GUI may be computed. Finally, a module for chemical structure curation tasks was incorporated, taking into account Tropsha's guidelines [47]. Table 5 shows a comparison between the old [48] ToMoCOMD software and the present one (QuBiLS-MAS module), highlighting the numerous functionalities incorporated. Table 6 compares the characteristics for common molecular descriptor calculating software and including the QuBiLS-MAS program, specifying the respective strengths and weaknesses.

Information content analysis based on Shannon's entropy
Shannon's entropy (SE) quantifies the information content codified by molecular indices, according to the principle that variables that effectively discriminate all molecules in a dataset possess high entropy values, while redundant variables have low entropy values. To perform this study, the Spectrum dataset (http://www.msdiscovery.com/spectrum.html) comprised by 1963 structures was used. The highest SE for this dataset is equal to 10.93 bits (log 2 N, where N is the number of compounds). In the following subsections the novel QuBiLS-MAS 2D-MDs are analyzed taking into account the proposed internal generalizations, as well as with respect to well-known MDs computed by other software. For this study, the IMMAN software was used [49].

Comparative variability analysis according to the matrix formalisms
The four matrix schemes defined in the present report are analyzed. To this end, 880 MDs are calculated for each matrix. Figure 2 shows similar entropy distributions for the non-, double-and simple-stochastic matrix approaches, while the best behavior is obtained with the mutual probability approach. The superior performance of the mutual probability formalism with respect to the other three matrix transformations justifies the theoretical contribution of this scheme in the computation of the QuBiLS-MAS 2D-MDs.

Analysis of variability according to the aggregation operators
The aim of this section is to evaluate the variability of the QuBiLS-MAS 2D-indices according to the mathematical operators used over the vector of LOVEIs. In this study, the aggregation operators classified as norms, means and statistical invariants are compared. To this end, 110 atom-based linear indices for each operator were calculated and the results are shown Fig. 3. As it can be noted, the best results are achieved by the Potential Mean, Quadratic Mean and Standard Deviation operators with 71, 67, 66 and 65% of the total variables having entropy values greater than 9.0 bits (82% of the maximum entropy), respectively. Moreover, the indices based on the Manhattan (sum of LOVEIs) and Minimum operators present the worst performance, while the remaining distributions have similar behavior. This result suggests that the generalization of the linear combination of LOVIEs to consider other aggregation operators yields variables with greater information content, and thus, it should contribute to a greater modeling capacity for the QuBiLS-MAS MDs.

Variability analysis of QuBiLS-MAS 2D-indices versus DRAGON descriptor families
The purpose of this analysis is to compare the entropy of the QuBiLS-MAS 2D-MDs with the DRAGON descriptor families. To perform this study some DRAGON descriptor-blocks were clustered into bigger families: (1) 0D_others for molecular properties, constitutional and charge descriptors, (2) 1D-fragment for functional group counts and atom-centered fragments, (3) 2D-conn_autocorr_inf for 2D autocorrelations, connectivity and information indices, (4) 2D-edge_walk for edge adjacency indices, walk and path counts, (5) 2D-eigenvalues for Burden eigenvalues, topological charge and eigenvalues-based indices, and (6) 3D-Ran-dic_geometrical for Randic molecular profiles and geometrical descriptors. The remaining DRAGON families were kept with the same denominations. The maximum number of descriptors considered for each family is 91, which corresponds to the 0D_others family that has the least number of MDs.
As it can be observed in Fig. 4, the QuBiLS-MAS 2D-MDs show the best overall performance with all the considered indices presenting entropy values

Features
Computer   . Chirality trigonometric 3D-chirality correction factor (c). Topological constraints (cut-offs) (1) keeping only the diagonal elements of the matrix, denoted as "Self-Returning Walks" (SRW), (2) keeping only the offdiagonal elements of the matrix, denoted as "Non-Self-Returning Walks" (NSRW), and (3) keeping only the elements within a given interval, based on the topological distance for a path cut-off, denoted as Lag p

Variability comparison for QuBiLS-MAS 2D-indices with respect to other descriptor computing software
The variability distribution of the QuBiLS-MAS MDs was computed and compared to MDs calculated with other programs used in cheminformatics tasks, such as: DRAGON [3], MOLD2 [4], PADEL [7], _ENREF_70 CDK Descriptor Calculator [9], MODESLAB [50], BLUECAL [51] and POWER MV [52]. To this end, the DRAGON's example data comprising 42 structurally diverse chemicals was used. The cut-off number of variables for this study was 170 MDs, determined by the BLUECAL software as it possesses the least number of indices. As it can be observed in Fig. 5, the QuBiLS-MAS topological indices achieve superior performance than other software considered, with the former presenting all its values above 4.62 bits [86% of the maximum entropy (log 2 41 = 5.35)], while the indices of the remaining approaches practically have all their indices inferior to this threshold. The high entropy distribution obtained for the QuBiLS-MAS topological indices demonstrates the relevance of these MDs, in the sense that they are sensitive to progressive structural modifications and should therefore be valuable in different cheminformatics tasks.

Linear independence of the QuBiLS-MAS algebraic descriptors
In this section, the possible orthogonality of the QuBiLS-MAS 2D-Indices with respect to the DRAGON 0D-2D MDs is examined, using the Principal Component Analysis (PCA) [53,54]. The PCA is a mathematical technique that converts several correlated variables into a reduced number of non-correlated variables, called principal components. The extracted components have the following features: (1) the first component will explain the highest possible variance of all determined components, (2) the successive components will explain the variance that the previous components did not explain, and (3)  In this analysis, 12 principal components were selected, which explain approximately 74.60% of the cumulative variance (see Additional file 1: SI6 and Additional file 1: SI7). As it can be observed, Factors 1 (27.83%), 2 (13.06%), 8 (2.47%) and 9 (1.99%) exhibit strong loadings for some QuBiLS-MAS indices and some 0D-2D descriptors of the DRAGON software. On the other hand, exclusive loadings are obtained for the QuBiLS-MAS descriptors in the Factors 3 (8.6%), 4 (6.26%), 5 (3.86%), 6 (3.51%), 7 (2.71%), 11 (1.42%) and 12 (1.20%), explaining 27% of the total variance. Factor 10 (1.62%) is important for some 0-2D DRAGON MDs as these are exclusively loaded in this factor, and these indices include: TI2 (B02 2D Topological), PW2 (B02 2D Topological), RBF (0D-others) and EEig01r (2D-edge_walk) [for details on these descriptors, see Additional file 1: SI8]. On the whole, much of the information codified by the 0D-2D DRAGON MDs is equally captured by the QuBiLS-MAS indices, considering that negligible variance (1.62%) is explained by the factor exclusive for the former (F10). Moreover, the numerous factors (i.e. F3, F4, F5, F6, F7, F11 and F12) exclusive for the QuBiLS-MAS MDs suggest that orthogonal information is codified and thus demonstrating the theoretical contribution of the generalization schemes adopted in this framework.

QSAR modeling of the binding affinity to corticosteroid binding globulin (CBG) of Cramer's steroid dataset
In what follows, the predictive ability of the QuBiLS-MAS approach is assessed. To accomplish this objective, QSAR models for predicting the "binding affinity to the corticosteroid-binding globulin (CBG) of the popular Cramer's steroid database" (see Additional file 1: SI9 for names and CGB values of compounds) were built. This dataset has been used as a "benchmark" to evaluate the quality of novel procedures. A total of 1455 variables were computed for each algebraic form (quadratic, bilinear and linear maps). The prediction models were built using Multiple Linear Regression (MLR) as the fitting method, coupled with the Genetic Algorithm (GA) as variable subset selection strategy and the statistical parameter Q 2 loo ("leave-one-out" cross validation) as the fitness function. Throughout the study, regression models of 2-6 variables were developed and the best model in each case retained for posterior validation. The GA was setup with the following configurations: population size-100, crossover/mutation rate-0.7, selection operator was fixed at 60 and the number of iterations-500,000. In addition, the tabu list option was configured to remove those MDs with correlation equal or greater than 0.95. The MLR-GA based model building was performed using the MobyDigs [55] computer program. The best models built were also assessed with the bootstrapping [56] (Q 2 boot ) and Y-scrambling [57] (a(Q 2 )) validation methods in order to assess the predictive power and the possible chance correlation with respect to the activity modeled.

Examination of matrix formalisms
In order to assess the performance of the NS, SS, DS and MP matrix-based approaches in QSAR modeling, 46 variables for each formalism were calculated. Figure 6a shows the statistical parameters achieved in this experiment, where the SS approach (Q 2 loo = 81.85%, Q 2 boot = 77.89%) presents the best behavior, followed by MP (Q 2 loo = 79.05%, Q 2 boot = 74.85%). The indices based on NS (Q 2 loo = 73.48%, Q 2 boot = 68.09%) and DS (Q 2 loo = 72.01%, Q 2 boot = 65.4%) matrices present a much lower performance. This result is in agreement with the variability analysis, where the highest entropy indices involved the SS and MP matrix formalisms.

Analysis of the aggregation operators
The following study evaluates the predictive power of the aggregation operators proposed as a generalization scheme for the linear combination of LOVEIs as method for obtaining global (or local) indices. As it can be observed in Fig. 6b, all Q 2 loo values are superior to 50%, with the best performances corresponding to the statistical operators, followed by the mean operators and lastly by the norms. Regarding the evaluation of the operators classified as "classical algorithms" (Fig. 6c) it is observed that Kier-Hall (KH), Total Sum (TS), Gravitational (GV) and Autocorrelation (AC) algorithms yield comparable to superior performance with respect to the remaining operators. It may therefore be concluded that the incorporation of the aforementioned generalization scheme improves the performance of the QuBiLS-MAS indices in modeling tasks and thus demonstrating its practical contribution.

The QuBiLS-MAS MDs versus literature reports
To evaluate the earnest contribution of the QuBiLs-MAS approach, it is necessary to assess its performance in correlation studies with determined molecular properties and compare the results with the existing methods. Different QSAR models for predicting the binding affinity to CBG of the 31 structures of Cramer's steroid database (1-31 or also 1-30 with compound 31 as outlier) have been reported in the literature, which will be compared here with the models obtained using the QuBiLs-MAS 2D-MDs. In this experiment, the best 3-5 variable models were selected according to the quality of the statistical parameters Q 2 loo and Q 2 boot . Table 7 shows the best regression models and their corresponding statistical parameters, based on the QuBiLs-MAS 2D-indices. Comparisons with other QSAR methodologies reported in the literature are presented in Table 8 according to the Q 2 loo statistic.
In general, when the 31 steroids are taken into account as training set, the models based on QuBiLS-MAS indices yield comparable-to-superior performance relative to other methods reported in the literature according to the Q 2 loo statistic. Up to now, the best model reported  Combined electrostatic and shape similarity matrix 5 PLS 0.70 [59] SAMFA-RF -RF 0.69 [81] SAMFA-PLS 4-5 PLS 0.69 [81] 4D-QSAR 2 PLS 0.69 [69] CoMMA (ab initio) 6 PLS 0.689 [82] has been the one based on the "Combined Electrostatic and Shape Similarity Matrix" (Q 2 loo = 0.941, var = 6), which is an alignment-and grid-based method known to be computationally expensive. Additionally, this model employs the Genetic Neural Network (GNN) as the fitting method, which generally yields more robust and better optimized models compared to other linear methods. Even then, comparable performance is obtained with QuBiLs-MAS models [(Q 2 loo = 0.937 (compound 31 excluded), var = 6), (Q 2 loo = 0.914 (compound 31 included), var = 6)] based on the MLR-GA, which is a much simpler technique than GNN. Therefore, based on the results obtained in this study, it can be claimed that the QuBiLs-MAS MDs proposed offer a considerable advantage over well-known traditional methodologies.

Conclusions
The QuBiLs-MAS approach for atom-pair relations, in its diverse generalizations and extensions, seems to renew the prospect of achieving 2D-QSAR models with good predictive power. Inspired by the "No Free Lunch" theorem [58], which postulates that there is no unique best alternative for tackling optimization problems, the different extensions constitute an innovative undertaking to suitably characterize the different phenomena that affect the molecular configuration and intermolecular interactions, and thus affecting their biological activity. Variability and Principal Component analyses of the QuBiLs-MAS indices demonstrated that the proposed generalizations yield indices with superior variability compared to other indices defined in the literature and capture chemical information not codified by the DRAGON MD families. Also, it was demonstrated that suitable gains are obtained in the predictive ability of the QSAR models with the QuBiLs-MAS approach. Therefore, the QuBiLs-MAS 2D-indices constitute a relevant tool for the diversity analysis of compound datasets and high-throughput screening of structure-activity data.

Futures outlooks
Future tasks include the development of a version of the QuBiLs-MAS module to compute molecular indices on a distributed computing system for high-throughput calculation, as well as, a version to use the Graphical Processing Units (GPU) present in several personal computers nowadays. Moreover, various (dis-)similarity multi-metrics to consider relations for more than two atoms (multilinear forms) are to be introduced, in addition to a new set of multi-metrics based cut-offs.
Authors' contributions YMP proposed the theory of the QuBiLS-MAS indices, supervised the chemical applications, the design of the GUI and prepared the manuscript. JRVM worked in the definition of the QuBiLS-MAS indices, in the computational implementation of API and GUI interfaces, performed the QSAR and other statistical analysis and prepared the manuscript. YSVA worked in the computational implementation of QuBiLs-MAS software. KMM, SJB, HLT and FPG worked in the QSAR modeling and performed the statistical analysis. CRGJ and CAM lead the informatics (program design) research related with this manuscript. All authors read and approved the final manuscript.

Additional file
Additional file 1. The mathematical definitions of the norms, means and statistical invariants as generalizations of the linear combination of LOVIs as global (and/or local) MDs aggregation operator, as well as classical algorithms which generalize the first three groups are presented as Figure SI1-Table S12. The UML diagram (Figure SI3), a debug report file content ( Figure SI4), a batch process manager dialog window ( Figure  SI5) are also listed. Some results of the factor analysis by the principal component method are shown as Table SI6-Table SI8, and finally, the names of structures for Cramer's steroid database and their corresponding values for the binding affinity to the corticosteroid-binding globulin (CBG) is in Table SI9. Wagener's (AMSP Method) -k-NN and FNN 0.630 [84] SAMFA-SVM -SVM 0.60 [81] ALPHA 2 PLS 0.57 [67] Italic values indicate the results of QuBiLS-MAS approach a When it is applicable, specifies the number of components (PCs) b 1.0 A models c Compound 31 excluded, taken as outlier, is not taken into account in the training set † Logarithm of the binding affinity to the corticosteroid-binding globulin (CBG)