Skip to main content

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

If I have seen further it is by standing on the shoulders of giants.

Isaac Newton, 1676.



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.


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.


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.



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 [39]. 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 [1013]. 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 [2730]), 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 [2428]. 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 graph-theoretic electronic-density Matrices and Atomic weightingS) indices using the so-called aggregation operators. Moreover, several new atom-based properties, chemical local-fragments (e.g. terminal methyl groups), distance-based 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 present-day cheminformatics tasks.

Theoretical scaffold: past and present

Brief history of algebraic maps-based indices

The algebraic forms-based topological MDs (0–2.5D) are divided into three main families: quadratic, bilinear and linear indices [12, 31, 32]. They are distinguished in atom-based [33] and bond-based indices [10] depending on whether they are derived from the atom-based or bond-based matrix, respectively. The main diagonal elements for the atom-based matrix [denominated as non-stochastic (NS) when it doesn’t involve any normalization procedure] describe the presence of loops on graph vertices, which are used to characterize atoms in conjugated systems having more than one canonical structure [31, 34]. Thus, the elements for the kth non-stochastic pseudograph-theoretic electronic-density matrix \(({\boldsymbol{\mathcal{M}}}_{ns}^{k} )\) are labeled as \({}^{k}m_{ij}\) and defined as follows:

