Annotation-based feature extraction from sets of SBMLmodels

Chebi GO SBO depth # of anno depth # of anno depth # of anno 1 0 0 1 5 5 1 0 0 2 0 0 2 14651 13 0 0 13 0 0 14 47 658 14 0 0 14 0 0 15 167 2505 15 0 0 15 0 0 1 6 6 9 6 1 6 0 0 1 6 0 0 avg 9,85913772 avg 5,653648226 avg 5,204734 sum 5729 10882 13012 CC Chebi GO SBO depth # of anno depth # of anno depth # of anno 1 0 0 1 0 0 1 0 0 2 0 0 2 2 3 4 6 2 2 4 3 0 0 3 39 117 3 1 3 4 0 0 4 76 304 4 69 276 5 0 0 5 194 970 5 48 240 6 0 0 6 333 1998 6 147 882 7 0 0 7 168 1176 7 16 112 8 2 16 8 105 840 8 4 32 9 2 18 9 15 135 9 0 0 10 11 110 10 1 10 10 0 0 1 1 4 4 4 1 1 0 0 1 1 0 0 12 10 120 12 0 0 12 0 0 1 3 7 9 1 1 3 0 0 1 3 0 0 14 0 0 14 0 0 14 0 0 15 0 0 15 0 0 15 0 0 1 6 1 1 6 1 6 0 0 1 6 0 0


Introduction
Thanks to standardization efforts in Systems Biology [1], modelers today have access to high-quality, curated models in standard formats. The Systems Biology Markup Language (SBML) [2] is an XML-based standard format to encode models as interactions between biological entities. The emerging networks are furthermore enriched with semantic annotations [3] which link model parts to external knowledge in domain-specific ontologies (bioontologies) [4]. Many SBML models live in open model repositories such as BioModels Database [5], the Physiome Model Repository [6], or JWS Online [7]. These repositories distribute computational models and associated data in standard formats. They support necessary management tasks, including curation, annotation, search, version control, data visualization etc. to different extents.
BioModels Database implements a native, SQL-based search [5]. An alternative search is the ranked model retrieval [8]. Here, models and their annotations are mapped on pre-defined model features (e. g., model organism, author, biological entity), leading to a characteristic term vector for each model. The properties of this vector are numeric values mostly describing term frequency and inverse document frequency (TF-IDF) [9]. The ranking is determined by the comparison of search terms (i. e. provided keywords) with the extracted characteristic term vector per model. Current approaches are solely capable of comparing a set of keywords against an indexed corpus of models and retrieve matching models. In addition, it is possible to create a characteristic term vector directly from a model and, subsequently, query a corpus by example.
For example, a standard search for the keywords "cell cycle" in BioModels Database retrieves all models in the corpus that are relevant to the term "cell cycle". Together, all models returned by this search can be seen as a new, cell cycle focused, model set (or corpus). The same is possible for keywords such as "apoptosis", "calcium oscillation" or "NF-κB". At this point, we end up with different sets of thematically related models. To characterize such a set and, later on, compare them, features describing this specific model set will be helpful. However, it is problematic to identify suitable characteristics for arbitrary or thematically focused sets of models.
In this paper we present four methods for annotationbased feature extraction from arbitrary sets of SBML models. Our methods build on combinations of existing approaches for feature extraction [10][11][12][13]. We exemplify our methods by comparing the characteristic features of thematic sets to the features of arbitrary sets of SBML models. The thematic sets were extracted from BioModels Database and represent the cell cycle, apoptosis, calcium oscillation, and NF-κB. Concepts, i. e. terms in the ontology, were extracted from three major bio-ontologies used to semantically enrich models (GO, ChEBI, SBO). We argue that our methods contribute to the determination of similarity between sets of SBML models. They also provide statistics on the use of ontology terms in SBML models, and on the relation between ontology terms and models.

