An efficient biological pathway layout algorithm combining grid-layout and spring embedder for complicated cellular location information

For biological pathways such as signal transduction pathways, gene regulatory networks, and metabolic pathways, one of the crucial techniques for understanding their characteristics is to use graph visualization. Both publicly [1] and commercially available pathway databases [2] display retrieved pathways in the form of graphs to enable users to understand them easily. Usually, in these databases, a large number of pathways are retrieved with various types of criteria according to biologists' purposes. However, it is laborious to manually draw graphs for each request, and hence, automatic layout algorithms specialized for biological pathways are strongly desired.

Thus far, several types of drawing algorithms have been designed for biological pathways and they have been integrated in biological modeling and/or simulation software, e.g., Cell Illustrator [3, 4], Pajek [5], PATIKA [6, 7], and CADLIVE [8, 9].

Karp and Paley extracted biological topologies such as linear, cyclic, and branching pathways and used them as the backbone of the layout [10]. For chemical reaction networks, Becker and Rojas proposed a method [11] that uses the longest directed cycle as the backbone of the layout to capture the flow of reactions. On the other hand, Wegner and Kummer used recursively extracted small cycles as the backbone of the layout [12], because such cycles are known to participate in important recycling processes. Their work has been implemented as an SBML application [13].

Several biological properties are considered in spring embedder approaches. The use of edge directions and simple positional constraints has been proposed for more general metabolic pathways [14, 15]. In GOlorize [16], an additional attractive force is applied to nodes belonging to the same Gene Ontology class. An SBML layout extension, SBWAutoLayout [17], employs the spring embedder approach as its layout algorithm. Schreiber et al. [18] proposed to a generic layout algorithm where the spring force cost is optimized independently in horizontal and vertical directions. Due to the optimization strategy, the algorithm can handle placement constraints for the biological comprehension such as horizontal or vertical aligning of nodes, non-overlapping of nodes, and keeping the same network motifs to some formation. Spring embedder approaches are very popular in the field of Bioinformatics because of the harmonized appearance of their results. However, Li and Kurata noted that spring embedder approaches are not suitable for generating compact layouts of complex pathways [19]. In addition, such approaches have a difficulty in handling complicated positional constraints such as arranging some nodes only on a tours-shaped region, which corresponds to cellular membranes, e.g., nuclear inner membrane, nuclear outer membrane, and plasma membrane.

A grid layout algorithm for biological networks was first proposed by Li and Kurata, in which nodes of the given graph are mapped to grid points and the locally optimal mapping of nodes in terms of the defined cost function is searched over all possible mappings [19]. The cost function is defined by the weighted sum of several components: node distances weighted according to the graph structure and Manhattan distance [19], edge-edge and node-edge crossings [20], rewarding scores for the aligned nodes possessing the same biological attributes [21], and negative inner product between directions of in-edges and out-edges that induces the traceability of the flows [22]. Because finding optimal mapping is NP-hard [23] even when only edge-edge crossings are considered in the cost function, the basic grid layout algorithm repeatedly updates the layout by moving nodes one by one under a greedy search strategy, and a locally optimal layout is obtained after convergence. For efficient calculation, the cost differences calculated by checking movements of a node to a grid point at the current step are cached for calculating the cost differences at the next step [19]. Swapping the positions of two nodes is additionally considered at update steps for the better local optimum without increasing the time complexity [21], whereas Barsky et al. restricted movements of a node to stochastically selected grid points at update steps [24]. Further, the reduction of time complexity is accomplished by using sweep calculation algorithm [22], which can efficiently calculate edge-edge and node-edge crossings and the Manhattan distance. In addition, grid layout algorithms can deal with complicated positional constraints that are often assumed in biological networks as sub-cellular localization information, and thus, they succeed in generating compact and biologically comprehensible layouts.

We propose a new grid-layout-based spring embedder algorithm that considers the spring force cost as a component of the cost function of a grid-layout algorithm. Hence, the cost function consists of the spring force cost and edge-edge and node-edge crossings. As stated above, the sweep calculation can be used to efficiently count the number of edge-edge and node-edge crossings and calculate the Manhattan distance. However, the sweep calculation cannot be used to calculate the spring force. Therefore, we devise a new caching approach for calculating the spring force and propose a new layout algorithm having the same time complexity.

The remainder of this paper is organized as follows. In the Results and Discussion section, we discuss the performance of the proposed algorithm by comparing it with that of the conventional spring embedder approach on three biological networks. The conclusions of our work are presented in the Conclusions section. Finally, the Methods section describes the procedure of the proposed algorithm and its time complexity.

The first model shown in Figure 1 is a famous signal transduction pathway, apoptosis, which is known to participate in various biological processes such as development, maintenance of tissue homeostasis, and elimination of cancer cells. Malfunctions of apoptosis have been implicated in many forms of human diseases such as neurodegenerative diseases, AIDS, and ischemic stroke. Apoptosis is reportedly caused by various inducers such as chemical compounds, proteins, or removal of NGF. The biochemical pathways of apoptosis are complex and depend on both the cells and the inducers. In particular Fas-induced apoptosis has been studied in detail and its simulation model has been proposed [27]. Fas ligands, which usually exist as trimmers in the extracellular region, bind and activate their receptors by inducing receptor trimerization in the cytoplasm membrane region. Activated receptors recruit adaptor molecules such as Fas-associating protein with death domain (FADD), which recruit procaspase-8 to the receptor complex, where it undergoes autocatalytic activation in the cytoplasm. Activated caspase-8 activates caspase-3 through two pathways. In the complex pathway, caspase-8 cleaves the Bcl-2 interacting protein and its COOH-terminal part translocates to the mitochondria where it triggers the release of cytochrome c. The cytochrome c released from the mitochondria binds to apoptotic protease activating factor-1 (Apaf-1) together with dATP and procaspase-9 and activates caspase-9 in the cytoplasm. Caspase-9 cleaves procaspase-3 and activates caspase-3. In other pathway, caspase-8 cleaves procaspase-3 directly and activates it. In the nucleus, caspase-3 cleaves DNA fragmentation factor (DFF) 45 in a heterodimeric factor of DFF40 and DFF45. The cleaved DFF45 dissociates from DFF40, inducing the oligomerization of DFF40 that has DNase activity. The active DFF40 oligomer causes internucleosomal DNA fragmentation, which is an apoptotic hallmark indicative of chromatin condensation. As stated above, these reaction events are strictly regulated in specific cellular locations, and therefore, the corresponding location information cannot be ignored in the resulting layout. Figure 1(a) clearly shows the regulation of these events in each cellular location, i.e., plasma membrane, cytoplasm, mitochondria, and nucleus. In contrast, Figure 1(b) shows the two different flows by Fas-induced apoptosis; however, it is difficult to capture the location information of each event.

Figure 2 shows the cell fate determination model of two gustatory neurons of C. elegans - ASE left (ASEL) and ASE right (ASER) [28]. These neurons are morphologically bilaterally symmetric but physically asymmetric in function, and their fates are strictly regulated by the double negative feedback loop (DNFL), the main path of which consists of four steps: (i) activation of DIE-1 protein leads to the activation of lsy-6 miRNA in the nucleus; (ii) lsy-6 miRNA is transported from the nucleus and inhibits the translation of cog-1 mRNA into COG-1 protein; (iii) if the COG-1 protein is not suppressed, then it is translocated into the nucleus and activates the transcription of mir-273 miRNA in the nucleus; and (iv) mir-273 miRNA is transported from the nucleus and inhibits the translation of die-1 mRNA; this completes the loop to (i). In a manner similar to apoptosis, these DNFL reaction events are strictly regulated in specific cellular locations, and therefore, the corresponding location information cannot be ignored in the resulting layout. Although Figure 2(b) shows these steps, it is difficult to capture the location information of each step. For instance, most of the proteins and miRNAs in this model, e.g., COG-1 protein, LIM-6 protein, DIE-1 protein, mir-273 miRNA, and lsy-6 miRNA, translocate between the nucleus and the cytoplasm. However, such information cannot be inferred from this figure. In contrast, Figure 2(a) shows the regulations of these steps while keeping the cellular location of each step, i.e., cytoplasm and nucleus.

These differences can be observed more clearly in larger models. Figure 3 shows the responses of endothelial cells to the tumor necrosis factor, with an emphasis on the induction of endothelial leukocyte adhesion molecules with more elements than in the other two models [26]. Since adhesion molecules are usually localized on the plasma membrane, many molecules should be on the plasma membrane domain. Usually, external signals are received by these adhesion molecules and transferred into the molecules located in the cytoplasm. Finally, these signals trigger the translation of mRNAs in the nucleus. From the viewpoint of the density of nodes, Figure 3(b) appears to suitably keep a uniform density of nodes. Unfortunately, in terms of the understanding of cascading events, Figure 3(b) does not provide completely useful information because, due to the lack of location information, it is difficult to interpret the network as the response model by the tumor necrosis factor from the external region to the nucleus via the cytoplasm. On the other hand, as shown in Figure 3(a), our layout requires no difficulty in tracing the flow of biological cascading events in our layout, i.e., a reader can easily interpret the network as the response model of a cell from the external region to the nucleus via the cytoplasm as a signal flow.

We propose a grid-layout-based spring embedder algorithm that exploits the advantages on both methods, i.e., consideration of location information and harmonized layouts. Not only the harmonized appearance of resulting layouts, spring force also contributes the reduction of crossings, which is verified by comparing two cases of grid layout algorithms: (i) without considering distance cost and (ii) considering only spring force. (Section 2 of Additional file 1). Although only spring force is considered as the distance cost in this study, we can incorporate other distance cost such as the Manhattan distance cost in the cost function simultaneously.