$${}^{k}m_{ij} = \left\{ {\begin{array}{*{20}ll} {P_{ij} } \hfill & if\,i\,\ne\,j\,{\bigwedge }\,\exists\,e_{ij}{:}\;e_{ij} \in \varvec{E} \hfill \\ {L_{ij} } \hfill & if\,i = j \,{\bigwedge }\,\exists\,e_{ij}{:}\;e_{ij} \in \varvec{E} \hfill \\ 0 \hfill & {otherwise} \hfill \\ \end{array} } \right.$$

where, \(i\) and \(j\) represent two vertices (atoms) of the molecular pseudograph G, k is the matrix power, \(\varvec{E}\) is the set of edges of G, \(P_{ij}\) is the number of edges \((e_{ij} )\) between the atoms \(i\) and \(j\) (e.g. \(P_{ij} = 3\) for a triple covalent bond between i and j), and \(L_{ij}\) is the number of loops in \(v_{i}\) [12, 13, 31, 33, 35, 36]. Likewise, the coefficients corresponding to the bond-based matrix, \({\boldsymbol{\mathcal{E}}}_{ns}^{\varvec{k}} ,\) may be defined. In this way, the entries \(e_{vw}\) belonging to \({\boldsymbol{\mathcal{E}}}_{ns}^{\varvec{k}}\) are equal to 1 if the edge \(v\) shares a common vertex with the edge \(w\) [37, 38]. Moreover, the NS matrix may be normalized by means of the simple stochastic (SS) procedure [10], yielding matrices whose row or column coefficients are non-negative real numbers which sum up to 1. This mathematical procedure has been explained in detail elsewhere [13, 18, 39]. Let us take a simple example of the isonicotinic acid structure, and consider its corresponding labeled molecular pseudograph and atom-based matrix [31]. Table 1 shows the non-stochastic (NS) matrix for the isonicotinic acid structure for k = 0, 1, 2.

Table 1 The molecular structure and the atom adjacency stochastic (ss) and non-stochastic (ns) matrices for the k = 0, 1, 2 corresponding to the Isonicotinic Acid

To compute the algebraic form-based indices, the molecular vector concept is employed, which uses atom-based properties as weighting schemes. Thus, atomic properties (e.g. mass, polarizability, electronegativity according to Pauling’s scale and Van der Waals volume) may be considered [11, 12]. In this way, the molecular structures may be represented as vectors. For instance, the Isonicotinic Acid molecule may be represented by the molecular vector \(\bar{x} = \left[ {x_{N1} ,x_{C2} ,x_{C3} ,x_{C4} ,x_{C5} ,x_{C6} ,x_{C7} ,x_{O8} ,x_{O9} } \right]\), where \(\bar{x} \in {\mathbb{R}}^{9}\) (i.e. considering an H-atoms suppressed molecular graph). Table 1 shows the Pauling electronegativity-based molecular vector for Isonicotinic acid. The weighting scheme for the bond-based molecular vector is built with values computed from the properties corresponding to the atoms that each bond connects [10, 13, 20, 40]:

$$w_{ij} = \frac{{w_{i} }}{{\delta_{i} }} + \frac{{w_{j} }}{{\delta_{j} }}$$

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} )\), \(\delta_{i}\) and \(\delta_{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 [4042].

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 bond-based (see Eq. 4) indices, respectively:

$${}^{ns, ss}b^{k} \left( {\bar{x},\bar{y}} \right) = \mathop \sum \limits_{i = 1}^{n} \mathop \sum \limits_{j = 1}^{n} m_{ij}^{k} x^{i} y^{j} = \left( {\bar{x}} \right)^{T} \times {\boldsymbol{\mathcal{M}}}_{ns, ss}^{\varvec{k}} \times \bar{y}\quad \forall k = 1,2, \ldots ,15$$
$${}_{e}^{ns,ss} b^{k} \left( {\bar{x},\bar{y}} \right) = \mathop \sum \limits_{i = 1}^{m} \mathop \sum \limits_{j = 1}^{m} e_{ij}^{k} x^{i} y^{j} = \left( {\bar{x}} \right)^{T} \times {\boldsymbol{\mathcal{E}}}_{ns, ss}^{\varvec{k}} \times \bar{y}\quad \forall k = 1,2, \ldots ,15$$

where, n (or m) is the number of atoms (or bonds) in the molecule, k = 1, 2, …15 is the matrix power, \(m_{ij}^{k}\) (or \(e_{ij}^{k}\)) represents the elements of the \({\boldsymbol{\mathcal{M}}}_{{\varvec{ns},\varvec{ss}}}^{\varvec{k}}\) (or \({\boldsymbol{\mathcal{E}}}_{{\varvec{ns},\varvec{ss}}}^{\varvec{k}}\)) non-stochastic (ns) and simple stochastic (ss) matrices, and \(x^{i}\) and \(y^{j}\) are the elements of the \(\bar{x}\) and \(\bar{y}\) atom-based (or bond-based) property vectors. On one hand, when the vectors \(\bar{x}\) and \(\bar{y}\) encode the same atomic property (i.e. \(\bar{x} = \bar{y}\)), the Eqs. 3 and 4 define the NS and SS total atom-based and bond-based quadratic indices, respectively. On the other hand, if \(\bar{x}\) is a vector with all entries equal to 1 and \(\bar{y}\) 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}{\boldsymbol{\mathcal{M}}}^{k}\) \(({}_{F}{\boldsymbol{\mathcal{E}}}^{k} )\), which is computed from the corresponding kth total matrix \({\boldsymbol{\mathcal{M}}}^{k}\) (\({\boldsymbol{\mathcal{E}}}^{\varvec{k}}\)) considering only those vertices (or edges) belonging to the selected molecular fragment. These fragments F may be heteroatoms (X), halogens (G) and H-bond donors (N or O atoms sharing a bond with an H-atom, labeled as D) [10, 34, 36]. Thus, NS and SS local-fragment atom/bond-based bilinear, quadratic and linear indices can be computed using the \({}_{F}{\boldsymbol{\mathcal{M}}}^{k}\) and \(_{F} {\boldsymbol{\mathcal{E}}}^{k}\) local-fragment matrices instead of the corresponding total matrices in the Eqs. 3 and 4.

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 \({\boldsymbol{\mathcal{M}}}_{{\varvec{ns},\varvec{ss}}}^{k}\) (or \({\boldsymbol{\mathcal{E}}}_{{\varvec{ns},\varvec{ss}}}^{k}\)) 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/bond-based algebraic indices have been computed as whole-molecule (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 \({\boldsymbol{\mathcal{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 \({\mathbb{R}}^{n}\), in a canonical basis set, and whose values are components (entries) of the vector \({\boldsymbol{\mathcal{L}}}\) denoted as \({\boldsymbol{\mathcal{L}}}_{a}\) and \({\boldsymbol{\mathcal{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:

$${{}_{b}}{\boldsymbol{\mathcal{L}}}_{a} = b^{a,k} \left({\bar{x},\bar{y}} \right) = \mathop \sum \limits_{i = 1}^{n} \mathop \sum \limits_{j = 1}^{n} \fancyscript{m}_{ij}^{a,k} x^{i} y^{j} = \left({\bar{x}} \right)^{T} \times {\boldsymbol{\mathcal{M}}}^{a,k} \times \bar{y}\quad \forall a = 1,2, \ldots,n$$
$${}_{b}{\boldsymbol{\mathcal{L}}}_{e} = b^{e,k} \left( {\bar{x},\bar{y}} \right) = \mathop \sum \limits_{i = 1}^{m} \mathop \sum \limits_{j = 1}^{m} {\text{e}}_{ij}^{e,k} x^{i} y^{j} = \left( {\bar{x}} \right)^{T} \times {\boldsymbol{\mathcal{E}}}^{e,k} \times \bar{y}\quad \forall e = 1,2, \ldots ,m$$

where x 1, …, x n(m) and y 1, …, y n(m) are the coordinates or components of the molecular vectors \(\bar{x}\) and \(\bar{y}\) [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 \(\bar{x} = \bar{y}\) atom- and bond-level quadratic indices are obtained [i.e. \({}_{q}{\boldsymbol{\mathcal{L}}}_{a} = q^{a,k} \left( {\bar{x},\bar{x}} \right)\) and \({}_{q}{\boldsymbol{\mathcal{L}}}_{e} = q^{e,k} \left( {\bar{x},\bar{x}} \right)\)], while if all coefficients of \(\bar{x}\) are equal to 1 then linear indices for atoms (or bonds) may be obtained [i.e. \({}_{f}{\boldsymbol{\mathcal{L}}}_{a} = f^{a,k} \left( {\bar{u},\bar{y}} \right)\) and \({}_{f}{\boldsymbol{\mathcal{L}}}_{e} = f^{e,k} \left( {\bar{u},\bar{y}} \right)\)].

The coefficients \({\fancyscript{m}}\) a,k ij (see Eq. 5) are the elements corresponding to the kth NS (or SS) total atom-level pseudograph-theoretic electronic-density matrix [NS(SS)-GEDM] \({\boldsymbol{\mathcal{M}}}^{a,k}\) for atom “a”, while the entries \(e_{ij}^{e,k}\) (see Eq. 6) belonging to kth NS (or SS) total bond-level edge-adjacency matrix [NS(SS)-EAM] \({\boldsymbol{\mathcal{E}}}^{e,k}\) for bond “e”. These atom/bond-level coefficients are obtained from the entries \({\fancyscript{m}}\) k ij of the \({\boldsymbol{\mathcal{M}}}^{k}\) total matrix and \(e_{ij}^{k}\) of the \({\boldsymbol{\mathcal{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 \({\fancyscript{m}}\) 1 ii 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_{ii}^{1}\) 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} {\boldsymbol{\mathcal{M}}}^{a,k}\) for atom “a” and the kth NS (or SS) local-fragment bond-level edge-adjacency matrices \(_{F} {\boldsymbol{\mathcal{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} {\boldsymbol{\mathcal{M}}}^{a,k}\) and \(_{F} {\boldsymbol{\mathcal{E}}}^{a,k}\) as matrix forms, respectively. Note that the coefficients \(_{F} \fancyscript{m}_{ij}^{a,k} \in\,{_{F} {\boldsymbol{\mathcal{M}}}^{a,k}}\) and \(_{F} e_{ij}^{e,k} \in\,{_{F} {\boldsymbol{\mathcal{E}}}^{e,k}}\) are calculated from the elements \(_{F} \fancyscript{m}_{ij}^{k} \in\,{_{F} {\boldsymbol{\mathcal{M}}}^{k}}\) and \(_{F} e_{ij}^{k} \in\,{_{F} {\boldsymbol{\mathcal{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 electronic-density (DS-GEDM, \(_{(F)} {\boldsymbol{\mathcal{M}}}_{{\varvec{ds}}}^{k}\)) and edge-adjacency (DS-EAM, \(_{(F)} {\boldsymbol{\mathcal{E}}}_{{\varvec{ds}}}^{k}\)) matrix approaches can be calculated from the corresponding \({\boldsymbol{\mathcal{M}}}_{ns}^{\varvec{k}}\) and \({\boldsymbol{\mathcal{E}}}_{ns}^{\varvec{k}}\) 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)} {\boldsymbol{\mathcal{M}}}_{mp}^{k}\)) and edge-adjacency matrix (MP-EAM, \(_{(F)} {\boldsymbol{\mathcal{E}}}_{{\varvec{mp}}}^{k}\)) 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.

Schema 1
scheme 1

The stages involved in the computation of the NS-, SS-, DS-, and MP-pseudograph-theoretical electronic-density matrices

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
Table 3 The zero, first and second powers of the molecular pseudograph’s atom adjacency double stochastic and mutual probability matrices for Isonicotinic Acid

Lastly, in order to obtain the global kth total (or local-fragment) bilinear, quadratic and linear indices from the corresponding atom-level (\({\boldsymbol{\mathcal{L}}}_{a}\)) or bond-level (\({\boldsymbol{\mathcal{L}}}_{e}\)) definitions, the summation operator is used. The global indices obtained using this operator over components of vector \({\boldsymbol{\mathcal{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 \({\boldsymbol{\mathcal{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 [2428] 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)} {\boldsymbol{\mathcal{M}}}^{k}\) and \(_{(F)} {\boldsymbol{\mathcal{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 \(_{(F)} {\boldsymbol{\mathcal{M}}}^{k}\) and \(_{(F)} {\boldsymbol{\mathcal{E}}}^{k}\) yields the following representations: the Self-Returning Walks matrices (i.e. \({}_{{\left( \varvec{F} \right)}}^{srw} {\boldsymbol{\mathcal{M}}}^{k}\) and \({}_{{\left( \varvec{F} \right)}}^{srw} {\boldsymbol{\mathcal{E}}}^{k}\)), the non-Self-Returning Walks matrices (i.e. \({}_{{\left( \varvec{F} \right)}}^{nsrw} {\boldsymbol{\mathcal{M}}}_{{}}^{k}\) and \(\varvec{ }{}_{{\left( \varvec{F} \right)}}^{nsrw} {\boldsymbol{\mathcal{E}}}^{k}\)), and the topological path cut-off matrices (i.e. \({}_{{\left( \varvec{F} \right)}}^{p} {\boldsymbol{\mathcal{M}}}^{k}\) and \({}_{{\varvec{ }\left( \varvec{F} \right)}}^{p} {\boldsymbol{\mathcal{E}}}^{k}\)), respectively. The coefficients \({}_{{\left( \varvec{F} \right)}}^{\varvec{p}} m^{1}\) and \({}_{{\left( \varvec{F} \right)}}^{\varvec{p}} e^{1}\) belonging to these last matrices, respectively, are defined as follows:

$${}_{{\left(\varvec{F} \right)}}^{\varvec{p}} \fancyscript{m}^{k} \left[{{}_{{\left(\varvec{F} \right)}}^{\varvec{p}} e^{k}} \right] = \left\{{\begin{array}{*{20}l} {_{{(\varvec{F})}} \fancyscript{m}^{k} \left[{_{{(\varvec{F})}} e^{k}} \right]} \hfill &\quad if\,p_{min}\,\le\,p_{ij}\,\le\,p_{max} \hfill \\ 0 \hfill &\quad {otherwise} \hfill \\ \end{array}} \right.$$

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 indices could be calculated for atoms separated by 1 step (covalent interactions) or for those atoms separated by more than 1 step (\(p \ge 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.

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 [2–5])
Schema 2
scheme 2

Workflow followed in the computation of the ToMoCoMD-CARDD QuBiLS-MAS MDs

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 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/bond-level 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 different parameters used, such as: algebraic forms, matrix approaches, atomic properties, topological cut-offs 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.

Fig. 1
figure 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)

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) 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 tab-separated 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.

Table 5 Comparison between the old software (TOMOCOMD) and the new one proposed in this report (QuBiLS-MAS)
Table 6 Main features of commonly used tools for molecular descriptors (MDs) calculations

Assessment of the performance of the QuBiLS-MAS descriptors

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 ( comprised by 1963 structures was used. The highest SE for this dataset is equal to 10.93 bits (log2N, 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.

Fig. 2
figure 2

In-house comparison of Shannon’s entropy distribution for the QuBiLS-MAS 2D-Indices considering the non-stochastic, simple stochastic, double-stochastic and mutual probability matrix formalisms

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.

Fig. 3
figure 3

In-house comparison of Shannon’s entropy distribution for the QuBiLS-MAS 2D-Indices considering the norms, the statistical operators of central tendency and the operators for dispersion and form

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-Randic_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 above 9.55 bits (87% of the maximum entropy). As for the DRAGON MD families, the 2D-edge_walk, 3D-GETAWAY and 2D-conn_autocorr_inf indices show the best behavior with 63, 21 and 15 variables presenting SE values greater than 8.70 bits (80% of the maximum entropy), respectively, although all these distributions are inferior to the one corresponding to the QuBiLS-MAS 2D-indices. This is a promising result bearing in mind that the DRAGON MD families are obtained from a diverse range of theoretical and practical considerations, encompassing over 30 years of research.

Fig. 4
figure 4

Shannon’s entropy distribution for DRAGON MDs families versus bilinear, linear and quadratic QuBiLS-MAS 2D-Indices

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.

Fig. 5
figure 5

Shannon’s entropy distribution for QuBiLS-MAS topological indices and other descriptors computed by well-known software used in cheminformatics studies

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) variables loaded in each component are linearly independent to the ones loaded in the remaining components. For all the studies performed in this section, the curated Spectrum Collection dataset (1963 molecules) was employed.

To perform this analysis, two sets of descriptors were calculated using QuBiLS-MAS MDs and the DRAGON (824 MDs) software, respectively, with the latter comprising of the following families: 0D-others (B01 Constitutional, B19 Charge and B20 Molecular Properties) with 91 indices, 1D-fragment (B17 Functional Groups Counts and B18 Atom-centered Fragments) with 274 indices, 2D-conn_autocorr_inf (B04 Connectivity, B05 Information and B06 2D-AutoCorrelations) with 176 indices, 2D-edge_walk (B03 Walk-Path Counts and B07 Edge Adjacency) with 154 indices, 2D-eigenvalues (B08 Burden, B10 Eigenvalue-based and B09 Topological Charge) with 129 indices, and finally the B02 2D Topological with 119 indices.

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 2loo (“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_{boot}^{2} )\) 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 2loo  = 81.85%, Q 2boot  = 77.89%) presents the best behavior, followed by MP (Q 2loo  = 79.05%, Q 2boot  = 74.85%). The indices based on NS (Q 2loo  = 73.48%, Q 2boot  = 68.09%) and DS (Q 2loo  = 72.01%, Q 2boot  = 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.

Fig. 6
figure 6

Comparison of the performance of some inner features of the QuBiLS-MAS software in QSAR modeling: a the matrix formalisms, b the aggregation operators and c the classical algorithms

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 2loo 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 2loo and Q 2boot . 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 2loo statistic.

Table 7 Statistical parameters for the best models for 2–6 variables for the physicochemical property log K, considering the 31 structures as the training set
Table 8 Comparison of Q 2loo statistics of nD-QSAR methods for the property log K (CGB) for 31 (or 30)

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 2loo statistic. Up to now, the best model reported has been the one based on the “Combined Electrostatic and Shape Similarity Matrix” (Q 2loo  = 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 2loo  = 0.937 (compound 31 excluded), var = 6), (Q 2loo  = 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.


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 (multi-linear forms) are to be introduced, in addition to a new set of multi-metrics based cut-offs.


  1. Todeschini R, Consonni V (2009) Molecular descriptors for chemoinformatics. In: Mannhold R, Kubinyi H, Folkers G (2009) Methods and principles in medicinal chemistry, Second, Revised and Enlarged ed. vol 1. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, p 2125

  2. Brown FK (1998) Chapter 35. Chemoinformatics: what is it and how does it impact drug discovery. In: James AB (ed) Annual reports in medicinal chemistry. Academic Press, New York, pp 375–384

    Google Scholar 

  3. Todeschini R et al (2006) DRAGON software: an easy approach to molecular descriptor calculations. MATCH Commun Math Comput Chem 56(2):237–248

    Google Scholar 

  4. Hong H et al (2008) Mold2, molecular descriptors from 2D structures for chemoinformatics and toxicoinformatics. J Chem Inf Comput Sci 48(7):1337–1344

    Article  CAS  Google Scholar 

  5. García-Jacas CR et al (2014) QuBiLS-MIDAS: a parallel free-software for molecular descriptors computation based on multilinear algebraic maps. J Comput Chem 35(18):1395–1409

    Article  Google Scholar 

  6. García-Jacas CR et al (2015) Multi-server approach for high-throughput molecular descriptors calculation based on multi-linear algebraic maps. Mol Inform 34(1):60–69

    Article  Google Scholar 

  7. Yap CW (2011) PaDEL-descriptor: an open source software to calculate molecular descriptors and fingerprints. J Comput Chem 32(7):1466–1474

    Article  CAS  Google Scholar 

  8. Cao D-S et al (2013) ChemoPy: freely available python package for computational biology and chemoinformatics. Bioinformatics 29(8):1092–1094

    Article  CAS  Google Scholar 

  9. Steinbeck C et al (2006) Recent developments of the chemistry development kit (CDK)—an open-source java library for chemo- and bioinformatics. Curr Pharm Des 12(17):2111–2120

    Article  CAS  Google Scholar 

  10. Marrero-Ponce Y et al (2006) Bond-Based global and local (bond and bond-type) quadratic indices and their applications to computer-aided molecular design. 1. QSPR studies of octane isomers. J Comput Aided Mol Des 20(10–11):685–701

    Article  CAS  Google Scholar 

  11. Castillo-Garit JA, Marrero-Ponce Y, Torrens F (2006) Atom-based 3D-chiral quadratic indices. Part 2: prediction of the corticosteroid-binding globulinbinding affinity of the 31 benchmark steroids data set. Bioorg Med Chem 14(7):2398–2408

    Article  CAS  Google Scholar 

  12. Marrero-Ponce Y et al (2008) Novel 2D TOMOCOMD-CARDD molecular descriptors: atom-based stochastic and non-stochastic bilinear indices and their QSPR applications. J Math Chem 44(3):650–673

    Article  CAS  Google Scholar 

  13. Marrero-Ponce Y et al (2010) Bond-based linear indices of the non-stochastic and stochastic edge-adjacency matrix. 1. Theory and modeling of ChemPhysical properties of organic molecules. Mol Divers 14(4):731–753

    Article  CAS  Google Scholar 

  14. Marrero-Ponce Y et al (2005) Ligand-based virtual screening and in silico design of new antimalarial compounds using nonstochastic and stochastic total and atom-type quadratic maps. J Chem Inf Model 45(4):1082–1100

    Article  CAS  Google Scholar 

  15. Marrero-Ponce Y et al (2006) Predicting antitrichomonal activity: a computational screening using atom-based bilinear indices and experimental proofs. Bioorg Med Chem 14(19):6502–6524

    Article  CAS  Google Scholar 

  16. Meneses-Marcel A et al (2005) A linear discrimination analysis based virtual screening of trichomonacidal lead-like compounds: outcomes of in silico studies supported by experimental results. Bioorg Med Chem Lett 15(17):3838–3843

    Article  CAS  Google Scholar 

  17. Montero-Torres A et al (2005) A novel non-stochastic quadratic fingerprints-based approach for the ‘in silico’ discovery of new antitrypanosomal compounds. Bioorg Med Chem 13(22):6264–6275

    Article  CAS  Google Scholar 

  18. Marrero-Ponce Y, Huesca-Guillén A, Ibarra-Velarde F (2005) Quadratic indices of the molecular pseudograph’s atom adjacency matrix and their stochastic forms: a novel approach for virtual screening and in silico discovery of new lead paramphistomicide drugs-like compounds. J Mol Struct 717(1–3):67–79

    Article  CAS  Google Scholar 

  19. Marrero-Ponce Y et al (2005) Atom, atom-type, and total nonstochastic and stochastic quadratic fingerprints: a promising approach for modeling of antibacterial activity. Bioorg Med Chem 13(8):2881–2899

    Article  CAS  Google Scholar 

  20. Casanola-Martın GM et al (2007) TOMOCOMD-CARDD descriptors-based virtual screening of tyrosinase inhibitors: evaluation of different classification model combinations using bond-based linear indices. Bioorg Med Chem 15(3):1483–1503

    Article  Google Scholar 

  21. Casañola-Martín GM et al (2006) New tyrosinase inhibitors selected by atomic linear indices-based classification models. Bioorg Med Chem 16(2):324–330

    Article  Google Scholar 

  22. Castillo-Garit JA et al (2008) Estimation of ADME properties in drug discovery: predicting Caco-2 cell permeability using atom-based stochastic and non-stochastic linear indices. J Pharm Sci 97(5):1946–1976

    Article  CAS  Google Scholar 

  23. Marrero-Ponce Y et al (2003) Total and local quadratic indices of the “molecular pseudograph’s atom adjacency matrix”. Application to prediction of Caco-2 permeability of drugs. Int J Mol Sci 4(8):512–536

    Article  Google Scholar 

  24. Barigye SJ et al (2013) Shannon’s mutual, conditional and joint entropy information indices: generalization of global indices defined from local vertex invariants. Curr Comput Aided Drug Des 9(2):164–183

    Article  CAS  Google Scholar 

  25. Barigye SJ et al (2013) Relations frequency hypermatrices in mutual, conditional and joint entropy-based information indices. J Comput Chem 34:259–274

    Article  CAS  Google Scholar 

  26. Marrero-Ponce Y et al (2012) Derivatives in discrete mathematics: a novel graph-theoretical invariant for generating new 2/3D molecular descriptors. I. Theory and QSPR application. J Comput Aided Mol Des 26(11):1229–1246

    Article  CAS  Google Scholar 

  27. Marrero-Ponce Y et al (2015) Optimum search strategies or novel 3D molecular descriptors: is there a stalemate? Curr Bioinform 10(5):533–564

    Article  CAS  Google Scholar 

  28. Garcia-Jacas CR et al (2014) N-linear algebraic maps for chemical structure codification: a suitable generalization for atom-pair approaches? Curr Drug Metab 15(4):441–469

    Article  CAS  Google Scholar 

  29. García-Jacas CR et al (2016) Examining the predictive accuracy of the novel 3D N-linear algebraic molecular codifications on benchmark datasets. J Cheminform 8(10):1–16

    Google Scholar 

  30. García-Jacas CR et al (2016) N-tuple topological/geometric cutoffs for 3D N-linear algebraic molecular codifications: variability, linear independence and QSAR analysis. SAR QSAR Environ Res 27(12):949–975

    Article  Google Scholar 

  31. Marrero-Ponce Y (2003) Total and local quadratic indices of the molecular pseudograph’s atom adjacency matrix: applications to the prediction of physical properties of organic compounds. Molecules 8(9):687–726

    Article  Google Scholar 

  32. Marrero-Ponce Y (2004) Linear Indices of the “molecular pseudograph’s atom adjacency matrix”: definition, significance-interpretation, and application to QSAR analysis of flavone derivatives as HIV-1 integrase inhibitors. J Chem Inf Comput Sci 44(6):2010–2026

    Article  CAS  Google Scholar 

  33. Marrero-Ponce Y et al (2004) Atom, atom-type, and total linear indices of the “molecular pseudograph’s atom adjacency matrix”: application to QSPR/QSAR studies of organic compounds. Molecules 9(12):1100–1123

    Article  Google Scholar 

  34. Marrero Ponce Y (2004) Total and local (atom and atom type) molecular quadratic indices: significance interpretation, comparison to other molecular descriptors, and QSPR/QSAR applications. Bioorg Med Chem 12(24):6351–6369

    Article  CAS  Google Scholar 

  35. Marrero-Ponce Y et al (2004) Tomocomd-Cardd, a novel approach for computer-aided ‘rational’ drug design: I. Theoretical and experimental assessment of a promising method for computational screening and in silico design of new anthelmintic compounds. J Comput Aided Mol Des 18(10):615–634

    Article  CAS  Google Scholar 

  36. Marrero-Ponce Y et al (2005) Atom, atom-type and total molecular linear indices as a promising approach for bioorganic and medicinal chemistry: theoretical and experimental assessment of a novel method for virtual screening and rational design of new lead anthelmintic. Bioorg Med Chem 13(4):1005–1020

    Article  CAS  Google Scholar 

  37. Todeschini R, Consonni V (2000) Handbook of molecular descriptors. In: Mannhold R, Kubinyi H, Timmerman H (eds) Methods and principles in medicinal chemistry, vol 11, 1st edn. WILEY-VCH Verlag GmbH, Weinheim, p 667

    Google Scholar 

  38. Estrada E, Molina E (2001) Novel local (fragment-based) topological molecular descriptors for QSPR/QSAR and molecular design. J Mol Graph Model 20(1):54–64

    Article  CAS  Google Scholar 

  39. Marrero-Ponce Y et al (2005) Non-stochastic and stochastic linear indices of the molecular pseudographs atom adjacency matrix: application to in silico studies for the rational discovery of new antimalarial compounds. Bioorg Med Chem 13(4):1293–1304

    Article  CAS  Google Scholar 

  40. Castillo-Garit JA et al (2008) Bond-based 3D-chiral linear indices: theory and QSAR applications to central chirality codification. J Comput Chem 29(15):2500–2512

    Article  CAS  Google Scholar 

  41. Marrero-Ponce Y et al (2008) 3D-chiral (2.5) atom-based TOMOCOMD-CARDD descriptors: theory and QSAR applications to central chirality codification. J Math Chem 44(3):755–786

    Article  CAS  Google Scholar 

  42. Marrero-Ponce Y et al (2006) Non-stochastic and stochastic linear indices of the molecular pseudograph’s atom-adjacency matrix: a novel approach for computational in silico screening and “rational” selection of new lead antibacterial agents. J Mol Model 12(3):255–271

    Article  CAS  Google Scholar 

  43. Castillo-Garit JA et al (2007) Atom-based stochastic and non-stochastic 3D-chiral bilinear indices and their applications to central chirality codification. J Mol Graph Model 26(1):32–47

    Article  CAS  Google Scholar 

  44. Castillo-Garit JA et al (2008) Atom-based non-stochastic and stochastic bilinear indices: application to QSPR/QSAR studies of organic compounds. Chem Phys Lett 464(1–3):107–112

    Article  CAS  Google Scholar 

  45. Axler SJ (2015) Linear algebra done right. In: Axler S, Ribet K (eds) Undergraduate texts in mathematics, vol 2, 3rd edn. Springer, New York

    Google Scholar 

  46. Sinkhorn R (1964) A relationship between arbitrary positive matrices and doubly stochastic matrices. Ann Math Stat 35(2):876–879

    Article  Google Scholar 

  47. Fourches D, Muratov E, Tropsha A (2010) Trust, but verify: on the importance of chemical structure curation in cheminformatics and QSAR modeling research. J Chem Inf Comput Sci 50(7):1189–1204

    Article  CAS  Google Scholar 

  48. Marrero-Ponce Y, Romero V (2002) TOMO-COMD (TOpological MOlecular COMputer Design) for Windows version 1.0. In: Preliminary version, may be obtained by email request to Marrero-Ponce ( Central University of Las Villas, Santa Clara

  49. Urias RP et al (2015) IMMAN: free software for information theory-based chemometric analysis. Mol Divers 19(2):305–319

    Article  CAS  Google Scholar 

  50. Gutiérrez Y, Estrada E (2002–2004) MODESLAB, v1.5 (MOlecular DEScriptors LABoratory) for windows. Universidad de Santiago de Compostela, España

  51. Georg H (2008) BlueDesc-molecular descriptor calculator. University of Tübingen, Tübingen

    Google Scholar 

  52. Liu J et al (2005) PowerMV: a software environment for molecular viewing, descriptor generation, data analysis and hit evaluation. J Chem Inf Model 45:515–522

    Article  CAS  Google Scholar 

  53. Massey WF (1965) Principal components regression in exploratory statistical research. J Am Stat Assoc 60(309):234–256

    Article  Google Scholar 

  54. Mardia KV, Kent JT, Bibby JM (1979) Multivariate analysis. Academic Press, London

    Google Scholar 

  55. Todeschini R et al (2003) MobyDigs: software for regression and classification models by genetic algorithms. In: Leardi R (ed) Data handling in science and technology. Elsevier, Amsterdam, pp 141–167

    Google Scholar 

  56. Wu CFJ (1986) Jackknife, bootstrap and other resampling methods in regression analysis. Ann Stat 14(4):1261–1295

    Article  Google Scholar 

  57. Lindgren F et al (1996) Model validation by permutation tests: applications to variable selection. J Chemom 10(5–6):521–532

    Article  CAS  Google Scholar 

  58. Wolpert DH, Macready WG (1997) No free lunch theorems for optimization. IEEE Trans Evol Comput 1(1):67–82

    Article  Google Scholar 

  59. So SS, Karplus M (1997) Three-dimensional quantitative structure-activity relationships from molecular similarity matrices and genetic neural networks. 1. Method and validations. J Med Chem 40(26):4347–4359

    Article  CAS  Google Scholar 

  60. Amat L, Besalu E, Carbo-Dorca R (2001) Identification of active molecular sites using quantum-self-similarity measures. J Chem Inf Comput Sci 41(4):978–991

    Article  CAS  Google Scholar 

  61. Shu-Shen L, Chun-Sheng L, Lian-Sheng W (2002) Combined MEDV-GA-MLR method for QSAR of three panels of steroids, dipeptides, and COX-2 inhibitors. J Chem Inf Comput Sci 42(3):749–756

    Article  Google Scholar 

  62. Beger RD, Harris SH, Xie Q (2004) Models of steroid binding based on the minimum deviation of structurally assigned 13C NMR spectra analysis (MiDSASA). J Chem Inf Comput Sci 44(4):1489–1496

    Article  CAS  Google Scholar 

  63. Polanski J (1997) The receptor-like neural network for modeling corticosteroid and testosterone binding globulins. J Chem Inf Comput Sci 37(3):553–561

    Article  CAS  Google Scholar 

  64. Robert D, Amat L, Carbo-Dorca R (1999) Three-dimensional quantitative–activity relationships from tuned molecular quantum similarity measures: prediction of the corticosteroid-binding globulin binding affinity for a steroid family. J Chem Inf Comput Sci 39(2):333–344

    Article  CAS  Google Scholar 

  65. Parretti MF et al (1997) Alignment of molecules by the Monte Carlo optimization of molecular similarity indices. J Comput Chem 18(11):1344–1353

    Article  CAS  Google Scholar 

  66. Silverman BD, Platt DE (1996) Comparative molecular moment analysis (CoMMA): 3D-QSAR without molecular superposition. J Med Chem 39(11):2129–2140

    Article  CAS  Google Scholar 

  67. Tuppurainen K et al (2004) Ligand intramolecular motions in ligand-protein interaction: ALPHA, a novel dynamic descriptor and a QSAR study with extended steroid benchmark dataset. J. Comput Aided Mol Des 18(3):175–187

    Article  CAS  Google Scholar 

  68. Tuppurainen K et al (2002) Evaluation of a novel electronic eigenvalue (EEVA) molecular descriptor for QSAR/QSPR studies: validation using a benchmark steroid data set. J Chem Inf Comput Sci 42(3):607–613

    Article  CAS  Google Scholar 

  69. Polanski J, Bak A (2003) Modeling steric and electronic effects in 3D- and 4D-QSAR schemes: predicting benzoic pKa values and steroid CBG binding affinities. J Chem Inf Comput Sci 43(6):2081–2092

    Article  CAS  Google Scholar 

  70. De K, Sengupta C, Roy K (2004) QSAR modeling of globulin binding affinity of corticosteroids using AM1 calculations. Bioorg Med Chem 12(12):3323–3332

    Article  CAS  Google Scholar 

  71. Kellogg GE et al (1996) E-state fields: applications to 3D QSAR. J. Comput Aided Mol Des 10(6):513–520

    Article  CAS  Google Scholar 

  72. Beger RD, Wilkes JE (2001) Developing 13C NMR quantitative spectrometric data-activity relationship (QSDAR) models of steroid binding to the corticosteroid binding globulin. J Comput Aided Mol Des 15(7):659–669

    Article  CAS  Google Scholar 

  73. Gregorio CD, Kier LB, Hall LH (1998) QSAR modeling with electrotopological state indices: corticosteroids. J Comput Aided Mol Des 12(6):557–561

    Article  Google Scholar 

  74. Turner DB et al (1999) Evaluation of a novel molecular vibration-based descriptor (EVA) for QSAR studies: 2. Model validation using a benchmark steroid dataset. J Comput Aided Mol Des 13(3):271–296

    Article  CAS  Google Scholar 

  75. Polanski J, Walczak B (2000) The comparative molecular surface analysis (COMSA): a novel tool for molecular design. J Comput Chem 24(5):615–625

    Article  CAS  Google Scholar 

  76. Pastor M et al (2000) GRid-INdependent descriptors (GRIND): a novel class of alignment-independent three-dimensional molecular descriptors. J Med Chem 43(17):3233–3243

    Article  CAS  Google Scholar 

  77. Kubinyi H, Hamprecht FA, Mietzner T (1998) Three-dimensional quantitative similarity–activity relationships (3D QSiAR) from SEAL similarity matrices. J Med Chem 41(14):2553–2564

    Article  CAS  Google Scholar 

  78. Beger RD et al (2002) Comparative structural connectivity spectra analysis (CoSCoSA) models of steroid binding to the corticosteroid binding globulin. J Chem Inf Comput Sci 42(5):1123–1131

    Article  CAS  Google Scholar 

  79. Maw HH, Hall LH (2001) E-state modeling of corticosteroids binding affinity validation of model for small data set. J Chem Inf Comput Sci 41(5):1248–1254

    Article  CAS  Google Scholar 

  80. Marín RM, Aguirre NF, Daza EE (2008) Graph theoretical similarity approach to compare molecular electrostatic potentials. J Chem Inf Model 48(1):109–118

    Article  Google Scholar 

  81. Manchester J, Czerminski R (2008) SAMFA: simplifying molecular description for 3D-QSAR. J Chem Inf Model 48(6):1167–1173

    Article  CAS  Google Scholar 

  82. Silverman BD et al (eds) (1998) Comparative molecular moment analysis (COMMA). In: Kubinyi H, Folkers G, Martin YC (eds) 3D QSAR in drug design, vol 3. Kluwer, Dordrecht, pp 183–196

  83. Good AC, So SS, Richards WG (1993) Structure-activity relationships from molecular similarity matrices. J Med Chem 36(4):433–438

    Article  CAS  Google Scholar 

  84. Wagener M, Sadowski J, Gasteiger J (1995) Autocorrelation of molecular surface properties for modeling corticosteroid binding globulin and cytosolic Ah receptor. J Am Chem Soc 117(29):7769–7775

    Article  CAS  Google Scholar 

Download references

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.


YM-P give thanks to support from USFQ with partial finance of Project ID5400 “Chancellor Grant 2016”. CRGJ acknowledges the support from “Dirección General de Asuntos del Personal Académico” (DGAPA) for the postdoctoral fellowship at “Instituto de Química, Universidad Nacional Autónoma de México (UNAM)” in 2016–2017. Work supported by “Programa de Apoyo a la Investigación y el Posgrado (PAIP) 5000-9163” and “Instituto de Química, UNAM” (KMM).

Authors’ information

Professor Yovani Marrero-Ponce received the BS degree in Pharmaceutical Sciences (summa cum laude) from the Central University of Las Villas (UCLV), Santa Clara, Cuba, in 2001, the M.S. degree in Biochemistry from Medical University “Dr. Serafin Ruiz-de Zarate Ruiz”, Santa Clara, Cuba, in 2004, and the Ph.D. degree in Chemistry from Havana University, Havana City, Cuba, in 2005. After post-doctoral fellowships at the University of Valencia, Spain, he founded the Unit of Computer-Aided Molecular “Biosilico” Discovery and Bioinformatic Research (CAMD-BIR Unit, today is known as CAMD-BIR International Network) as a spin-off of the Department of Pharmacy at UCLV. At present, he is an Full Professor/Research of Molecular Pharmacology and Pharmacotherapy at the Universidad San Francisco de Quito (USFQ), and Head of “Grupo de Medicina Molecular y Traslacional (MeM&T)”, Colegio de Ciencias de la Salud (COCSA), Escuela de Medicina, Edificio de Especialidades Médicas 170157, Pichincha, Ecuador. His research interests include molecular modelling and drug discovery, chem-bio-med-informatics, chemometrics, molecular descriptor, chemogenomics, and mathematical, theoretical and computational chemistry. Scopus Author ID: 55665599200. ResearcherID: H-5724-2011. ResearchGate:, Google scholar:, Facebook:

Availability of data and materials

The QuBiLS-MAS software and the respective user manual are freely available online at

Availability and requirements

Project name: QuBiLs Suite project. Project home page: Operating system(s): Platform independent. Programming language: Java. Other requirements: Java 1.8. License: Open source.

Competing interests

The authors declare that they have no competing interests.


This work was partially supported from the USFQ (Project ID5400 “Chancellor Grant 2016”). Dr. CRGJ was further supported by a specific DGAPA’s postdoctoral fellowship to work at “Instituto de Química”, UNAM.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Yovani Marrero-Ponce.

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.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Valdés-Martiní, J.R., Marrero-Ponce, Y., García-Jacas, C.R. et al. QuBiLS-MAS, open source multi-platform software for atom- and bond-based topological (2D) and chiral (2.5D) algebraic molecular descriptors computations. J Cheminform 9, 35 (2017).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: