Background A major challenge in systems biology is to understand the intricate network of interacting molecules. The complexity in biological systems arises not only from various individual protein molecules but also from their organization into systems with numerous interacting partners. In fact, most cellular processes are carried out by multi-protein complexes, groups of proteins that bind together to perform a specific task. Some proteins form stable complexes, such as the ribosomal complex that consists of more than 50 proteins and three RNA molecules, while other proteins form transient associations and are part of several complexes at different stages of a cellular process. A better understanding of this higher-order organization of proteins into overlapping complexes is an important step towards unveiling functional and evolutionary mechanisms behind biological networks. Data on protein complexes are collected from the study of individual systems, and more recently through high-throughput experiments, such as yeast two-hybrid (Y2H) [1,2] and tandem affinity purification followed by mass spectrometry (TAP/MS) [3,4]. The TAP/MS approach helps pinpoint proteins that interact with a tagged bait protein, either directly or indirectly, and are thus suited to identify multi-protein complexes. In fact, several research groups have systematically applied TAP/MS technology to study protein complexes involved in different signaling pathways [5]. Protein interactions are routinely represented as graphs, with proteins as nodes and interactions as edges (links). Therefore, it is not surprising that analysis of protein interaction networks reach out for a variety of graph-theoretical tools. Following the observation that protein interaction networks display a characteristic power-law like node degree distribution [6], a substantial body of research focused on statistical properties of protein interaction networks [7,8]. In 1999, Hartwell et al. [9] introduced a notion of a functional module, a group of cellular components and their interaction that can be attributed a specific biological function. The authors also suggested the modular organization of molecular interaction networks, where each functional module involves a small number of cellular components and is autonomous, i.e., its interaction with other modules is limited to a few cellular components. Subsequently, this assumption was used in several computational methods to identify protein complexes and functional modules in high-throughput protein interaction networks [10-15]. Some methods [10-13] look for densely connected subgraphs within a protein interaction network, either cliques or "cliquish" components. For example, Spirin et al. [13] use the term functional module to denote groups of proteins which are densely connected within themselves but sparsely connected with the rest of the network. Other methods [14,15] combine protein interaction with other information to identify functional modules, such as signal transduction pathways, that do not necessarily correspond to densely connected regions of the network. In a recent paper, Gagneur et al. applied modular decomposition to elucidate the organization of protein complexes [16]. The basic principle behind modular decomposition is to iteratively identify and contract nodes that are in a certain sense equivalent, until no more equivalent nodes can be found in the graph. A graph is called prime if it cannot be decomposed any further. Only graphs that belong to a very special graph family called cographs can be completely decomposed (that is, the iterative reduction process does not halt with a non-trivial prime graph). While the modular decomposition provides an excellent description of combinatorial variants within a family of complexes, it does not impose any order on the complexes within the family. As such it lacks the description power to represent the dynamics of complex formation, i.e., the manner in which proteins form transient interactions to participate in the complexes within the family. The order imposed on protein complexes within the family is particularly interesting if the family corresponds to a functional module where biological function is achieved through a dynamic formation of protein complexes and the order reflects this formation. In this work, we model a functional module as a union of overlapping dense subnetworks called here functional groups. A functional group is either a maximal clique (typically representing a protein complex) or a set of alternative variants of such complexes/cliques. As components of a larger functional module, functional groups are not assumed to be well separated and can have significant overlaps. Intuitively, if a functional module performs a function that requires a sequence of steps (like in the case of a signaling pathway) then we would like functional groups to be snapshots of protein associations at these steps. We propose a new method for identifying and representing overlapping functional groups in a functional module. Furthermore, if the module corresponds to a dynamic process that requires certain complexes (or more generally functional groups) come into contact in a specific order, our method attempts to discover this order. Our method is motivated by a fundamental result for chordal graphs [17], which states that every chordal graph has the so called clique tree representation. However, not every protein interaction network is chordal and not every functional group is a clique. Therefore, we developed a graph-theoretical framework that enables automatic construction of a tree-like representation, analogous to the clique tree representation, for much broader family of graphs. We call this representation the Tree of Complexes representation. The nodes in the tree are functional groups, and for every protein, the set of functional groups that contain this protein forms a single subtree. The "single subtree" requirement restricts significantly the way in which the nodes of the tree can be interconnected. As a consequence, this representation shows a smooth transition between functional groups and allows for tracking a protein's path through a cascade of functional groups. Therefore, depending on the nature of the network, the representation may be capable of elucidating temporal relations between functional groups. We developed a new method, Complex Overlap Decomposition (COD), that given a protein interaction network identifies its functional groups and constructs the Tree of Complexes representation. Our method requires that the network satisfies certain mathematical properties. We applied the COD method to several protein interaction networks, such as the TNFα/NF-κB signaling pathway and the pheromone signaling pathway. The corresponding subnetworks for all interaction networks are extracted from high throughput experimental data. Our results show that the COD method opens a new avenue for the analysis of protein interaction networks.