In addition, to explicitly consider the reduction of crossings, edge-edge and node-edge crossing costs are included in the cost function. To calculate spring forces among nodes, we proposed an efficient calculation method for spring force cost with O(|V|²·h·w), and to calculate other costs, we employ the sweep calculation [22], which can count the crossings for all the possible movements of a node at once. By applying the proposed algorithm and spring embedder algorithm to three biological networks, we verified that the consideration of location information significantly improves the understandability of a network from a biological viewpoint.

In order to realize better biological pathway layouts, under the framework of grid layout, several useful cost functions were proposed, e.g., (i) rewarding score for aligned nodes in one line with the same attribute and (ii) negative inner product of directions of in-edge and out-edges. Feature (i) is very important for biologists because the nodes in a biological pathway usually have biological attributes, e.g., a node is mRNA, protein, modified protein, or complex of proteins, and they explicitly distinguish these components. Feature (ii) is also very important for biochemists because it helps in understanding the reaction flows of the biological pathway. These cost functions can be easily plugged in to our grid layout algorithm without increasing the time complexity. Furthermore, to obtain a better resulting layout, we can also introduce the swapping operation of nodes at each step of moving a node to a vacant grid point to increase the search space while keeping the time order.

Our proposed algorithm succeeded in realizing the required features for biological pathway layouts; however, several enhancements are still required. For example, usually, the combination and order of some biological reactions can be grouped, and thus, this set of reactions, e.g., phosphorylation and dephosphorylation, and related biological elements, e.g., protein and modified protein, can be considered as subgraphs that consume several grid points. Although in our search strategy the final result might fall into bad local optima, for the better local optimum, we can use simulated annealing or other techniques which enables escape from the bad local optima although more computational time is required for the final result. If we could extend the current grid layout algorithm to allow the movement of multiple fixed structured nodes at once, then the required feature would be realized. Our layout framework assumes that compartments representing sub-cellular localizations are allocated by users in advance and then the layout algorithm is applied, but we also considered dynamic adjustment of sizes and positions of these compartments. In this work, the initial state of compartments are given in advance. For automatically providing their initial state, the following approach can be considered as an example. The size of compartment can be determined by the square root of the number of nodes that localize in the same compartment. For the positions of the compartments, we put pair of compartments in close positions if many edges are bridging them.

In addition, the bending of edges that enables bypassing edge-edge and node-edge crossings has not been considered in the current grid layout algorithms. This could be achieved by considering bends as virtual nodes and handling them in a manner similar to normal nodes in search steps.

Given a graph G = (V, E) and a grid of h rows and w columns, we define a cost function for mappings of nodes to grid points and show an algorithm that finds the mapping of nodes, minimizing the cost function in a greedy manner. The cost function is defined by the weighted sum of four components:

(a) Attraction force F _a (d(P (v), P (u))) between pairs of adjacent nodes v and u in the graph G, where P (v) and P (u) are grid points to which v and u are mapped, respectively, and d(P (v), P (u)) is the distance between two grid points P (v) and P (u).

(b) Repulsion force F _r (d(P (v), P (u))) between any pairs of nodes v and u.

(c) Number of edge-edge crossings ∑_e.f∈EI _e (e, f), where I _e (e, f ) is a binary function that returns 1 if e and f cross with each other and 0 otherwise.

(d) Number of node-edge crossings ∑_u∈V,e∈EI _n (v, e) where I _n (v, e) is a binary function that returns 1 if v and e cross with each other and 0 otherwise.

Formally, the cost function is given by

(1)

where (v) is the set of adjacent nodes of v, and w _a , w _r , w _e , and w _n ∈ R⁺ are weights for the components.

Efficient calculation of spring force

Repulsion force for a node v is given by

where the function P (i) returns the grid point to which i ∈ is mapped. Checking the movement of a node to all the vacant grid points requires |V|·h·w calculations, and hence, O(|V|²·h·w) time is required in total at each step.

Although the above naïve calculation has a higher time complexity than existing grid layout algorithms, we propose an efficient calculation. When v is moved from P (v) to q, the repulsion force for v is given by:

Because the term ∑_u∈VF _r (d(q, P(u))) in the above equation depends on q, but not on v, by calculating ∑_u∈VF _r (d(q, P(u))) for all the vacant points q initially, the calculation of c _r (v) requires a constant time. The term ∑_u∈VF _r (d(q, P(u))) for all the vacant points requires O(|V|·h·w) time, and |V|·h·w movements are considered at each step. Therefore, in total, O(|V|·h·w) time is required at each step to calculate the repulsion force.

For the attraction force, the delta cost Δ_{v, p}induced by the movement of a node v to grid point p can be calculated by considering the attraction force between v and its adjacent nodes. In addition, the movement of a node v influences the delta costs only for v and its adjacent nodes, i.e., the delta costs for its non-adjacent nodes at the previous and current steps are the same. Thus, by using the cached delta costs obtained at the previous step, we can calculate the delta costs efficiently. If v is moved from p to q at the previous step, the delta cost for the movement of v to r can be updated by

and for a node u in (v) to r,