Bio-ontologies
SBML is an XML format. It uses an RDF scheme to add semantic annotations to model parts [14]. Among the ontologies that are used to enrich SBML models, we chose here the following three ontologies, which we believe are the most relevant in model annotation: An ontology of gene and gene product attributes, the Gene Ontology (GO) [15]; an ontology of chemical entities, the Chemical Entities in BIology (ChEBI) [16]; and an ontology for modeling in biology, the Systems Biology Ontology (SBO) [3].
The GO is proposed and maintained by the Gene Ontology Consortium. It aims at standardizing the representation of gene and gene product attributes across species and databases by a structured, precisely defined, common, controlled vocabulary. GO covers three domains. The most important relationships within each domain are is-a and part-of. Additionally, each concept is linked to other kinds of information, including many gene and protein keyword databases.
ChEBI is an ontology of chemical entities of biological interest. All database entries are is_a linked within the ontology. Chemical classifications of ChEBI are aligned with the classification of chemical processes in the GO, and the majority of chemical processes in GO are defined in terms of the ChEBI entities that participate in them.
The SBO provides a set of controlled vocabularies of terms commonly used in Systems Biology. It consists of seven orthogonal branches. Terms within each branch are linked by standard is_a relationships. Formal ties to SBO have been developed for several representation formats in Systems Biology. SBML elements a , for example, carry an optional sboTerm attribute, which allows for a precise definition of the meaning of encoded model entities and their relationships.

Feature extraction from ontologies
For feature extraction it is important to group similar items and to find categories that represent the content of the objects.
Several techniques to determine similarity use distance measures as a basis. Common techniques are euclidian or cosinus distance in vector space [17] or the editing distance for text [9,[17][18][19]. In the context of this work the techniques to distances in ontologies and tree structures are of significance.
The hierarchical structure of the ontology can be used to determine the (semantic) similarity between objects [17]. A distinction is made between two approaches; the graphtheoretic and information-theoretic approach.
Examples for the graph-theoretic approach are the works of Bernstein et al. [17] and Wang et al. [20]. They describe the traditional approach for distance determination in ontologies using the number of edges between the nodes. The inheritance structure is represented in a directed acyclic graph in which the specialization of objects increases with each level. In such a graph the ontology distance can be described as the shortest path between two nodes. The shorter the distance between two nodes, the more similar they are. The problem with this approach is the assumption that the edges represent uniform distances within a taxonomy; i.e. the semantic connections are of equal weight. Li et al. therefore investigate in [21] how path length, depth and local semantic density influence the quality of the similarity function. They come to the conclusion, that for a semantic knowledge base especially path length and depth are important to get similarity results that compare to the human perception of similarity. The similarity values are used in cluster analysis approaches for hierarchical clustering [22]. Applied to the feature extraction task, we group concepts based on their distance in the ontology graph for one bioontology at a time. The top-down approach starts with a cluster containing all concepts and then splits this cluster into smaller groups. The bottom-up approach starts with clusters only containing one concept. Those clusters are merged into larger clusters.
The most prominent representative of the informationtheoretic approach is Resnik [12,13]. This approach exploits the information content of objects to compare. The more information two objects have in common, the more similar they are. The information content of a concept c is dependent on the concept's probability. The probability p(c) is calculated by the frequency freq(c) of the concept and the count N of all concepts of the ontology. It is formally defined by Resnik [12]: If all concepts in an ontology are subordinate to one item, then this item has the greatest probability of 1, because its classification always applies. However, the smaller the probability of a concept is, the higher is its information content. The information content IC can be calculated by the negative logarithm of the likelihood: For example, the root term of the Gene Ontology summarizes all concepts of the ontology and consequently has an information content of zero. A child concept such as establishment of localization (GO_0051234) that summarizes 1408 concepts has a higher information content of 3.34 and a leaf concept such as natural killer cell mediated cytotoxicity directed against tumor cell target (GO_0002420) has the highest information content of 10.59.
In order to determine the common information content of two objects, one considers the deepest element that classifies both objects together. The information content of this element is the degree of mutual information content.
The Information Content can be used to address the problem of overgeneralization when using parent concepts as representatives for child concepts [23]. The challenge of feature extraction in ontologies is to find summarizing features that do not generalize too strongly. Concepts further up in the ontology are less specific than concepts further down in the ontology and, thus, have less "information content". Counting the number of references of a concept and its successor concepts would rank the general concept always highest, as it has more references. The counting approach does not consider the loss of specificity when moving up the ontology. Trißl et al. propose a similarity-based scoring function where a general concept must be supported by more references to yield a good score of representativeness.
For our work we identified the information-theoretic approach and especially the notion of the information content to be of interest. Furthermore, we considered existing approaches for feature extraction in other areas, such as text classification, and selected the document frequency to be to some extend applicable in extracting a pre-defined number of features from sets of SBML models.
The document frequency describes the number of documents in which a term occurs [10,11]. It is used to reduce a vocabulary by removing to rare or common words, respectively. In text classification, common words are removed, because they are not discriminating for any particular class. Rare words are eliminated because they are considered non-informative for category prediction and not influential in global performance. In our specific application, common concepts from bio-ontologies are kept because they are very convenient as features. The discriminating power of a concept is given by the feature value that is saved for each model. However, rarely used concepts are removed during the feature extraction process.
For example, the Gene Ontology Term mRNA catabolic process (GO_0006402) is referenced in over 40 documents, terms of the branch establishment of localization (GO_0051234) are contained in over 200 documents, while terms of cell killing (GO_0001906) are rarely annotated. While the first two terms could be suitable as features, cell killing is not suitable at all, because only a few annotated documents could be found by this term.

Implementation
As a proof of concept, we implemented the four different methods described in Section "Results and discussion" in a prototype application b . We then tested all methods on seven different model sets, which we extracted from BioModels Database.

Prototype
The prototype implementation incorporates two major technologies. First, ontologies are imported using the OWL API [24] and the JFact [25] reasoner. The Web Ontology Language (OWL) is a specification of the World Wide Web Consortium (W3C) to create, publish and to distribute ontologies based on a formal description language [26]. Most bio-ontologies are available in OWL format.
Second, all relevant information about the models and the ontologies is stored in a graph database [27]. A graph database is well suited for models in SBML structure and ontologies alike. It supports links between ontology concepts and SBML models, and it allows for efficient queries [28]. For evaluation purposes, we imported the ontology concepts and their taxonomic relationships and counted the number of annotations referring from a model to a particular ontology concept. The storage approach has been described in detail in an earlier publication [29].

Test sets
We generated seven different test sets containing SBML models from BioModels Database [30]. Two model sets contain arbitrary models, four model sets have a certain biological focus, and one model set contains the complete BioModels Database (Additional file 1: Table S1).
The cell cycle set (CC) contains only models from the curated branch. This ensures ground truth in model annotation as annotations in the curated branch are manually reviewed [5]. In addition to the cell cycle set, the two random sets (RS1 and RS2), the thematic test sets for apoptosis (APOP), calcium oscillation (CA) and NF-κB (NFKB), and the set containing all 490 curated models (BMDB) were assembled from the curated branch. In contrast to the CC set (containing 30 models) the thematic test sets APOP, CA and NFKB only contain about 13 models each.
Consequently, we rely on the cell cycle set in our analysis of methods, we use the three other thematic sets for evaluation purposes.
The models for all model sets were pre-selected using our previously developed retrieval algorithm [8]. For example, the first test set is a thematic set containing SBML encodings of published cell cycle models. We used the term "cell cycle" for a keyword based search to retrieve a list of relevant models. To exclude possible false positive search results we manually validated the retrieved models based on their reference publications, resulting in the 34 given models for cell cycle. The model sets APOP, CA and NFKB were compiled in the same way.
From the biological point of view, the test sets CC, APOP, and NFKB are thematically similar. NFKB, which is one of the most prominent transcription factors, is able to manipulate cyclins that drive the cell cycle [31] and additionally has stimulus dependent pro-or anti-apoptotic functions [32]. Moreover, the connection between cell cycle and apoptosis is presented by many cells starting their apoptotic cell fate decision from the cell cycle arrest (G1/S checkpoint), i. e. after caspase activation [33]. More recently, calcium oscillations were shown to influence NF-κB activity depending on the calcium spike duration [34]. We deliberately introduced the NFKB set with strong relations to the CC and APOP sets to evaluate if our methods reflect these relations in terms of similarity of extracted features. The assumption is that biologically similar model sets share semantic annotations.

Results and discussion
Our main hypothesis is that it should be possible to extract characteristic features from semantic annotations, both for thematic sets of models and for arbitrary ones. The following subsections explain our four methods for feature identification, based on the aforementioned feature extraction methods (Section "Implemented feature extraction methods"); discuss their applicability to feature extraction from model sets (Section "Applicability of methods"); show the distribution of model annotations in BioModels Database (Section "Distribution of SBO concepts in SBML models"); and discuss the results obtained from two selected methods when applied to the abovementioned test sets (Section "Feature extraction from arbitrary model sets"). We conclude that it is indeed possible to identify characteristic features. These features can, for example, help with model retrieval, comparison and clustering.

Implemented feature extraction methods
Our methods are designed to identify a predefined, maximum number of features for each compiled set of models. All methods incorporate the structure of the underlying ontology when grouping the concepts within it. Parent concepts represent the group containing their child concepts. Consequently, the developed methods are only applicable to taxonomy-shaped ontologies. Method 1 depends only on the chosen ontology, but not on the input set of models. All other methods additionally consider the annotations in the given set of models.
Method 1 is a top-down clustering. To decide on the suitability of a concept for characterization, the probability p of each concept in the ontology is determined, following Resnik's definition (Equation 1). In the context of this work, the frequency freq(c) refers to the number of all concepts that are summarized by a parent concept c.
Method 2 is a top-down clustering that considers both the ontology structure and the annotations used in models of the given set. Consequently, the real distribution of references to ontology concepts used in models is regarded. Selected features depend on the given set of models. For each concept in the ontology, we count the number of annotations that refer to it. We call this number entity frequency. Additionally, we store the sum of a concept's entity frequency and its descendants' entity frequencies as aggregated entity frequency EF. All concepts with EF > 0 provide the basis for feature extraction. Method 2 re-uses the algorithm of Method 1. The algorithm is adjusted to the dynamic setting by using the entity frequency metric instead of the probability p(c). To better compare the balance of the branches, we will normalize EF as entity probability ep(c): Method 3 is a bottom-up clustering relying on the same input as Method 2. It also uses the entity probability ep(c) but begins with the individual concepts, which are gradually merged to form greater clusters. The results of this method are nearly identical to the ones of Method 2, but the performance of Method 2 is much better.
Method 4 is a bottom-up clustering that addresses the problem of overgeneralization. It uses an adaptation of the scoring function as described in [23]: The Score T (c) for a grouping represented by the concept c considers the information content and the aggregated entity frequency. The information content is calculated depending on the probability of c (see Equations 1 and 2). A group is formed by merging concepts with the ancestor that reaches the highest possible score.

Applicability of methods
We tested the applicability of all described methods on sets of SBML models taken from BioModels Database. Method 1 calculates the probability to hit a certain node in an ontology with a model entity. It condenses a given ontology to a defined number of features, based on the probability of a concept in the ontology only. Thus, the results obtained from Method 1 do not depend on the actual ontology concepts that are referenced in the model set. Consequently, it does not adapt to the specifics of the corpus under study. Therefore, Method 1 is only suitable to provide a static set of features, solely based on the underlying ontology. As a result we dismissed Method 1 for the problem of finding characteristics for arbitrary model sets. However, Method 1 calculates the distribution of concepts in bio-ontologies, as shown in Section "Distribution of SBO concepts in SBML models". Method 2 and Method 3 rely on entity probabilities. Our evaluations show that Method 2 (top-down) and Method 3 (bottom-up) produce almost identical results.
The direction is only relevant in the rare constellation that two concepts are subsumed to the same score. In the following, we consider Method 2 for further evaluations. Method 4 is a dynamic approach that calculates the score value by entity frequency and information content. Based on the unique scoring and the absence of splits, Method 4 generally finds fewer features than the prior methods. It also selects more specific features (located further down in the ontology tree) that are still representative for the model sets. In Section "Feature extraction from arbitrary model sets" we use Method 2 and Method 4 to discuss the specificity and distinctness of extracted features.

Distribution of SBO concepts in SBML models
Using Method 1, we compare the distributions of concepts in the SBO with the frequency of annotations as they occur in all models from BioModels Database. It becomes obvious that the concepts are unequally distributed across seven top-level branches (Figure 1, top). This is explained by the design of the SBO and its orthogonal branches. For example, the branch modeling framework (SBO:0000004) lists a "set of assumptions that underlay a mathematical description" whereas the branch mathematical expression (SBO:0000064) contains "formal representation of a calculus linking parameters and variables of a model".
Consequently, one expects more entries for mathematical expression than for modeling framework.   In conjunction with the application of SBO in model annotation, concepts of some branches are annotated more frequently (Figure 1, bottom). For example, the branch physical entity representation (SBO:0000236), which is a "representation of an entity that may participate in an interaction, a process or relationship of significance", contains only 10% of SBO concepts, but 47% of the model annotations link to that branch. We expect that the characteristic features follow the distribution of the model annotations as seen in the lower part of the figure. Indeed, after applying Method 4, the selected SBO features show a distribution (66.6% physical entity representation (SBO:0000236), 6.6% participant role (SBO:00000003), 13.3% occurring entity representation (SBO:0000231), 6.6% mathematical expression (SBO:0000064), and 6.6% systems description parameter (SBO:0000545) that is closer to Figure 1 (bottom) than before; please refer to Table 1, Method 4, SBO, 15 features).
We also investigated for each model set the distribution of the depth of annotated concepts in the ontology tree.
This knowledge helps us to decide on how specific a model annotation is. Figure 2 shows the distribution for model annotations using ChEBI, GO and SBO (Additional file 2).
Here, we plotted the distribution of annotations for the CC and the BMDB sets. As one would expect, both test sets show normal distributions. The 30 models contained in the CC set make up 6% of the 490 models in the BMDB set. However, the number of annotations in the CC set that refer to ChEBI is less than 1% compared to the number of annotations in the BMDB set. It should be considered that very sparsely annotated model set may be inferior in terms of specificity and distinctness. This information helps us later on in Section "Feature extraction from arbitrary model sets" to decide on the value of the extracted features.

Feature extraction from arbitrary model sets
We hypothesize that the vast property space of a set of models can be condensed into a smaller, but still descriptive, number of features. To establish such "characteristic features", we collect the models' annotations and analyze the semantics behind the linked ontology terms. We focus on the semantics behind the model elements because we believe that this information will be most influential. All our methods require setting a maximum number of features.
Here, we chose to run our extraction methods with five and 15 features as an upper limit. The resulting sets of features for all feature extraction algorithms, models, and ontologies are shown in Table 1.
Specificity of selected ontology concepts. Table 1 shows the average depth of concepts in all three ontologies for all identified features in the CC and BMDB sets. Additionally, Figure 2 contains the average depth of annotation for the CC and BMDB sets before applying the feature extraction methods. The data confirms that the average depth of annotations decreases for Methods 2 and 4 (for all three ontologies and both model sets). Thus, selected concepts are higher up in the ontology, and more generic. This behavior is expected as the feature extraction process also involves generalization. However, the features extracted by Method 4 are more specific than the features extracted by Method 2. This is in accordance with the design of Method 4 to prevent overgeneralization. Moreover, the average annotation depth for the CC set is higher than for the corresponding BMDB set. This supports our assumption that thematically similar models share more annotations, and consequently the extracted features are more specific. For example, the concepts that were selected from ChEBI by Method 2 with a maximum of 15 features for the CC set have an average annotation depth of 5.9. In contrast, the concepts that were selected for the BMDB set only have an average depth of 4.8. According to our obtained data we infer that Method 4, in general, provides features that correspond to deeper concepts in the ontology than the features obtained from Method 2. We conclude from our test data that the depth of chosen concepts decreases with the increased randomness in the sets of models. This is not unexpected, as a broader data basis should not be characterizable by very specific ontology concepts. Rather, an arbitrary model set should cover many different semantic concepts, leading to more generic features being extracted. This behavior is also reflected in our data. In summary, both methods extract features that are specific to the model set. However, features extracted by Method 4 are mostly more specific than those extracted by Method 2. An exemption where the average depth slightly increases is Method 4 for SBO and 15 features. SBO is relatively small compared to GO or ChEBI. As Method 2 is required to select 15 features and Method 4 is only required to select up to 15 features, Method 4 selects only the most relevant features whereas Method 2 selects exactly 15 features. Due to the size of SBO, Method 2 adds features that are not best matches, nevertheless have a higher depth within SBO.
This phenomenon did not occur for the lager ontologies, GO and ChEBI.
Distinctness of feature sets. Another important question is how distinct the obtained features are for our test sets. If the methods retrieved similar concepts for the four test sets, then the extracted features could not be regarded specific to the set of models. Consequently, we measure overlap of concepts between the different characteristic features that we calculated with Method 2 and Method 4. Ideally, there would be almost no overlap of features selected for the CC set with any other selected set, whereas an overlap between BMDB and  the random sets is expectable. Our results are shown in Figure 3. A good result is achieved for Method 4 using 15 features and GO. Here, the cell cycle features have almost no overlap. The result achieved for Method 2 using 15 features and GO is not satisfiable. Here, the cell cycle features largely overlap with at least two other sets. However, the Venn diagrams, in general, confirm that both methods determine features that are specific to the model sets. They contain higher numbers of overlapping features at the intersection between arbitrary sets and very few overlapping features at the intersection between the CC and the BMDB sets. This is particularly visible for the results obtained from Method 4.
Similarity of model sets. We are also interested in how characteristic the sets of extracted features are for a given set of models. We first calculate the similarity of two concepts within the same ontology, as described by Li et al. [21]: The variable h is the depth of the least common subsumer of the concepts c 1 and c 2 , and the variable l is the length of the shortest path between both concepts. Following [21], the parameters are set to α = 0.2 and β = 0.6. We calculate this similarity value for each possible combination of features from two sets of models. Afterwards we apply an adaptation of the Hungarian method [35] to the matrix resulting from the above calculations. The Hungarian method, a solution for the assignment problem, aligns pairs of features, in a way that ensures a global maximum similarity. Based on this similarity of features, we then calculate the total similarity of two sets of features, which corresponds to the similarity of the associated sets of models. The results are shown in Table 2.
Firstly, we discuss specificity of extracted features for the cell cycle set versus the set containing all curated models from Biomodels Database, and one random set. Desirable are low similarities for BMDB vs CC as well as CC vs RS1. As CC is a thematic set, its extracted features should differ from the features extracted from the BMDB and arbitrary model sets. A higher similarity is expected for BMDB vs RS1, as both sets represent a wide range of model topics. The results in Table 2 reflect our expectations. Particularly, the similarity values for Method 4 using 15 features clearly distinguish the extracted features of two sets. Method 2 using five features still shows the desired result, but due to the limited number of features the selected ones are more general and not very distinguishable. Even though results of Method 2 show the expected behavior, we conclude that the results of Method 4 are superior.
Secondly, we discuss the specificity for all thematic sets. Here, we narrow our scope to the Gene Ontology. As Table 3 indicates only the number of distinct annotations using GO is sufficient for all four thematic sets. In addition, we manually reviewed the extracted features and deduced that the features extracted for GO have the highest biological meaning. We use the aforementioned approach to calculate similarity between extracted features of six sets (BMDB, RS1, CC, APOP, CA, NFKB), as shown in Tables 1 and 4. Results for five and 15 selected sets are shown in Table 5. It becomes obvious that the similarities for Method 2 are to high in general, this supports our previous assumption of Method 2 over-generalizing the extracted features. An example of over-generalization is Method 2 using 5 features and the sets RS2 and NFKB. Both sets perfectly match. The reason for this match is, that Method 2 selected only top and second level representatives (both sets have an average depth of 1.8).
Desirable are low similarities for each thematic set versus the BMDB, RS1 or RS2 set, respectively. Both result tables show according similarity values for Method 4. We expect NFKB to have a slightly higher similarity to the other three thematic sets as NF-κB has a regulatory effect on cell cycle, apoptosis and calcium oscillation. For five and 15 selected features Method 4 fits our expectation. The relation between CC and APOP is also visible as many cells start apoptosis from the cell cycle arrest. This is also supported by Method 4 for five and 15 features, respectively. In contrast, we predict CA to be distinct from CC and APOP as calcium oscillation has low overlap with cell cycle or apoptosis. Again, Method 4 advocates our prediction. In conclusion, Method 4 was able to support all our assumptions, even if only five characteristic features are provided per set.