Information Flow in Interaction Networks

Aleksandar Stojmirović and Yi-Kuo Yu***to whom correspondence should be
addressed

National Center for Biotechnology Information
National Library of Medicine
National Institutes of Health
Bethesda, MD 20894
United States

Interaction networks, consisting of agents linked by their interactions, are
ubiquitous accross many disciplines of modern science. Many methods of
analysis of interaction networks have been proposed, mainly concentrating on
node degree distribution or aiming to discover clusters of agents that are
very strongly connected between themselves. These methods are principally
based on graph-theory or machine learning.

We present a mathematically simple formalism for modelling context-specific
information propagation in interaction networks based on random walks. The
context is provided by selection of sources and destinations of information
and by use of potential functions that direct the flow towards the
destinations. We also use the concept of dissipation to model the aging of
information as it diffuses from its source.

Using examples from yeast protein-protein interaction networks and some of the
histone acetyltransferases involved in control of transcription, we
demonstrate the utility of the concepts and the mathematical constructs
introduced in this paper.

##  1 Introduction

Interaction networks are abundant and have recently gained significant
publicity in many diverse modern disciplines such as electronics (Cancho _et
al._ , 2001), sociology (Wasserman and Faust, 1994; Newman, 2004) and
epidemiology (Barthelemy _et al._ , 2005). In its simplest form, an
interaction network consists of a collection of entities (or agents), where
two agents are linked if they interact in some way. For example, in an
acquaintance network (Amaral _et al._ , 2000), the agents represent persons
and two persons are linked together if they know each other while the Woldwide
Web network consists of web pages with links between pages (Broder _et al._ ,
2000). Mathematically, networks correspond exactly to graphs (or multigraphs),
with agents as vertices and links as edges, which can be weighted and/or
directed depending on the exact application being modeled. The key to analysis
of interaction networks is the assumption of information transitivity:
information can flow through or can be exchanged via paths of interactions.

Biology in post-genomic era also contains numerous examples of molecular
networks (Galitski, 2004). Metabolic networks have been modeled by
representing metabolites as nodes and chemical reactions as links: two
metabolites are linked if they participate in the same reaction (Ma and Zeng,
2003). Genetic networks have genes as nodes with two genes being linked if
they interact through directed transcriptional regulation (Guelzim _et al._ ,
2002). Protein-protein interaction networks have proteins as nodes, with the
links representing physical interactions (binding) between proteins
(Pellegrini _et al._ , 2004). Large scale high-throughput studies in model
organisms such as Saccharomyces cerevisiae (baker’s yeast) (Ito _et al._ ,
2001; Uetz _et al._ , 2000), Drosophilla melanogaster (fruit-fly) (Giot _et
al._ , 2003), Caenorhabditis elegans (roundworm) (Li _et al._ , 2004) and
humans (Stelzl _et al._ , 2005; Rual _et al._ , 2005), provided extensive
datasets of protein-protein interactions, stored in publicly-available
databases such as the Database of Interacting Proteins (DIP) (Xenarios _et
al._ , 2002; Salwinski _et al._ , 2004). Unfortunately, there is very little
consistency between the protein-protein interaction data coming from different
high-throughput experiments (Sprinzak _et al._ , 2003) and significant effort
has been expended in devising ways to discover false positives and false
negatives (Suthram _et al._ , 2006). This problem is not restricted to
protein-protein interactions: microarray data also contains non-negligible
inconsistencies (Miklos and Maleszka, 2004) .

Numerous approaches have been proposed for analysis of biological and, in
particular, protein-protein interaction networks (Aittokallio and Schwikowski,
2006). However, due to space restrictions, we will refer to just a few. Most
algorithms aim to discover ‘functional modules’ (Hartwell _et al._ , 1999),
representing well connected clusters of nodes with the same or similar
function, by using clustering techniques from graph theory and/or machine
learning (Steffen _et al._ , 2002; Spirin and Mirny, 2003; Rives and Galitski,
2003; Pereira-Leal _et al._ , 2004; Nabieva _et al._ , 2005; Xiong _et al._ ,
2005; Chua _et al._ , 2006; Chen and Yuan, 2006; Hwang _et al._ , 2006). Very
frequently, these techniques make use of additional experimental data which is
not present in the network structure itself. For example, methods for
discovery of complexes from protein-protein interaction networks often refer
to the data from dataset from different species (Kelley _et al._ , 2003;
Sharan _et al._ , 2005a, b), microarray expression studies (Steffen _et al._ ,
2002; Chen and Yuan, 2006), or human-curated functional classifications
(Nabieva _et al._ , 2005; Chua _et al._ , 2006).

Our approach to analyzing interaction networks is very different, relying
solely on the network structure. We model diffusion of information through the
network by discrete-time random walks moving from the nodes representing the
sources of information to their destinations. The choice of sources and
destinations provides the _context of analysis_ with the nodes most affected
by information flow being called _Information Transduction Modules_. We use
two modes of diffusion, dual to each other, which we call absorbing and
emitting, with our absorbing mode directly corresponding to deeply
investigated absorbing Markov chains (Kemeny and Snell, 1976). Random walks
and corresponding Markov chains are one of the subjects of spectral graph
theory (Chung, 1997) but we do not use eigenspace decomposition in our work,
instead relying on a basic matrix algebra approach similar to that of Kemeny
and Snell (1976).

The algorithm Functional Flow by Nabieva _et al._ (2005), also modeling
diffusion of information from sources, is closest to our emitting model.
However, to delineate a certain biological context, we additionally direct the
flow from sources to selected destinations using potential functions and allow
the information content to dissipate (evaporate) from the network at each time
step, thus modeling natural ‘aging’ of information.

Our models allow investigation of several types of biological questions from
protein-protein interaction networks. Many proteins perform their function in
cooperation with other proteins through, often large, protein complexes. Thus,
to elucidate the function of a given protein, it is useful to know the most
likely members of complexes it may belong to and their relations to each
other. Additionally, if two proteins are known to have similar function, what,
if any, are the proteins they share in their respective complexes? To help
answer such questions, we employ our absorbing diffusion mode.

The answers to the above questions can provide the general interaction
environment of one or more proteins. It is also very instructive to identify
specific modules mediating interactions between distant (in network terms)
proteins. Our emitting diffusion mode can be used to find possible candidates
for members of such modules. Furthermore, analysis of interaction modules
obtained from considering different proteins in the same biological context
may lead to discovery of fundamental units of information transduction. To
achieve this we developed the concept of information interference. More
concrete definitions will be presented in the body of the text.

This paper is organized as follows. Section 2 outlines the theory behind our
models of information diffusion in networks. For better readability, all the
theorems and proofs, using mainly the basic concepts and results from the
matrix algebra are given in Appendix (the reader may wish to consult the
standard linear algebra textbooks such as (Hoffman and Kunze, 1971) or (Bapat
and Raghavan, 1997) for background). Section 3 introduces the methods of
analysis of results obtained using the concepts of Section 2, while Section 4
presents concrete examples centered around yeast histone acetyltransferases.
We finish with discussion and conclusion in Section 5.

##  2 Theory

###  2.1 Preliminaries

We represent an interaction network as a weighted directed graph
$\Gamma=(V,E,w)$ where $V$ is a finite set of vertices of size $n$,
$E\subseteq V\times V$ is a set of edges and $w$ is a non-negative real-valued
function on $V\times V$ that is positive on $E$, giving the weight of each
edge (the weight of non-existing edge is defined to be $0$). Assuming an
ordering of vertices in $V$, we represent a real-valued function on $V$ as a
state (column) vector $\mathbf{\boldsymbol{\varphi}}\in\mathbb{R}^{n}$ and the
connectivity of $\Gamma$ by the _weight_ matrix $\mathbf{W}$ where
$W_{ij}=w(i,j)$ (the weight of an edge from $i$ to $j$). If $\Gamma$ is an
unweighted undirected graph, $\mathbf{W}$ is the adjacency matrix of $\Gamma$
where

$W_{ij}=\begin{cases}2&\text{if $i=j$ and $(i,i)\in E$},\\\ 1&\text{if $i\neq
j$ and $(i,j)\in E$},\\\ 0&\text{if $(i,j)\not\in E$.}\end{cases}$ (1)

Throughout this paper, we will not make distinction between a vertex $v\in V$
and its corresponding state given by a particular ordering of vertices.

Let $\mathbf{P}$ denote the $n\times n$ _transition_ matrix of $\Gamma$ where

$P_{ij}=\frac{W_{ij}}{\sum_{k}W_{ik}},$ (2)

that is, $\mathbf{P}$ is the weight matrix of $\Gamma$ normalized by row. The
matrix $\mathbf{P}$ can be used to model random walks on $\Gamma$: for any
pair of vertices $i$ and $j$, $P_{ij}$ gives the probability of the random
walk moving from vertex $i$ to vertex $j$ in one time step, which is
proportional to the weight $W_{ij}$. Since the matrix $\mathbf{P}$ is
stochastic (all rows sum to unity), it can also be interpreted as the
transition matrix for Markov chain on the set $V$. In the following sections
we will model information diffusion as a random walk on $\Gamma$ with
particular starting and terminating points.

###  2.2 Constrained diffusion

In this section we select certain vertices as sources or sinks of information
and solve for the number of times a vertex is visited. Let $S$ denote the set
of selected vertices, let $T=V\setminus S$ and let $m=\left|T\right|$.
Assuming that the first $n-m$ states correspond to vertices in $S$, we write
the matrix $\mathbf{P}$ in the canonical form:

$\mathbf{P}=\left[\begin{array}[]{cc}\mathbf{P}_{SS}&\mathbf{P}_{ST}\\\
\mathbf{P}_{TS}&\mathbf{P}_{TT}\end{array}\right].$ (3)

Here $\mathbf{P}_{AB}$ denotes a matrix giving probabilities of moving from
$A$ to $B$ where $A,B$ stand for either $S$ or $T$. The states (vertices)
belonging to the set $T$ are called _transient_.

####  2.2.1 Absorption in sinks

Suppose now that the set $S$ represents the set of _sinks_ of information: any
information reaching a sink vertex is absorbed and cannot not leave it. Let
$\mathbf{F}(t)$ denote an $m\times(n-m)$ matrix such that $F_{ij}(t)$ is the
probability that the information originating at $i\in T$ is absorbed at $j\in
S$ in $t$ or fewer steps. Since information can only be absorbed once in any
state $s\in S$, it follows that the information reaching $j$ avoided all other
sinks. For the same reason, $F_{ij}(t)$ can be interpreted as the expected
number of visits to the state $j$ of a random walk starting at $i$ for all
times up to $t$.

Absorption at $j$ after not more than $t$ steps can be achieved in two ways:
either the content reached vertex $j$ in the first step, with probability
$P_{ij}$ or it moved to some transient vertex $k$ in the first step and was
absorbed by $j$ from there in at most $t-1$ steps, with probability
$P_{ik}F_{kj}(t-1)$. Therefore, we have for all $t=1,2,\ldots$,

$F_{ij}(t+1)=P_{ij}+\sum_{k\in T}P_{ik}F_{kj}(t),$ (4)

or in the matrix form

$\mathbf{F}(t+1)=\mathbf{P}_{TS}+\mathbf{P}_{TT}\mathbf{F}(t).$ (5)

We solve for the long-term or equilibrium state, where
$\mathbf{F}(t+1)=\mathbf{F}(t)=\mathbf{F}$. In this case, Equation (5) becomes

$\mathbf{F}=\mathbf{P}_{TS}+\mathbf{P}_{TT}\mathbf{F},$ (6)

or

$(\mathbb{I}-\mathbf{P}_{TT})\mathbf{F}=\mathbf{P}_{TS},$ (7)

where $\mathbb{I}$ denotes the identity matrix. If
$\mathbb{I}-\mathbf{P}_{TT}$ is invertible, let
$\mathbf{G}=(\mathbb{I}-\mathbf{P}_{TT})^{-1}$. Equation (7) then has a unique
solution

$\mathbf{F}=\mathbf{G}\mathbf{P}_{TS}.$ (8)

####  2.2.2 Diffusion from sources

Now consider the dual problem where $S$ is a set of sources of information.
Each source emits a unit of information at each time step and no information
can enter any source: we assume any information entering a source vanishes.
Let $\mathbf{H}(t)$ denote an $(n-m)\times m$ matrix such that $H_{ij}(t)$ is
the total expected number of times the transient vertex $j$ is visited by a
random walk emitted from source $i$ for the time up to $t$.

The information emitted from $i$ can arrive at $j$ at time $t$ in two
different ways: either the content was emitted from $i$ at time $t$ and
reached $j$ directly, or it was emitted at an earlier time step, was located
at some transient vertex at time $t-1$ and moved from there to $j$ at time
$t$. The former option contributes $P_{ij}$ while the latter contributes
$H_{ik}(t-1)P_{kj}$ for all $k\in T$ towards $H_{ij}$. Therefore, we have for
all $t=1,2,\ldots$,

$H_{ij}(t+1)=P_{ij}+\sum_{k\in T}H_{ik}(t)P_{kj},$ (9)

or in the matrix form

$\mathbf{H}(t+1)=\mathbf{P}_{ST}+\mathbf{H}(t)\mathbf{P}_{TT}.$ (10)

Similarly to the previous case, we are interested in the steady state,
representing the total expected number of visits, where
$\mathbf{H}(t+1)=\mathbf{H}(t)=\mathbf{H}$. In this case, Equation (10)
becomes

$\mathbf{H}=\mathbf{P}_{ST}+\mathbf{H}\mathbf{P}_{TT},$ (11)

or

$\mathbf{H}(\mathbb{I}-\mathbf{P}_{TT})=\mathbf{P}_{ST}.$ (12)

If $\mathbb{I}-\mathbf{P}_{TT}$ is invertible, Equation (12) has a unique
solution

$\mathbf{H}=\mathbf{P}_{ST}\mathbf{G}.$ (13)

####  2.2.3 Existence and interpretation of solutions

It can immediately be observed that existence of solutions to Equation (12)
and Equation (7) are equivalent: they both depend on the existence of the
inverse of $\mathbb{I}-\mathbf{P}_{TT}$. Specifically, they are special cases
of the discrete Laplace equation on $T$ with the Dirichlet boundary condition
on $S$ (Chung, 1997; Chung and Yau, 2000).

Given a square matrix $\mathbf{M}$, the matrix $\mathbb{I}-\mathbf{M}$ is
often called the _discrete Laplace operator_ of $\mathbf{M}$. Let
$\Delta=\mathbb{I}-\mathbf{P}_{TT}$ ($\Delta$ is the discrete Laplace operator
of $\mathbf{P}$ restricted to $T$). Equation (7) can then be written as

$\Delta\mathbf{F}=\mathbf{P}_{TS}.$ (14)

Denote by $\mathbf{e}_{k}$ the $k$-th standard basis (column) vector of length
$n-m$ where $(\mathbf{e}_{k})_{j}=\delta_{kj}$ ($\delta$ here is the
Kronecker’s delta). Let $\mathbf{f}_{k}=\mathbf{F}\mathbf{e}_{k}$ denote the
$k$-th column of $\mathbf{F}$ and let
$\mathbf{p}_{k}=\mathbf{P}_{TS}\mathbf{e}_{k}$. Then, solving Equation (14) is
equivalent to solving the discrete Laplace equation

$\Delta\mathbf{f}_{k}=\mathbf{p}_{k}$ (15)

for all $k\in S$. The standard basis vectors $\mathbf{e}_{k}$ provide exactly
the _Dirichlet boundary conditions_ on the set $S$ (the set $S$ can be assumed
to be a boundary of $T$).

It is also easy to see that Equation (12) can be written as

$\mathbf{H}\Delta=\mathbf{P}_{ST}.$ (16)

Hence, the solution to (16) is obtained by solving the discrete Laplace
equation in terms of the discrete Laplace operator of the transpose of
$\mathbf{P}$.

The _Green’s function_ is defined to be the inverse of the Laplacian. In our
case the inverse of $\Delta$ is exactly the matrix
$\mathbf{G}=(\mathbb{I}-\mathbf{P}_{TT})^{-1}$ and hence the existence of
solutions to Equations (12) and (7) is equivalent to existence of the Green’s
functions to the corresponding Laplacian. In the absorbing Markov chain theory
(Kemeny and Snell, 1976), the matrix $\mathbf{G}$ is known as the _Fundamental
matrix_ of the corresponding absorbing Markov chain. The entry $G_{ij}$
represents the mean number of times the random walk reaches vertex $j\in T$
having started in state $i\in T$.

We now present some elementary sufficient conditions for existence of the
Green’s functions of the discrete Laplacians of the graphs. The full proofs
are given in Appendix A. For the development of the discrete Green’s functions
(for undirected graphs) in terms of the eigenvalues and eigenfunctions of the
Laplacian, we refer the reader to the paper by Chung and Yau (2000).

######  Proposition 2.1.

Suppose that $\Gamma$ is a weighted directed graph such that for every $p\in
T$ there exists $s\in S$ such that there exists a directed path from $p$ to
$s$. Then, the matrix $\mathbb{I}-\mathbf{P}_{TT}$ is invertible and

$(\mathbb{I}-\mathbf{P}_{TT})^{-1}=\sum_{k=0}^{\infty}(\mathbf{P}_{TT})^{k}.$
(17)

Proposition 2.1 thus guarantees existence of the Green’s functions if every
transient vertex can be connected to a source or sink via a directed path. If
the underlying graph is undirected, this condition can be rephrased as
follows: every connected component of $V$ contains at least one vertex from
$S$.

In the context of information diffusion, the connectivity condition implies
that all information entering the transient set at any specific time must
eventually leave it, either by absorption into $S$ when $S$ is a set sinks, or
by dissipation when $S$ represents the set of sources. We will further discuss
the concept of dissipation in 2.3.

Assuming the Green’s function exists, the entries of the matrices $\mathbf{F}$
and $\mathbf{H}$ can be interpreted in several different ways. Fundamentally,
both $F_{ij}$ and $H_{ij}$ represent the total expected number of times the
vertex $j$ is visited by the information originating at the vertex $i$ while
avoiding all members of the boundary set $S$ (the proofs are given in Appendix
B.1). It is also clear, by Equation (17), that $\mathbf{F}$ and $\mathbf{H}$
are both non-negative matrices and that
$\mathbf{F}=\lim_{t\to\infty}\mathbf{F}(t)$ and
$\mathbf{H}=\lim_{t\to\infty}\mathbf{H}(t)$. In addition, the rows of
$\mathbf{F}$ all sum to $1$ (Lemma B.3 in Appendix B.2) and thus $F_{ij}$ is
the overall probability an information originating from transient vertex $i$
is absorbed at the sink $j$ while avoiding all other sinks.

If we assume that a random walk deposits a fixed amount of information content
each time it visits a node, we can interpret $H_{ij}$ is the overall amount of
information content originating from the source $i$ deposited at the transient
vertex $j$. If $\Gamma$ is an undirected graph with symmetric weight matrix
$\mathbf{W}$ and $S$ contains a single source, the value of $H_{ij}$ is
directly proportional to the degree of the transient vertex $j$ (Appendix
B.2). Hence, in this case, the total average number of times of visits for
each transient node is proportional to its degree. This is no longer true if
$\mathbf{W}$ is not symmetric.

Furthermore, we can interpret $F_{ij}$ as the sum of probabilities of paths
originating at the vertex $i\in T$ and terminating at the vertex $j\in S$ that
avoid all other nodes in the set $S$, and $H_{ij}$ as the sum of probabilities
of paths originating at the vertex $i\in S$ and terminating at the vertex
$j\in T$, also avoiding all other nodes in the set $S$. Each such path has a
finite but unbounded length. However, unlike $F_{ij}$, $H_{ij}$ does not
represent a probability because the events of the information being located at
$j$ at the times $t$ and $t^{\prime}$ are not mutually exclusive (a random
walk can be at $j$ at time $t$ and revisit it at time $t^{\prime}$). For
$F_{ij}$, the absorbing events at different times are mutually exclusive.

###  2.3 Information dissipation

It was mentioned previously that the requirement that every transient node is
connected to a node in the set $S$ is effectively equivalent to the property
that all information content entering the transient set leaves it at the nodes
in $S$. In the present section we extend our model to allow the information to
dissipate not only at those nodes but also at the transient nodes.

Let $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ be vectors of length $n$
such that for all $i\in V$, $\alpha_{i}>0$ and $\beta_{i}>0$. We form the
matrix $\mathbf{\tilde{P}}$ with entries

$\tilde{P}_{ij}=\alpha_{i}\beta_{j}P_{ij},$ (18)

and use the new matrix to compute the matrices $\mathbf{\tilde{F}}$ and
$\mathbf{\tilde{H}}$ by replacing the matrix $\mathbf{P}$ in the previous
section with $\tilde{P}$ so that.

$\mathbf{\tilde{F}}=\mathbf{\tilde{G}}\mathbf{\tilde{P}}_{TS}.$ (19)

and

$\mathbf{\tilde{H}}=\mathbf{\tilde{P}}_{ST}\mathbf{\tilde{G}}.$ (20)

where $\mathbf{\tilde{G}}=(\mathbb{I}-\mathbf{\tilde{P}}_{TT})^{-1}$, provided
$\mathbb{I}-\mathbf{\tilde{P}}_{TT}$ is invertible.

The entry $\alpha_{i}$ gives the proportion of the signal leaving the vertex
$i$ that is retained (we call the value of $1-\alpha_{i}$ the _outgoing
dissipation coefficient_ of the node $i$) while the entry $\beta_{j}$ gives
the proportion of the signal entering the vertex $j$ that is retained (the
value $1-\beta_{j}$ is called the _incoming dissipation coefficient_ of the
node $j$). The case where $\alpha_{i}=\beta_{i}=1$ for all $i\in V$ gives back
the original matrix $\mathbf{P}$. Note that our definition allows entries of
$\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ that are greater than $1$,
corresponding to negative dissipation coefficients. Such coefficients lead to
amplification of the signal. However, in order for the Green’s function
$\mathbf{\tilde{G}}$ to exist, any amplification should be balanced by
dissipation.

We now establish a sufficient condition for existence of $\mathbf{\tilde{G}}$.
The proof, as well as a discussion of its generalization, is given in Appendix
A.1.

######  Proposition 2.2.

Let $\alpha_{*}=\max\\{\alpha_{i}:i\in V\\}$ and
$\beta_{*}=\max\\{\beta_{i}:i\in V\\}$ and suppose $\alpha_{*}\beta_{*}<1$.
Then, the matrix $\mathbb{I}-\mathbf{\tilde{P}}_{TT}$ is invertible and

$(\mathbb{I}-\mathbf{\tilde{P}}_{TT})^{-1}=\sum_{k=0}^{\infty}(\mathbf{\tilde{P}}_{TT})^{k}.$
(21)

Proposition 2.2 makes no assumptions on the connectivity of the graph: the
equilibrium solutions exist regardless of the graph topology. The reason for
the removal of the connectivity conditions is that a unit of information
originating anywhere in the network has a nonzero probability of being
dissipated at each time step and therefore will disappear in the long term,
with a portion possibly reaching a sink in the absorbing model. The vectors of
coefficients $\alpha$ and $\beta$ provide us with the ability to consider
different rates of dissipation at different vertices. We demonstrate the
utility of the extended model in examples involving protein-protein
interaction networks (Section 4), where we use vertex specific dissipation to
construct ‘evaporating nodes’ that dissipate most of the information coming in
but allow unrestricted outward flow.

A possible further generalization of this model is for the entries of the
vectors $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ to be functions of the
state variable of the dynamical system instead of constants. The dynamical
system in this case would become non-linear, allowing us to model
amplification or dissipation of the information depending on the time specific
state of the system.

###  2.4 Potentials

Our models so far, including the dissipation modifications described above,
model ‘free diffusion’ of information through the network: the likelihood for
the signal to move from vertex $i$ to vertex $j$ is proportional to the
relative weight of the edge $(i,j)$ among all edges emanating from $i$
(dissipation only affects the total amount transmitted). In order to direct
the flow of information towards or away from selected nodes, we adjust the
weights of edges of our network graph $\Gamma$ using _potentials_ , real-
valued monotone functions defined on the nodes that depend on the distances
from selected points.

Let $\rho$ denote the path-metric on the weighted directed connected graph
$\Gamma=(V,E,w)$, where for all $i,j\in V$, $\rho(i,j)$ denotes the sum of the
reciprocals of the weights of the edges forming the shortest directed path
from $i$ to $j$. Suppose $R$ is a subset of $T$ such that for each $k\in R$
there exists a monotone potential function
$\theta_{k}:\mathbb{R}\to\mathbb{R}$. For each vertex $j\in V$ define the
_total potential_ at $j$, denoted $\Theta(j)$ by

$\Theta(j)=\sum_{k\in R}\theta_{k}(\rho(j,k)).$ (22)

Let $\hat{\Gamma}$ denote the new weighted directed graph $(V,E,\hat{w})$
where

$\hat{W}_{ij}=W_{ij}\exp\left(-\Theta(j)\right).$ (23)

The form of Equation (23) ensures that the signal preferentially diffuses from
each vertex towards the vertices adjacent to it that have lower potential
relative to other adjacent vertices.

A vertex $i\in V$ is called a _destination_ if $\Theta$ has a minimum at $i$.
There can be multiple destinations in a network. The natural candidates for
destinations are the members of the set $S$ since all information entering
them does not leave them. Some transient states, with the weights of their
outgoing edges adjusted to partially accumulate the signal, are also good
candidates for destinations.

Let $K$ be a subset of $T$ and let $0\leq\gamma\leq 1$. From the already
modified graph $\hat{\Gamma}$, we form the graph $\Gamma^{\prime}$ represented
by the weight matrix $\mathbf{W}^{\prime}$ where

$W^{\prime}_{ij}=\begin{cases}\hat{W}_{ij}&\text{if $i\not\in K$,}\\\
\gamma\hat{W}_{ij}&\text{if $i\in K$ and $i\neq j$,}\\\
\hat{W}_{ij}+(1-\gamma)\sum_{k\neq i}\hat{W}_{ik}&\text{if $j\in K$ and
$i=j$}.\end{cases}$ (24)

The effect of this modification is to turn each vertex $i\in K$, called a
_pseudosink_ , into a partial sink: some proportion of the weights of edges
emanating out of $i$ is transferred to the edge pointing back to $i$. The
parameter $\gamma$, representing the proportion of information allowed to
leave each pseudosink while the remainder is accumulated, is called the
_pseudosink leakage coefficient_. The value $\gamma=1$ implies no change in
edge weights.

The value $\gamma=0$ is a special case because no directed path exists between
pseudosinks and source nodes in the resulting graph $\Gamma^{\prime}$ and
Proposition 2.1 does not apply. In this case, there are two possibilities
leading to the existence of the Green’s function: either set the outgoing
dissipation coefficient of the pseudosinks to something less than $1$, or
treat the pseudosinks as parts of the boundary set $S$, as a ‘non-emitting
source’ defined in 3.2 below.

Note that, while dissipation is applied to the transition matrix $\mathbf{P}$,
potentials and pseudosinks are applied to the weight matrix $\mathbf{W}$ prior
to normalization. Since applications of potentials and pseudosinks do not
commute, potentials are applied before pseudosinks, although pseudosinks can
be potential centers (members of the set $R$).

##  3 Theoretical Methods for Analysis

In the previous section we introduced the basic concepts related to our models
of diffusion of information through networks as well as some modifications to
the underlying graph and the transition matrix that lead to biologically
realistic models. After all modifications are applied, we obtain the matrices
$\mathbf{\tilde{F}}$ and $\mathbf{\tilde{H}}$, the Green’s functions arising
where $S$ represents sinks and sources, respectively. Here we turn to the
practical interpretation of these results, which depend on the boundary
conditions imposed on the vertices in $S$.

###  3.1 Absorbing model

In the case where $S$ represents sinks of information (the _absorbing model_),
the entries of the matrix $\mathbf{\tilde{F}}$ have a clear probabilistic
interpretation: $\tilde{F}_{ij}$ is the probability that information starting
at transient vertex $i$ reaches the sink $j$ while avoiding all other sinks,
taking into account the dissipation as well as the new weights induced by the
potentials. Generally, each sink $j$ exerts a ‘region of influence’, including
the transient points with large $\tilde{F}_{ij}$. Depending on the
distributions of sinks within the network, some transient node may have a
$\tilde{F}_{ij}$ small for all $j$: information emerging from these points is
more likely to dissipate than to reach any of the sinks.

If $S^{\prime}\subset S$ is a selection of sink nodes, then $\sum_{j\in
S^{\prime}}F_{ij}$ gives the total probability of information reaching the set
$S^{\prime}$ from the vertex $i$, avoiding all other nodes in $S$. In this
context, we call the nodes in $S^{\prime}$ _explicit sinks_ (since we
investigate the probabilities of reaching them) and the remaining nodes in $S$
_implicit sinks_ , the points that serve as sinks of information but are not
considered. Furthermore, if the sinks are treated as general boundary points,
with boundary values not restricted to $0$ and $1$, the entries of
$\mathbf{\tilde{F}}$ can be interpreted as temperatures (Zhang _et al._ ,
2007).

###  3.2 Emitting model

Where $S$ represents sources (the _emitting model_), the entries of
$\mathbf{\tilde{H}}$ can be interpreted as visiting times or as information
contents: $\tilde{H}_{ij}$ is the total information content emitted from the
source $i$ deposited at the transient vertex $j$. Information is dissipated at
all sources and the value of $\tilde{H}_{ij}$ is dependent on transient
dissipation coefficients $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ and
the potentials. For biological applications, we will consider the case where
at least one pseudosink is present in addition to one or several sources, with
the potential directing the flow towards the pseudosinks. The distribution of
entries of the $i$-th row of $\mathbf{\tilde{H}}$ will then describe the
_information transduction module_ (ITM) involved in transfer of information
from $i$ to the pseudosinks, with the nodes with largest entries being most
significant.

Let $\boldsymbol{\xi}$ denote the vector of length $\left|S\right|$ such that
for all $i\in S$, $\xi_{i}\geq 0$. We call $\xi_{i}$ the _source strength_ of
the source $i$, representing the amount of information emitted from $i$ at
each time step. In this context, we call $i\in S$ an _emitting source_ if
$\xi_{i}>0$ and a _non-emitting source_ if $\xi_{i}=0$. Non-emitting sources
are essentially information ‘black holes’, dissipating any information coming
in and not emitting any.

####  3.2.1 Total content

For any $i\in S$, let $\boldsymbol{\epsilon}_{i}$ denote the standard $i$-th
row basis vector of length $n-m$, where
$(\boldsymbol{\epsilon}_{i})_{j}=\delta_{ij}$. For $x>0$ define the vector
$\boldsymbol{\phi}_{i}$ by

$\boldsymbol{\phi}_{i}=\xi_{i}\boldsymbol{\epsilon}\mathbf{\tilde{H}},$ (25)

that is, $\boldsymbol{\phi}_{i}$ denotes the $i$-th row of
$\mathbf{\tilde{H}}$ multiplied by $\xi_{i}$. Its entries give the amount of
information content originating from the source $i$ of strength $\xi_{i}$
deposited at transient vertices. The value of
$\left\|\boldsymbol{\phi}_{i}\right\|_{1}$ is then the total amount of content
originating at source $i$ deposited at the transient states. In our examples
in the following sections we choose the source strengths $\boldsymbol{\xi}$ so
that $\left\|\boldsymbol{\phi}_{i}\right\|_{1}$ is the same for all $i\in S$
(we call the resulting vectors $\phi_{i}$ normalized content vectors). The
_joint information content_ vector, denoted $\boldsymbol{\tau}$, is defined by

$\boldsymbol{\tau}=\sum_{i\in S}\boldsymbol{\phi}_{i}.$ (26)

The vector $\boldsymbol{\tau}$ implicitly depends on the matrix
$\mathbf{\tilde{H}}$ and the source strength vector $\boldsymbol{\xi}$: we
have $\boldsymbol{\tau}=\boldsymbol{\xi}\mathbf{\tilde{H}}$.

####  3.2.2 Participation ratio

Let $\mathbf{x}\in\mathbb{R}^{n}$ be any vector and recall that for any $0\leq
p<\infty$, the $\ell_{p}$-norm of $\mathbf{x}$, denoted
$\left\|\mathbf{x}\right\|_{p}$, is given by
$\left\|\mathbf{x}\right\|_{p}=\left(\sum_{k}\left|x_{k}\right|^{p}\right)^{1/p}$.
Define the _participation ratio_ of $\mathbf{x}$, denoted $\pi(\mathbf{x})$ by

$\pi(\mathbf{x})=\frac{\left\|\mathbf{x}\right\|_{1}^{2}}{\left\|\mathbf{x}\right\|_{2}^{2}}=\frac{\left(\sum_{k}\left|x_{k}\right|\right)^{2}}{\sum_{k}x_{k}^{2}}.$
(27)

Participation ratio is well known under a slightly different definition in the
physics literature (Thouless, 1974). It gives the number of components of
$\mathbf{x}$ whose magnitude is ‘significant’. Clearly, $\pi$ is independent
of the scale of $\mathbf{x}$: we have for any $\lambda>0$,
$\pi(\lambda\mathbf{x})=\pi(\mathbf{x})$. We illustrate the usage by examples.

######  Example 3.1.

Let $\mathbf{x}=[1,1,1,1,1]$. Then, $\pi(\mathbf{x})=\frac{5^{2}}{5}=5$. All
components are equally significant and this is reflected in the participation
ratio.

######  Example 3.2.

Now consider $\mathbf{x}=[1,1,0,0,0]$. We have,
$\pi(\mathbf{x})=\frac{2^{2}}{2}=2$. Only the first two components are non-
zero and are of equal magnitude.

######  Example 3.3.

Finally, let
$\mathbf{x}=\left[1,\frac{1}{2},\frac{1}{4},\frac{1}{8},\frac{1}{16}\right]$.
We obtain $\pi(\mathbf{x})\approx 2.8181$. Here all five components are non-
zero but their magnitudes differ significantly. The participation ratio here
implies that the first two components and to a large extent the third are
significant while the remaining two are much smaller.

In our biological examples, we use $\pi(\boldsymbol{\tau})$ to choose the
number of the transient vertices with largest total mass to display as a
‘significant’ subgraph, together with all sources and pseudosinks.

####  3.2.3 Interference

Given the vector of source strengths $\boldsymbol{\xi}$, the entry of
$\tau_{j}$ can be interpreted as providing the total amount of information
deposited at the vertex $j$. It is also possible to investigate the
interaction of the signals from different sources using the concept of
destructive interference.

For any vector $\mathbf{x}\in\mathbb{R}^{n}$, let $\mu$ denote an
_interference function_ such that
$0\leq\mu(\mathbf{x})\leq\left\|\mathbf{x}\right\|_{1}$. When applied to a
vector containing information content from different sources, interference
function is interpreted as removing some of the information present due to the
interaction of the various information types and returning the remaining
information content. Interference functions can take various forms depending
on the nature of the types of information in each application.

######  Example 3.4.

Suppose $\mathbf{x}$ consists of two components representing information types
that are assumed to completely cancel out each other. In this case, the
interference function takes the form
$\mu(\mathbf{x})=\left|x_{1}-x_{2}\right|$.

######  Example 3.5.

When $\mathbf{x}$ has more than two components, there are may possible ways to
generalize the above example. We distinguish two general modes of
interference: exclusive and partial. Exclusive interference mode represents
the case where simultaneous presence of all types of information is necessary
for destructive interference. For example, if each information type carries
the same weight, the interference function is:

$\mu(\mathbf{x})=\sum_{k}\left(x_{k}-\nu\right),$ (28)

where $\displaystyle\nu=\min_{k}x_{k}$.

######  Example 3.6.

We call the partial interference the case where presence of all types of
information is not necessary. It can be modeled in many ways depending on the
desired interpretation. For example, if there are three sources, we can use
complex numbers to set $\mu$ so that

$\mu(\mathbf{x})=\left|\sum_{k=1}^{3}x_{k}\exp\left(\frac{\iota
k\pi}{3}\right)\right|,$ (29)

where $\iota$ denotes the imaginary unit. In this case, some content is lost
when any two types of signal are present but all three must be present for
complete annihilation.

Given the interference function $\mu$, define the _interference strength
function_ $\psi:\mathbb{R}^{n}\to\mathbb{R}\cup\\{\infty\\}$ by

$\psi(\mathbf{x})=\begin{cases}\left\|\mathbf{x}\right\|_{1}\log\left(\frac{\left\|\mathbf{x}\right\|_{1}}{\mu(\mathbf{x})}\right)&\text{if
$\left\|\mathbf{x}\right\|_{1}>0$,}\\\ 0&\text{if
$\left\|\mathbf{x}\right\|_{1}=0$.}\end{cases}$ (30)

By the definition of $\mu$

Since $0\leq\mu(\mathbf{x})\leq\left\|\mathbf{x}\right\|_{1}$, it follows that
$\psi$ takes non-negative values (including $+\infty$). The value of $\psi$ is
infinite if $\mu(\mathbf{x})=0$ (perfect interference) and finite otherwise.
For an $m\times n$ matrix $\mathbf{X}$ define the vector
$\boldsymbol{\sigma}(\mathbf{X})$ of length $n$ having the components

$\sigma_{i}(\mathbf{X})=\psi(\mathbf{X}\mathbf{e}_{i})$ (31)

(recall that $\mathbf{e}_{i}$ is the standard column basis vector and hence
$\mathbf{X}\mathbf{e}_{i}$ represents the $i$-th column of $\mathbf{X}$). We
will call $\boldsymbol{\sigma}$ the _interference strength vector_.

For our applications, the entries of the matrix $\mathbf{X}$ above are
interpreted as information contents over some graph: $X_{ij}$ is the the
content of type $i$ at the vertex $j$. For each node $j$, the $\ell_{1}$-norm
in Equation (30) can be interpreted in this context as the total information
content at $j$ and the value of $\mu$ applied to the $j$-th column of
$\mathbf{X}$ as the information content remaining after interference. Hence,
interference strength of each node measures how much information content was
lost by interference, adjusted by the node’s joint information content.

The matrix $\mathbf{\tilde{H}}$ is therefore a natural input to $\psi$ and
$\boldsymbol{\sigma}$, however other derived matrices can be used such as
$\mathbf{\tilde{H}}$ adjusted for source strength by multiplying each row by
its corresponding source strength $\xi_{i}$. Furthermore, rows of $\mathbf{X}$
can come from different $\mathbf{\tilde{H}}$ matrices, using different
potentials or dissipation coefficients, as long as the underlying vertex set
is the same. The general purpose of interference strength is to measure the
amount of interaction or overlap between different ITMs.

##  4 Biological Examples

The theory and methods outlined in previous sections can be applied to any
interaction network. This section will present some examples using biological
networks, more specifically, yeast protein-protein interaction networks. Since
the interaction data obtained using many high-throughput methods is generally
inconsistent (Sprinzak _et al._ , 2003), we use the core yeast dataset from
DIP, version ScereCR20060402, consisting of 2554 proteins and 5952
interactions for all our examples. The core dataset, obtained using the
methods of Deane _et al._ (2002), contains only the most reliable interactions
from the DIP dataset of all yeast protein-protein interactions.

Our examples are restricted to investigation of information transduction
modules related to yeast histone acetyltransferases (HATs). Histones are
nuclear proteins that are major components of eukaryotic chromatin (Wolffe,
1992): eukaryotic DNA is organized as a repeating array of nucleosomes
consisting of 146 bp of DNA wound around a histone octamer consisting of two
of each of histone proteins H2A (Hta1, Hta2 in yeast), H2B (Htb1, Htb2 in
yeast), H3 (Hht1, Hht2 in yeast) and H4 (Hhf1, Hhf2 in yeast). It has been
repeatedly demonstrated that transcription is strongly influenced by the
chromatin structure and DNA-histone interactions in particular. The regions of
DNA that interact with histones are generally unavailable for transcription
and transcriptional activation and deactivation are connected with chromatin
alterations (Wolffe, 2001).

Histone acetyltransferases are enzymes that acetylate histones, leading to
weakening of the nucleosome structure and making the DNA involved accessible
to transcription factors (Struhl, 1998; Workman and Kingston, 1998).
Saccharomyces cerevisiae contains several HATs from two major classes with a
variety of biological functions and substrate specificities (Sterner and
Berger, 2000). The proteins Hat1, Gcn5, Elp3, Spt10 and Hpa2 belong to the
GNAT superfamily (Neuwald and Landsman, 1997), while Esa1, Sas2 and Sas3
belong to the MYST family (Borrow _et al._ , 1996; Smith _et al._ , 1998). The
proteins TAF1 (TATA-binding protein associated factor), a subunit of the TFIID
complex, and Nut1 (Med5), a subunit of the mediator complex (Biddick and
Young, 2005), have also been associated with histone acetyltransferase
activity (Mizzen _et al._ , 1996; Lorch _et al._ , 2000).

Unfortunately, the core dataset does not contain the relevant data for all
known HATs. The HATs Hpa2 and Spt10 are not present in the core while HAT1 has
interactions only with Hat2 and its substrate Hhf2. We chose to primarily
concentrate on HATs Gcn5, Esa1 and Elp3 because they are well researched and
the interaction data is abundant. They are all involved in transcriptional
activation, unlike Sas2, which promotes silencing (Osada _et al._ , 2001).

Gcn5 is the best characterized of all HATs, preferentially acetylating histone
H3 (Sternglanz and Schindelin, 1999). It forms the catalytic subunit of the
ADA and SAGA transcriptional activation complexes (Grant _et al._ , 1997). In
addition to Gcn5, the SAGA complex also contains the proteins Tra1, TAF5,
TAF6, TAF9, TAF10, TAF12, Hfi1 (Ada1), Ada2, Ngg1 (Ada3), Spt3, Spt7, Spt8 and
Spt20 (Ada5) (Timmers and Tora, 2005). The ADA complex contains a subset of
proteins from the SAGA complex, namely Gcn5, Hfi1, Ada2, Ngg1 and Spt20, plus
the adaptor protein Ahc1 (Eberharter _et al._ , 1999). The TAF proteins in
SAGA also belong to the TFIID complex, which overall consists of 15 subunits
including a TATA-binding protein and 14 TAFs (Sanders and Weil, 2000).

Esa1 is the catalytic subunit of the NuA4 histone acetyltransferase complex
essential for growth in yeast (Smith _et al._ , 1998; Allard _et al._ , 1999)
that catalyses acetlyaltion of the histone H4. It has been established that
the NuA4 complex, containing, in addition to Esa1, the proteins Tra1, Epl1
Yng2, Eaf1, Eaf2, Eaf3, Eaf5, Eaf6, Act1, Arp4 and Yaf9, is recruited by a
variety of transcriptional complexes as a transcriptional coactivator and is
involved in DNA repair (Doyon and Cote, 2004).

Elp3 is a part of the six component elongator complex , which is associated
with RNA polymerase II during transcript elongation (Wittschieben _et al._ ,
1999). The elongator complex also includes the proteins Iki3 (Elp1), Elp2–4,
Iki1 (Elp5) and Elp6 (Krogan and Greenblatt, 2001).

This section contains four examples of the application of our models, depicted
in Figures 1–5. Subsection 4.2 describes possible complexes associated with
the HATs Gcn5, Esa1 and Elp3, taken individually and in competition, that can
be inferred from the protein-protein interaction network using the absorbing
model. Subsection 4.3 investigates possible physical interaction interfaces
between the MADS box protein Mcm1 (Shore and Sharrocks, 1995) and the HATs
Esa1 and Gcn5. In this case, the emitting model is employed to discover the
pathways through which Mcm1 can recruit the above HATs and whether they are
recruited through the same interface. Before presenting our results we
describe the model parameters and computational techniques used.

###  4.1 Parameters and computation

####  4.1.1 Dissipation

For all our examples, we set $\alpha_{i}=1$ for every node $i$ in our
interaction network so that the outgoing flow from any node is not dissipated.
Modeling the incoming dissipation the coefficients $\beta_{i}$ can take two
values: one for ‘ordinary’ and one for _evaporating_ vertices. In our examples
that use the absorbing model (4.2), $\beta_{i}$ is set to $0.70$ for ordinary
nodes and $0.01$ for evaporating nodes while the examples using the emitting
model (4.3) set $0.87$ for ordinary nodes and $0.01$ for evaporating nodes.
The evaporating nodes consisted of cytoskeleton proteins Act1, Myo1, Myo2,
Myo3, Myo4, Myo5, Smy1, Smy2, Sla1, Arc40, Arp2, Rvs167, Tpm1, Tpm2, Aip1 and
Las17 and histones (Hta1, Hta2, Htb1, Htb2, Hht1, Hhf2, Htz1, Hho1).

The coefficients for the ordinary nodes were chosen using the following
reasoning. For the emitting model we considered the dissipation rate that
would allow the random walk emitted from the source to reach an ‘average’ node
along the shortest path to it with the probability slightly less than $0.5$,
say $0.49$. We found that the average length of the shortest path between two
points in the yeast core dataset is $5.23$ and hence our coefficient is
$0.49^{(1/5.23)}=0.872$, which is rounded to $0.87$. A different coefficient
was needed for the absorbing examples because we were interested in only the
immediate complexes containing our selected HATs: the coefficient
$\beta_{i}=0.87$ would lead to most of the members of the RNA polymerase II
holoenzyme to be retrieved as members of the resulting ITM. We chose to
consider the shortest paths of length $2$, rather than of the average length
$5.23$. Using the same calculation as above, we obtain $0.49^{(1/2)}=0.7$.

The reason for having evaporating nodes with larger dissipation rate is that
both the cytoskeleton proteins and the histones form extended structures in
the cell and the nucleus, respectively. In our physical interaction network,
we assume that information can flow from one protein to another through an
intermediate node if all three nodes are brought close together in space and
time. Information is not likely to flow through proteins that are parts of
extended structures because proteins with completely different biological
function may bind them at different locations and at different times.
Therefore, allowing significant information flow through such nodes would
yield biologically implausible results.

However, depending on the exact context of the investigation, such nodes may
have an important role to play and removing them completely from the
interaction networks or assigning them to the boundary set $S$ would not be
appropriate. Hence, we set a very high incoming dissipation rate at
evaporating nodes while allowing the information to originate from them. In
terms of our models, this approach means that the evaporating nodes will have
very small visiting times in the emitting models and hence will not be
components of any ITM. On the other hand, depending on the exact network
topology, they may be part of ITMs obtained by the emitting model. Note that
other proteins that bind their interacting partners in a non space and time
specific manner can be chosen as additional evaporating nodes; we chose
histones and cytoskeleton proteins due to their direct relevance to our
selected examples.

####  4.1.2 Potentials

All our examples use attracting potentials centered at each pseudosink or
sink. The potential function, heuristic in nature, is the same in every
example has the the form

$\theta_{k}(x)=\begin{cases}a_{1}x&\text{if $0<x\leq b$,}\\\
a_{1}x+a_{2}(x-b)^{2}&\text{if $x>b$,}\end{cases}$ (32)

where $a_{1}=0.8181$, $a_{2}=0.05$, $b=2$ and $k$ is any pseudosink or a sink.
The potential function shown above is long-range, affecting the whole graph,
with a linear portion for short ranges $0\leq x\leq 2$ and quadratic for
distances larger than $2$. We do not expect to see qualitative changes in the
results if the form of the potential function is modified as long as it has
the effect of attracting information towards the destination.

The sources (in the case of emitting models) and evaporating points were
excluded from the graph prior to calculating distances (their distances from
the centers were set to an arbitrary large number) in order to exclude the
paths passing through them from consideration. The reason for excluding the
paths passing through sources was that, by construction, the information never
enters a source from a transient vertex, while the evaporating points were
excluded because most of the signal entering them is dissipated.

####  4.1.3 Numerical implementation

The code for computation of the results was implemented in the Python
programming language, using the NumPy and SciPy packages (Jones _et al._ ,
2001–). In particular, the computation of the matrices $\mathbf{\tilde{F}}$
and $\mathbf{\tilde{G}}$ (Equations (19–20)) was performed by the embedded
FORTRAN code from the UMFPACK (Davis, 2004) solver of sparse systems of linear
equations, using the Automatically Tuned Linear Algebra Software (ATLAS)
(Whaley and Petitet, 2005) implementation of Basic Linear Algebra Subprograms
(BLAS). The graphical representations of the subgraphs of interest were
produced by the neato program from the Graphviz graph visualization suite
(Gansner and North, 2000).

###  4.2 HAT complexes: absorbing examples

(a) |
---|---
|
(b) |
|
(c) |
|
Figure 1: ITMs obtained by running the absorbing model with Esa1(a), Gcn5(b)
and Elp3(c) as a sink. The shades of grey at the nodes represent the
probability of the information originating at the corresponding protein being
absorbed at the sink, the darker nodes indicating higher probability.

Figure 1 shows the three subgraphs of the yeast core interaction graph
consisting of the top scoring nodes according to the absorbing model with
Esa1, Gcn5 and Elp3 as single sinks, respectively. The information orginating
at the proteins shown has more than $0.07$ probability of being absorbed by
the sink (under the influence of the potential centered at the sink) as
opposed to being dissipated. Hence, the subgraphs show the proteins that are
likely to be in the same complex with the HATs chosen as sinks.

Figure 1(a), with Esa1 as the sink, shows all the proteins from the NuA4
complex that are available in the core dataset as highly significant. Some of
the proteins from ADA and SAGA complexes can also be seen because Tra1 belongs
to these complexes as well as to NuA4. The four types of histones forming the
histone octamer can also be seen interacting with Arp4. The proteins Vps51–54
on the right of Figure 1(a) belong to the Vps Fifty-three thethering (VFT)
complex, involved in vesicle assembly (Reggiori _et al._ , 2003). The proteins
Tlg1 and Ypt6 are interacting partners of the VFT complex (Reggiori _et al._ ,
2003). The relation between VFT and NuA4 is not established as these two
complexes are localized in different cellular compartments: NuA4 in the
nucleus and VFT in golgi-vacuole transport vesicles. The relationship observed
in Figure 1(a) results exclusively from the Yng2–Vps51 interaction, which was
orginally observed in a yeast-two-hybrid screen by Ito _et al._ (2000, 2001).
Based on the above information, it appears that VFT and NuA4 complexes do not
interact _in vivo_. Note that the histones as well as actin, although selected
as evaporating points, can be seen in the figure because the outgoing flow
from evaporating nodes is allowed.

In a similar fashion, Figure 1(b), with Gcn5 as the sink, shows the members of
SAGA, ADA and TFIID transcriptional activator complexes as well as many other
transcription factors, mostly members of subcomplexes of the RNA polymerase II
holoenzyme. Also worth mentioning is Cti6, which bridges the Cyc8-Tup1
corepressor and the SAGA coactivator to overcome repression of the GAL1 gene
(Papamichos-Chronakis _et al._ , 2002). The Cyc8 protein is also shown while
Tup1 is not, most likely because it is involved in many other interactions
away from Gcn5, bringing down its relative significance. Figure 1(c), with
Elp3 as the sink, clearly outlines the elongator complex, as well as some
members of the core RNA polymerase II complex (Rbp2–5, Rbp7, Rpc10, Rpo26)
(Myer and Young, 1998).

Figure 2 shows the top scoring nodes according to the absorbing model with
Esa1, Gcn5 and Elp3 as simultaneous sinks with attracting potentials. In this
case, the information originating at the depicted nodes has more than $0.05$
total probability of being absorbed by any of the sinks as opposed to being
dissipated.

Fewer nodes can be seen in this figure as compared to Figure 1 because the
three attracting potentials are now involved that may cancel each other out.
It can be seen that the elongator complex centered around Elp3 is not
connected to the subgraph around Esa1 and Gcn5. Although all of the NuA4,
SAGA, ADA and elongator complexes belong to the RNA polymerase II holoenzyme,
they do so at different times. The NuA4, ADA and SAGA complexes have a role in
initiation of transcription while the elongator complex is involved in
transcript elongation (Martinez, 2002). The green (mixture of cyan and yellow)
color of Tra1 is indicative of the fact that it is a subunit of both
Esa1-containing NuA4 complex and the Gcn5-containing SAGA complex.

Figure 2: ITM obtained by running the absorbing model with Esa1, Gcn5 and Elp3
as simultaneous sinks. The strength of each of cyan, yellow and magenta color
component of the node shows the square root of the probability of absorption
at Esa1, Gcn5 and Elp3, respectively.

###  4.3 Transcription factor interaction interfaces; emitting examples

Mcm1 is a yeast transcription factor essential for cell viability. It controls
many cellular functions including cell cycle transition (Althoefer _et al._ ,
1995), mating (Mead _et al._ , 2002) and arginine metabolism (Messenguy and
Dubois, 1993), through interactions with different cofactors. It has been
determined that Mcm1 acts both as an activator and a repressor of
transcription (Bruhn _et al._ , 1992; Messenguy and Dubois, 1993) and here we
explore the possible ways it can interact with the NuA4 and SAGA HAT
complexes.

(a) |
---|---
|
(b) |
|
Figure 3: ITMs resulting from the emitting model with Mcm1 as a source and
Esa1 as a pseudosink using the original yeast core dataset (a) and the
modified dataset additionally including the edges Tra1–Gal4 and Tra1–Gcn4 (b).
The proteins containing the largest amounts of deposited information are
shown, with the information content indicated by shading (darkest nodes
contain the most information).

Figure 3(a) shows the subgraph consisting of the $22$ proteins with the
largest deposited information content obtained by running our emitting model
with Mcm1 as a source and Esa1 as a pseudosink. The number of proteins to
display ($20$ plus the source and the pseudosink) was chosen because the
participation ratio for the information content vector (excluding the source
and the pseudosink) was $20.33$.

The ITM shown in Figure 3(a) gives the likely pathways of physical interaction
from Mcm1 to Esa1, according to the yeast core interaction dataset. It can be
immediately observed that Esa1 is reached solely through Tra1, which is known
to be the general interaction domain of both NuA4 and SAGA HAT complexes
(Allard _et al._ , 1999; Grant _et al._ , 1998). Directly associated with Mcm1
are the proteins Arg80–Arg82, belonging to the ArgR complex involved in
regulation of arginine metabolism (Dubois and Messenguy, 1991). The majority
of the ITM is dominated by the members of the SRB mediator subcomplex of the
RNA polymerase II holoenzyme (Srb2, Srb4, Srb7) (Biddick and Young, 2005) and
the TFIID, SAGA and ADA complexes. Also prominent are transcriptional
activators Gal4 and Gcn4 (Hinnebusch, 2005; Traven _et al._ , 2006).

The subgraph image suggests two possible interaction pathways: the main (based
on the intensities of deposited information) through Srb4 and members of
SAGA/ADA complex and the alternative through Ume6–TAF10–Spt7. Ume6 is a DNA
binding protein that acts as a transcriptional repressor by recruiting histone
deacetylases, which have the catalytic activity opposite to the HATs (Kassir
_et al._ , 2003). While simultaneous existence of activating and repressing
pathways is biologically plausible, we do not anticipate both pathways to be
in action at the same time. On the other hand, interaction of Mcm1 with the
NuA4 through any of the above pathways _in vivo_ is doubtful because both
pathways lead through the interacting partners of Tra1 in the SAGA complex
that are not associated with it in the NuA4 complex (Doyon and Cote, 2004;
Timmers and Tora, 2005). Note that the direct physical interaction of the
ArgR/Mcm1 complex and the SAGA complex was hypothesized by Ricci _et al._
(2002) in relation to regulation of arginine metabolism.

Nevertheless, it is likely that the yeast core dataset does not contain all
the interactions of Tra1 and that the interactions not in the dataset may
provide us with the plausible explanation. Brown _et al._ (2001) have
indicated that HAT complexes are recruited through Tra1 by Gal4 and Gcn4
transcriptional activators. To investigate if adding the implied edges would
significantly change the resulting ITM we added the Gcn4–Tra1 and Gal4–Tra1
links to the core dataset and rerun the emitting model with all other
parameters unchanged. The resulting ITM, with participation ratio of $21.66$,
is shown in Figure 3(b). We observe few changes: the proteins Ssn3, Srb5, Srb6
and Gal11, belonging to the mediator complex, replaced Cti6 and Srb7, thus
placing more emphasis to the mediator complex.

In this example, our emitting model appears to be quite robust to changes in
the pseudosink leakage parameter $\gamma$. Using the original core dataset, in
addition to the original run with $\gamma=0.3$, we ran our model with
$\gamma=0$, $\gamma=0.5$ and $\gamma=1$, obtaining participation ratios of
$19.43$, $20.34$ and $20.75$ and very little change in constitution of the
ITMs. For example, when $\gamma=1$, the new ITM contains the NuA4 proteins
Arp4 and Yng2 in the place of Cti6 and Srb7. Hence, larger pseudosink leakage
coefficient allows exploration of the nodes surrounding the pseudosinks
without affecting the remainder of the ITM in a major way. Such exploration is
very desirable for protein-protein interaction networks because it reveals
more of the complexes around pseudosinks, thus giving some of the
characteristics of the absorbing model to the emitting model. Note that many
of the interacting partners of the sources are found in the ITM solely due to
proximity of the source.

(a) |
---|---
|
(b) |
|
Figure 4: ITM resulting from the emitting model with Esa1 and Gcn5 as sources
and Mcm1 as a pseudosink: (a) information content, (b) interference strength.

To explore the extend the HATs Esa1 and Gcn5 share their interaction interface
with Mcm1 we set Esa1 and Gcn5 as sources and Mcm1 as a pseudosink
destination. Figure 4 shows the ITM based on the total information content
(participation ratio $24.62$, 28 nodes shown), with the nodes shaded according
to total content and interference strength. The proteins shown as nodes in
Figure 4 have appeared in one of the previous figures, mostly forming parts of
NuA4, SAGA/ADA, TFIID and mediator complexes. The nodes with the largest total
content are Tra1, Ada2, Ngg1 and Srb4 and the latter three are also the nodes
with by far the largest interference strength. This fact does not surprise us
because although Tra1 is a member of both NuA4 and SAGA complexes, information
flowing from Gcn5 to Mcm1 largely avoids it.

The paths used by the information emitted from Esa1 and Gcn5 separately can
best be seen in a color figure (Figure 5(a)) where the information content
from Esa1 and Gcn5 is shown as cyan and yellow, respectively. The nodes
colored strongly cyan contain mostly information from Esa1 while those colored
yellow contain mostly the information from Gcn5. The nodes colored green
contain information from both sources. In this way it can be observed that
members of NuA4 contain the information solely from Esa1, some SAGA proteins
contain the information solely from Gcn5, while Ada2, Ngg1 and Srb4 contain a
significant amount of information from both sources.

(a) |
---|---
|
(b) |
|
Figure 5: Information content of members of the ITM arising from the emitting
model with Esa1 and Gcn5 as sources and Mcm1 as a pseudosink: (a) using the
yeast core dataset; (b) using the modified dataset additionally including the
edges Tra1–Gal4 and Tra1–Gcn4. The strength of the cyan and yellow color
component of the node corresponds to the information content originating from
Esa1 and Gcn5, respectively.

Using additional links based on Brown _et al._ (Figure 5(b)) produces effects
similar to Figure 3(b): the common interface through the mediator complex is
emphasized at the expense of the paths through the SAGA complex. For example,
note the difference in color of Spt7, Gcn4 and Gal4 between Figure 5(a) and
Figure 5(b). The common interface through the mediator complex appears
biologically more plausible than directly through members of the SAGA complex
but we are as yet unable to find direct evidence in the literature confirming
either possibility.

##  5 Discussion and conclusion

The proposed information diffusion models appear to capture some of the
essential features of the yeast protein-protein interaction network in our
examples. Our absorbing model performed well in identifying complexes related
to sinks while the emitting model with pseudosinks is able to illuminate the
possible interaction interfaces between sources and pseudosinks. Application
of the concept of destructive interference in this context provides a way to
assess the degree of overlap of different ITMs.

The salient feature of our models is a novel use of attraction potentials and
dissipation. While the entries of the Green’s function can be interpreted in
graph-theoretic terms as sums of weights of paths from a source to a transient
vertex (for the emitting model) or from a transient vertex to a sink (for the
absorbing model), the potentials, together with the choice of boundary,
provide a unique context for information diffusion in the network. The weights
of the edges and hence the nature of the underlying graphs are changed every
time a different potential is applied, thus bringing forward different aspects
of the network. The potential function used for our examples was heuristic in
nature and we hope that our work would generate interest in developing
theoretical foundations for directed information propagation through networks.

Dissipation coefficients provide a natural and extremely flexible way of
controlling the spread of information content through the network. While
Girvan and Newman (2002) proposed a similar formulation for penalizing longer
paths connecting two nodes in a network, they did so in the context of
hierarchical clustering and using a single dissipation rate. Node specific
dissipation rates are important because they allow construction of
‘evaporating nodes’ and possible integration of additional information to our
model. Having the dissipation rates dependent on the environment of the node
may lead to a more sophisticated model of information transduction.

When modelling physical cellular protein networks, the main limitation of our
approach is that the the publicly available representations of protein-protein
interaction networks contain a limited amount of information. Each interaction
is shown as either occurring or not occurring, without reference to the
dynamics, time-scale, or specificity of binding. Furthermore, the spatial
location of the interactions on the protein molecules is not available, so
that it cannot be determined if a protein known to belong to two separate
complexes, such as Tra1 in our examples, can belong to both at the same time
and therefore transmit information between them. Therefore, our model of
protein cellular networks is only metaphorical at this stage. However, our
diffusion paradigm can be adapted to account for additional information about
proteins, such as their concentrations, cellular compartment localizations,
post-translational modifications or rate constants for binding interactions,
as it becomes available. One way to do that is to associate each protein to a
vector instead of a scalar value and to construct an evolution operator that
reflects the nature of the additional information. In such circusmstances, the
dynamics of information flow could be as revealing as the steady state we use
at this stage.

The quality of the interaction dataset also has a strong influence to the
outcomes of our models. Addition or deletion of edges may make the results
more realistic, as in our emitting examples, but also may completely alter the
ITM produced, if a particular edge provides a shortcut towards the
destination. Hence, in order to obtain the results useful in field of
application, it is imperative to use datasets of interactions that precisely
reflect the network being investigated. In the case of yeast protein-protein
interactions, Collins _et al._ (2007) were recently able to derive a
significantly more reliable collection of interactions, primarily based on two
large-scale studies of protein complexes by tandem affinity purification of
complexes followed by mass spectroscopic identification of individual proteins
(Gavin _et al._ , 2006; Krogan _et al._ , 2006). It is interesting that the
same transcriptional complexes encountered in our examples are prominent in
the unified physical interactome map presented by Collins _et al._ (2007).

The problem of ‘shortcuts’ through the network was also observed by Steffen
_et al._ (2002), who completely eliminated certain nodes in their effort to
model signal transduction pathways using the yeast protein-protein
interactions. Our evaporating nodes, with a very large incoming dissipation
rate, have a similar role with an added advantage that they can be visible as
parts of complexes observed using the absorbing model. The list of evaporating
nodes used by us is not exhaustive and it would be necessary to add further
classes of proteins to it for large-scale investigations of the yeast protein
interactome using our methods.

In this paper, we introduced a flexible mathematical framework for analysis of
interaction networks and indicated its utility by examples. We believe that
the ability to select a particular context for information propagation by
setting various model parameters will be extremely useful for addressing
questions involving interaction networks in biology and many other
disciplines.

##  6 Acknowledgments

We thank Drs John Wootton and David Landsman for encouragement and comments.
This work was supported by the Intramural Research Program of the National
Library of Medicine at National Institutes of Health.

##  Appendix A Existence of Green’s Functions

In this appendix we provide the elementary proofs of the results about
existence of the Green’s functions stated in the main text. As before,
$\Gamma=(V,E,w)$ denotes a weighted directed graph with $N$ vertices, with the
weight matrix $\mathbf{W}$ and transition matrix $\mathbf{P}$. We also have
$T\subset V$ and $S=V\setminus T$.

Recall that for every matrix $\mathbf{M}$, the induced $\ell_{\infty}$ norm of
$\mathbf{M}$, written $\left\|\mathbf{M}\right\|_{\infty}$, is defined by

$\left\|\mathbf{M}\right\|_{\infty}=\sup_{\mathbf{x}\in\mathbb{R}^{n}}\frac{\left\|\mathbf{M}\mathbf{x}\right\|_{\infty}}{\left\|\mathbf{x}\right\|_{\infty}},$
(33)

where $\left\|\mathbf{x}\right\|_{\infty}=\max_{i}\left|x_{i}\right|$. One can
easily show that

$\left\|\mathbf{M}\right\|_{\infty}=\max_{i}\sum_{j}\left|M_{ij}\right|.$ (34)

Also recall that the spectral radius of a square matrix $\mathbf{M}$ is
defined to be the largest absolute value of its eigenvalues. It is well known
that that for every eigenvalue $\lambda$ of $\mathbf{M}$ and any
$k=1,2,\ldots$,

$\left|\lambda\right|\leq\left\|\mathbf{M}^{k}\right\|_{\infty}^{1/k}.$ (35)

######  Lemma A.1.

Let $\mathbf{M}$ be a square matrix with the spectral radius strictly less
than $1$. Then,

  1. (i)

$\mathbf{M}^{k}\to\mathbf{0}$ as $k\to\infty$,

  2. (ii)

The matrix $\mathbb{I}-\mathbf{M}$ is invertible and
$(\mathbb{I}-\mathbf{M})^{-1}=\sum_{k=0}^{\infty}\mathbf{M}^{k}$.

###### Proof.

By the Jordan matrix decomposition, we can write
$\mathbf{M}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}$ for some matrix
$\mathbf{V}$, where $\boldsymbol{\Lambda}$ is a block-diagonal matrix of the
form

$\boldsymbol{\Lambda}=\left[\begin{array}[]{cccc}\mathbf{B}_{1}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&\mathbf{B}_{2}&\cdots&\mathbf{0}\\\ \vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\cdots&\mathbf{B}_{N}\end{array}\right],$

with each of the sub-blocks $\mathbf{B}_{j}$, $1\leq j\leq N$, is of the form
$\mathbf{B}_{j}=\lambda_{j}\mathbb{I}+\mathbf{C}_{j}$ where

$\mathbf{C}_{j}=\left[\begin{array}[]{ccccc}0&1&0&\cdots&0\\\
0&0&1&\cdots&0\\\ \vdots&\vdots&\ddots&\ddots&\vdots\\\ 0&0&0&\cdots&1\\\
0&0&0&\cdots&0\end{array}\right]$

and $\lambda_{1},\ldots\lambda_{N}$ are eigenvalues of $\mathbf{M}$. Hence,
$\mathbf{M}^{k}=\mathbf{V}\boldsymbol{\Lambda}^{k}\mathbf{V}^{-1}$ and

$\boldsymbol{\Lambda}^{k}=\left[\begin{array}[]{cccc}\mathbf{B}_{1}^{k}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&\mathbf{B}_{2}^{k}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\cdots&\mathbf{B}_{N}^{k}\end{array}\right].$

For each eigenvalue $\lambda_{j}$ and each block $\mathbf{B}_{j}$, we can
write

$\mathbf{B}_{j}^{k}=(\lambda_{j}\mathbb{I}+\mathbf{C}_{j})^{k}=\sum_{p=0}^{k}\binom{k}{p}\lambda_{j}^{k-p}\mathbf{C}_{j}^{p}.$

It can easily be shown that for each $j$, $\mathbf{C}_{j}$ is a nilpotent
matrix, that is, if $\mathbf{C}_{j}$ is an $m\times m$ matrix, then
$\mathbf{C}^{m}=\mathbf{0}$. Therefore, for $k\geq m-1$,

$\mathbf{B}_{j}^{k}=\lambda_{j}^{k-m+1}\left(\sum_{p=0}^{m-1}\binom{k}{p}\lambda_{j}^{m-p-1}\mathbf{C}_{j}^{p}\right).$

Observe that the above expression in parenthesis gives an (upper triangular)
matrix whose entries are $m-1$-th degree polynomials in $k$ and hence, that
the whole expression for $\mathbf{B}_{j}^{k}$ is dominated by
$\lambda_{j}^{k-m+1}$. Since, by the spectral radius assumption,
$\left|\lambda_{j}\right|<1$ for each $i$, it follows that for each $j$,
$\mathbf{B}_{j}^{k}\to\mathbf{0}$ as $k\to\infty$ and hence
$\boldsymbol{\Lambda}^{k}\to\mathbf{0}$ as $k\to\infty$ by the block
structure. This proves the first statement.

For the second statement suppose that $\mathbb{I}-\mathbf{M}$ is singular.
Then $\mathbb{I}-\mathbf{M}$ has $0$ as an eigenvalue and hence $\lambda=1$ is
an eigenvalue of $\mathbf{M}$, contradicting our assumption about the spectral
radius of $\mathbf{M}$. Therefore, $\mathbb{I}-\mathbf{M}$ is invertible.
Furthermore, it can easily be obtained using the block diagonal structure of
$\boldsymbol{\Lambda}$ and the ratio test that the sum
$\sum_{k=0}^{\infty}\mathbf{M}^{k}$ converges, Hence,

$(\mathbb{I}-\mathbf{M})\sum_{k=0}^{\infty}\mathbf{M}^{k}=\sum_{k=0}^{\infty}\mathbf{M}^{k}-\sum_{k=0}^{\infty}\mathbf{M}^{k+1}=\mathbb{I}+\sum_{k=1}^{\infty}\mathbf{M}^{k}-\sum_{k=1}^{\infty}\mathbf{M}^{k}=\mathbb{I}.$

∎

Since the matrix $\mathbf{P}$ is stochastic, we have
$\left\|\mathbf{P}\right\|_{\infty}=1$ and hence the spectral radius of
$\mathbf{P}$ is bounded by $1$. Since $\mathbf{P}_{TT}$ is a submatrix of
$\mathbf{P}$, we have $\left\|\mathbf{P}_{TT}\right\|_{\infty}\leq 1$ and its
spectral radius is also bounded by $1$. To prove Proposition 2.1 (denoted
Proposition A.5 below) we will show that the spectral radius of
$\mathbf{P}_{TT}$ is strictly smaller than $1$ if there is some vertex in $S$
that can be reached from any transient node via a directed path. Before
presenting the main proof, we require several lemmas.

######  Lemma A.2.

Let $\mathbf{B}$ and $\mathbf{C}$ be $n\times n$ matrices with non-negative
entries such that $\left\|\mathbf{B}\right\|_{\infty}\leq 1$ and
$\left\|\mathbf{C}\right\|_{\infty}\leq 1$ and let
$\mathbf{D}=\mathbf{C}\mathbf{B}$. Suppose there exists $1\leq p\leq n$ such
that $0<\sum_{j}B_{pj}<1$. Then, for every $1\leq i\leq n$ such that
$C_{ip}>0$,

$\sum_{j}D_{ij}<1.$

###### Proof.

Let $K=\\{k:C_{ik}>0\\}$. Then $p\in K$ and

$\displaystyle\sum_{j}D_{ij}$ $\displaystyle=\sum_{j}\sum_{k}C_{ik}B_{kj}$
$\displaystyle=\sum_{k\in K}C_{ik}\sum_{j}B_{kj}$ $\displaystyle\leq\sum_{k\in
K\setminus\\{p\\}}C_{ik}\left\|\mathbf{B}\right\|_{\infty}+C_{ip}\sum_{j}B_{pj}$
$\displaystyle<\sum_{k\in K\setminus\\{p\\}}C_{ik}+C_{ip}$ $\displaystyle\leq
1.$

∎

######  Lemma A.3.

Let $\Gamma$ be a weighted directed graph with weight matrix $\mathbf{W}$. Let
$i$ and $j$ be distinct nodes of $\Gamma$ connected by a directed path from
$i$ to $j$ of length $n\geq 1$. Then $W^{n}_{ij}>0$.

###### Proof.

We use induction. If $i$ and $j$ are connected with a path of length $1$, then
there exists an edge $(i,j)\in E$ and hence $W_{ij}>0$. Assume that
$W^{m}_{ij}>0$ if $i$ and $j$ are connected by a directed path from $i$ to $j$
of length $m$. Suppose $i$ and $j$ are connected by a path of length $m+1$.
Then there exists a vertex $k$ such that $i$ and $k$ are connected by a
directed path from $i$ to $k$ of length $m$ and there exists a directed edge
$(k,j)$. Hence, by our assumption $W^{m}_{ik}>0$ and $W_{kj}>0$. Therefore,

$W^{m+1}_{ij}=\sum_{k^{\prime}\in V}W^{m}_{ik^{\prime}}W_{k^{\prime}j}\geq
W^{m}_{ik}W_{kj}>0.$

∎

######  Lemma A.4.

Let $\mathbf{M}=\mathbf{P}_{TT}$, let $i\in T$ and suppose there exists $s\in
S$ such that there exists a directed path from $i$ to $s$ of length $m$. Then
for all $n\geq m$,

$\sum_{k\in T}M^{n}_{ik}<1.$ (36)

###### Proof.

Let $i\in T$ and let $s\in S$ be a vertex such that there exists a directed
path from $i$ to $s$ of length $m$. Let $J$ be the set of vertices in $T$
directly adjacent to a vertex in $S$. Then, by our assumption, for every $i\in
T$ there exists $j\in J$ such that there exists a directed path from $i$ to
$j$ of length $m-1$. Since the matrix $\mathbf{P}_{TT}$ can be treated as the
weight matrix for the subgraph of $\Gamma$ restricted to vertices in $T$, it
follows by Lemma A.3 that $M^{m-1}_{ij}>0$.

Since every point in $J$ is adjacent to a point in $S$, it also follows that
$\sum_{k\in T}M_{jk}<1$. Clearly, $\left\|\mathbf{M}\right\|_{\infty}\leq 1$
and hence $\left\|\mathbf{M}^{m-1}\right\|_{\infty}\leq 1$. Applying Lemma A.2
to the matrices $\mathbf{M}$ and $\mathbf{M}^{m-1}$ we obtain that for every
$i\in T$, $\sum_{k\in T}M^{m}_{ik}<1$.

Let $t\geq m$ and assume $\sum_{k\in T}M^{t}_{ik}<1$. We have

$\sum_{k\in T}M^{t+1}_{ik}=\sum_{k\in T}\sum_{k^{\prime}\in
T}M^{t}_{ik^{\prime}}M_{k^{\prime}k}=\sum_{k^{\prime}\in
T}M^{t}_{ik^{\prime}}\sum_{k\in T}M_{k^{\prime}k}\leq\sum_{k^{\prime}\in
T}M^{t}_{ik^{\prime}}<1$

and our result follows by induction. ∎

######  Proposition A.5.

Suppose that for every $p\in T$ there exists $s\in S$ such that there exists a
directed path from $p$ to $s$. Then, the matrix $\mathbb{I}-\mathbf{P}_{TT}$
is invertible and

$(\mathbb{I}-\mathbf{P}_{TT})^{-1}=\sum_{k=0}^{\infty}(\mathbf{P}_{TT})^{k}.$
(37)

###### Proof.

Let $\mathbf{M}=\mathbf{P}_{TT}$. Observe that our assumption implies that for
every $i\in T$ there exists $s\in S$ such that there exists a directed path
from $i$ to $s$ of length at most $N$. By Lemma A.4, we have for every $i\in
T$, $\sum_{k\in T}M^{N}_{ik}<1$. Hence,
$\left\|\mathbf{M}^{N}\right\|_{\infty}<1$ and therefore the spectral radius
of $\mathbf{M}=\mathbf{P}_{TT}$ is strictly smaller than $1$. Our result
follows by Lemma A.1. ∎

###  A.1 Information dissipation

######  Proposition A.6.

Let $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ be vectors of length $N$
such that for all $i\in V$, $\alpha_{i}>0$ and $\beta_{i}>0$. Define the
$N\times N$ matrix $\mathbf{\tilde{P}}$ with entries

$\tilde{P}_{ij}=\alpha_{i}\beta_{j}P_{ij},$

Let $\alpha_{*}=\max\\{\alpha_{i}:i\in V\\}$ and
$\beta_{*}=\max\\{\beta_{i}:i\in V\\}$ and suppose $\alpha_{*}\beta_{*}<1$.
Then, the matrix $\mathbb{I}-\mathbf{\tilde{P}}_{TT}$ is invertible and

$(\mathbb{I}-\mathbf{\tilde{P}}_{TT})^{-1}=\sum_{k=0}^{\infty}(\mathbf{\tilde{P}}_{TT})^{k}.$
(38)

###### Proof.

Let $\mathbf{M}=\mathbf{\tilde{P}}_{TT}$ and let $i\in T$. Then,

$\displaystyle\sum_{j\in T}M_{ij}$ $\displaystyle=\sum_{j\in
T}\alpha_{i}\beta_{j}P_{ij}\leq\alpha_{*}\beta_{*}\sum_{j\in T}P_{ij}<1.$

Hence, $\left\|\mathbf{M}\right\|_{\infty}<1$ and thus the spectral radius of
$\mathbf{\tilde{P}}_{TT}$ is strictly smaller than $1$. Our result then
follows by Lemma A.1. ∎

More generally, it is possible to interpret dissipation in the light of
Proposition A.5 by constructing a new graph $\tilde{\Gamma}$ with the vertex
set $\tilde{V}=V\cup\\{v\\}$, where $v$ denotes an additional vertex. The
weight matrix of $\tilde{\Gamma}$, denoted $\mathbf{\tilde{W}}$, has entries

$\tilde{W}_{ij}=\begin{cases}\alpha_{i}\beta_{j}P_{ij}&\text{if $i\in V$ and
$j\in V$,}\\\ 1-\sum_{k\in V}\alpha_{i}\beta_{k}P_{ik}&\text{if $i\in V$ and
$j=v$,}\\\ 0&\text{if $i=v$.}\end{cases}$ (39)

Clearly, a random walk on $\tilde{\Gamma}$ is equivalent to a random walk on
$\Gamma$ with dissipation: the dissipated information is directed towards the
additional vertex $v$ and then disappears. If we place $v$ in the boundary set
$\tilde{S}$, by Proposition A.5, the necessary condition for existence of the
Green’s function $(\mathbb{I}-\mathbf{\tilde{P}}_{TT})^{-1}$ is that from
every transient node $i$ there exists a directed path to either a node $s\in
S$ or a node $j\in T$ such that $\sum_{k\in V}\alpha_{j}\beta_{k}P_{jk}<1$
(such node $j$ is adjacent to $v$ in the graph $\tilde{\Gamma}$. Proposition
A.6 then just represents the special case where every transient vertex is
adjacent to $v$ in $\tilde{\Gamma}$.

##  Appendix B Interpretations of the matrices $\mathbf{F}$ and $\mathbf{H}$

###  B.1 $\mathbf{F}$ and $\mathbf{H}$ as matrices of expected visiting times

We will show that both $F_{ij}$ and $H_{ij}$ can be interpreted as the
expected number of times a random walk originating at the vertex $i$ visits
the vertex $j$, while avoiding all vertices in the boundary set $S$. Note that
in the case of the matrix $\mathbf{F}$, we have $i\in T$ and $j\in S$ while
for the matrix $\mathbf{H}$, $i\in S$ and $j\in T$. We will use $\mathbb{E}$
to denote the expectation operator.

######  Lemma B.1.

Suppose the boundary set $S$ represents sinks and let $Z_{ij}$ be a random
variable denoting the total number of times a random walk starting at $i\in T$
is absorbed at $j\in S$. Then,

$\mathbb{E}(Z_{ij})=F_{ij}.$ (40)

###### Proof.

Let $Y_{ij}(t)$ be the random variable taking the value $1$ if the random walk
originating at $i\in T$ is absorbed at $j\in S$ at time $t$, with probability
$\sum_{k\in T}P_{ik}^{t-1}P_{kj}$, and taking the value $0$ otherwise. We have
$Z_{ij}=\sum_{t=1}^{\infty}Y_{ij}(t)$ and $\mathbb{E}(Y_{ij}(t))=\sum_{k\in
T}P_{ik}^{t-1}P_{kj}$. Thus,

$\displaystyle\mathbb{E}(Z_{ij})$
$\displaystyle=\mathbb{E}\left(\sum_{t=1}^{\infty}Y_{ij}(t)\right)$
$\displaystyle=\sum_{t=1}^{\infty}\mathbb{E}(Y_{ij}(t))$
$\displaystyle=\sum_{t=1}^{\infty}\sum_{k\in T}P_{ik}^{t-1}P_{kj}$
$\displaystyle=\sum_{k\in T}\sum_{t=0}^{\infty}P_{ik}^{t}P_{kj}$
$\displaystyle=\sum_{k\in T}G_{ik}P_{kj}$ $\displaystyle=F_{ij}.\qed$

######  Lemma B.2.

Suppose the boundary set $S$ represents sources and let $Z_{ij}$ be a random
variable denoting the total number of times a random walk starting at $i\in S$
visits the node $j\in T$. Then,

$\mathbb{E}(Z_{ij})=H_{ij}.$ (41)

###### Proof.

In the same fashion as above, let $Y_{ij}(t)$ be the random variable taking
the value $1$ if the random walk originating at $i\in S$ is at $j\in T$ at
time $t$, with probability $\sum_{k\in T}P_{ik}P^{t-1}_{kj}$, and taking the
value $0$ otherwise. We have $Z_{ij}=\sum_{t=1}^{\infty}Y_{ij}(t)$ and
$\mathbb{E}(Y_{ij}(t))=\sum_{k\in T}P_{ik}P_{kj}^{t-1}$. Thus,

$\displaystyle\mathbb{E}(Z_{ij})$
$\displaystyle=\mathbb{E}\left(\sum_{t=1}^{\infty}Y_{ij}(t)\right)$
$\displaystyle=\sum_{t=1}^{\infty}\mathbb{E}(Y_{ij}(t))$
$\displaystyle=\sum_{t=1}^{\infty}\sum_{k\in T}P_{ik}P^{t-1}_{kj}$
$\displaystyle=\sum_{k\in T}\sum_{t=0}^{\infty}P_{ik}P^{t}_{kj}$
$\displaystyle=\sum_{k\in T}P_{ik}G_{kj}$ $\displaystyle=H_{ij}.\qed$

###  B.2 Invariants of $\mathbf{F}$ and $\mathbf{H}$

Let $\mathbf{1}\in\mathbb{R}^{n}$ denote the vector whose entries are all
$1$’s. Since all rows of $\mathbf{P}$ sum to unity, it follows that
$\mathbf{P}\mathbf{1}=\mathbf{1}$ and hence $\mathbf{1}$ is a right
eigenvector of $\mathbf{P}$ for the eigenvalue $\lambda=1$. Define
$\mathbf{d}$ as a vector of length $n$ having entries $d_{i}=\sum_{j}W_{ij}$.
If $\Gamma$ is unweighted graph, $d_{i}$ gives the degree of the node $i$.
Assuming $\mathbf{W}$ is symmetric,

$\sum_{k}P_{kj}d_{k}=\sum_{k}W_{kj}=\sum_{k}W_{jk}=d_{j}$

and therefore $\mathbf{d}$ is a left eigenvector of $\mathbf{P}$ corresponding
to the eigenvalue $\lambda=1$. This leads to the following result.

######  Lemma B.3.

Suppose that the matrix $\mathbb{I}-\mathbf{P}_{TT}$ is invertible. Let
$\mathbf{u}$ and $\mathbf{v}$ be the left and right eigenvector of the matrix
$\mathbf{P}$ corresponding to the eigenvalue $\lambda=1$, respectively. Write
$\mathbf{u}=[\mathbf{u}_{S}\ \mathbf{u}_{T}]$ and
$\mathbf{v}=\left[\begin{array}[]{c}\mathbf{v}_{S}\\\
\mathbf{v}_{T}\end{array}\right]$. Then,

$\mathbf{u}_{T}=\mathbf{u}_{S}\mathbf{H},$ (42)

and

$\mathbf{v}_{T}=\mathbf{F}\mathbf{v}_{S}.$ (43)

###### Proof.

Using the canonical form of the matrix $\mathbf{P}$ (Equation (3)) and the
fact that $\mathbf{u}$ and $\mathbf{v}$ are left and right eigenvectors of
$\mathbf{P}$ respectively, we obtain

$\mathbf{u}_{T}=\mathbf{u}_{S}\mathbf{P}_{ST}+\mathbf{u}_{T}\mathbf{P}_{TT},$
(44)

and

$\mathbf{v}_{T}=\mathbf{P}_{TS}\mathbf{v}_{S}+\mathbf{P}_{TT}\mathbf{v}_{T}.$
(45)

Rearranging Equations (44) and (45) leads to

$\mathbf{u}_{T}(\mathbb{I}-\mathbf{P}_{TT})=\mathbf{u}_{S}\mathbf{P}_{ST},$
(46)

and

$(\mathbb{I}-\mathbf{P}_{TT})\mathbf{v}_{T}=\mathbf{P}_{TS}\mathbf{v}_{S}.$
(47)

Our result then follows as the consequence of invertibility of
$\mathbb{I}-\mathbf{P}_{TT}$. ∎

Since $\mathbf{1}$ is a right eigenvector of $\mathbf{P}$, it follows from
(43) that for all $i$, $\sum_{j\in S}F_{ij}=1$. Furthemore, recall that if
$\Gamma$ is an undirected graph, $\mathbf{W}$ is symmetric and $\mathbf{d}$ is
a left eigenvector of $\mathbf{P}$ for $\lambda=1$. Assuming the matrix
$\mathbf{H}$ exists, we obtain from Lemma B.3 that, if $S$ contains a single
point, the matrix $\mathbf{H}$ is a row vector, which is a multiple of
$\mathbf{d}_{T}$.

## References

  * Aittokallio and Schwikowski (2006) Aittokallio, T. and Schwikowski, B. 2006.  Graph-based methods for analysing networks in cell biology.  _Brief Bioinform_ 7, 243–255.
  * Allard _et al._ (1999) Allard, S., Utley, R. T., Savard, J. _et al._ 1999\.  NuA4, an essential transcription adaptor/histone H4 acetyltransferase complex containing Esa1p and the ATM-related cofactor Tra1p.  _EMBO J_ 18, 5108–5119.
  * Althoefer _et al._ (1995) Althoefer, H., Schleiffer, A., Wassmann, K. _et al._ 1995\.  Mcm1 is required to coordinate G2-specific transcription in Saccharomyces cerevisiae.  _Mol Cell Biol_ 15, 5917–5928.
  * Amaral _et al._ (2000) Amaral, L. A., Scala, A., Barthelemy, M. _et al._ 2000\.  Classes of small-world networks.  _Proc Natl Acad Sci U S A_ 97, 11149–11152.
  * Bapat and Raghavan (1997) Bapat, R. B. and Raghavan, T. E. S. 1997.  _Nonnegative matrices and applications_ , _Encyclopedia of Mathematics and its Applications_ , volume 64.  Cambridge University Press, Cambridge.
  * Barthelemy _et al._ (2005) Barthelemy, M., Barrat, A., Pastor-Satorras, R. _et al._ 2005\.  Dynamical patterns of epidemic outbreaks in complex heterogeneous networks.  _J Theor Biol_ 235, 275–288.
  * Biddick and Young (2005) Biddick, R. and Young, E. 2005.  Yeast mediator and its role in transcriptional regulation.  _C R Biol_ 328, 773–82.
  * Borrow _et al._ (1996) Borrow, J., Stanton, V. P. J., Andresen, J. M. _et al._ 1996\.  The translocation t(8;16)(p11;p13) of acute myeloid leukaemia fuses a putative acetyltransferase to the CREB-binding protein.  _Nat Genet_ 14, 33–41.
  * Broder _et al._ (2000) Broder, A., Kumar, R., Maghoul, F. _et al._ 2000\.  Graph structure in the web.  _Comput. Networks_ 33, 309–320.
  * Brown _et al._ (2001) Brown, C. E., Howe, L., Sousa, K. _et al._ 2001\.  Recruitment of HAT complexes by direct activator interactions with the ATM-related Tra1 subunit.  _Science_ 292, 2333–2337.
  * Bruhn _et al._ (1992) Bruhn, L., Hwang-Shum, J. J., and Sprague, G. F. J. 1992.  The N-terminal 96 residues of MCM1, a regulator of cell type-specific genes in Saccharomyces cerevisiae, are sufficient for DNA binding, transcription activation, and interaction with alpha 1.  _Mol Cell Biol_ 12, 3563–3572.
  * Cancho _et al._ (2001) Cancho, R. F. i., Janssen, C., and Solé, R. V. 2001.  Topology of technology graphs: Small world patterns in electronic circuits.  _Phys. Rev. E_ 64, 046119.
  * Chen and Yuan (2006) Chen, J. and Yuan, B. 2006.  Detecting functional modules in the yeast protein–protein interaction network.  _Bioinformatics_ 22, 2283–2290.
  * Chua _et al._ (2006) Chua, H. N., Sung, W.-K., and Wong, L. 2006.  Exploiting indirect neighbours and topological weight to predict protein function from protein-protein interactions.  _Bioinformatics_ 22, 1623–1630.
  * Chung and Yau (2000) Chung, F. and Yau, S.-T. 2000.  Discrete Green’s functions.  _J. Combin. Theory Ser. A_ 91, 191–214.
  * Chung (1997) Chung, F. R. K. 1997.  _Spectral graph theory_ , _CBMS Regional Conference Series in Mathematics_ , volume 92.  Published for the Conference Board of the Mathematical Sciences, Washington, DC.
  * Collins _et al._ (2007) Collins, S. R., Kemmeren, P., Zhao, X.-C. _et al._ 2007\.  Toward a Comprehensive Atlas of the Physical Interactome of Saccharomyces cerevisiae.  _Mol Cell Proteomics_ 6, 439–450.
  * Davis (2004) Davis, T. A. 2004.  Algorithm 832: Umfpack v4.3—an unsymmetric-pattern multifrontal method.  _ACM Trans. Math. Softw._ 30\.
  * Deane _et al._ (2002) Deane, C. M., Salwinski, L., Xenarios, I. _et al._ 2002\.  Protein interactions: two methods for assessment of the reliability of high throughput observations.  _Mol Cell Proteomics_ 1, 349–356.
  * Doyon and Cote (2004) Doyon, Y. and Cote, J. 2004.  The highly conserved and multifunctional NuA4 HAT complex.  _Curr Opin Genet Dev_ 14, 147–154.
  * Dubois and Messenguy (1991) Dubois, E. and Messenguy, F. 1991.  In vitro studies of the binding of the ARGR proteins to the ARG5,6 promoter.  _Mol Cell Biol_ 11, 2162–2168.
  * Eberharter _et al._ (1999) Eberharter, A., Sterner, D. E., Schieltz, D. _et al._ 1999\.  The ADA complex is a distinct histone acetyltransferase complex in Saccharomyces cerevisiae.  _Mol Cell Biol_ 19, 6621–6631.
  * Galitski (2004) Galitski, T. 2004.  Molecular networks in model systems.  _Annu Rev Genomics Hum Genet_ 5, 177–187.
  * Gansner and North (2000) Gansner, E. R. and North, S. C. 2000.  An open graph visualization system and its applications to software engineering.  _Software — Practice and Experience_ 30, 1203–1233.
  * Gavin _et al._ (2006) Gavin, A.-C., Aloy, P., Grandi, P. _et al._ 2006\.  Proteome survey reveals modularity of the yeast cell machinery.  _Nature_ 440, 631–636.
  * Giot _et al._ (2003) Giot, L., Bader, J. S., Brouwer, C. _et al._ 2003\.  A protein interaction map of Drosophila melanogaster.  _Science_ 302, 1727–1736.
  * Girvan and Newman (2002) Girvan, M. and Newman, M. E. J. 2002.  Community structure in social and biological networks.  _Proc Natl Acad Sci U S A_ 99, 7821–7826.
  * Grant _et al._ (1998) Grant, P., Schieltz, D., Pray-Grant, M. _et al._ 1998\.  The ATM-related cofactor Tra1 is a component of the purified SAGA complex.  _Mol Cell_ 2, 863–7.
  * Grant _et al._ (1997) Grant, P. A., Duggan, L., Cote, J. _et al._ 1997\.  Yeast Gcn5 functions in two multisubunit complexes to acetylate nucleosomal histones: characterization of an Ada complex and the SAGA (Spt/Ada) complex.  _Genes Dev_ 11, 1640–1650.
  * Guelzim _et al._ (2002) Guelzim, N., Bottani, S., Bourgine, P. _et al._ 2002\.  Topological and causal structure of the yeast transcriptional regulatory network.  _Nat Genet_ 31, 60–63.
  * Hartwell _et al._ (1999) Hartwell, L. H., Hopfield, J. J., Leibler, S. _et al._ 1999\.  From molecular to modular cell biology.  _Nature_ 402, 47–52.
  * Hinnebusch (2005) Hinnebusch, A. 2005.  Translational regulation of GCN4 and the general amino acid control of yeast.  _Annu Rev Microbiol_ 59, 407–50.
  * Hoffman and Kunze (1971) Hoffman, K. and Kunze, R. 1971.  _Linear algebra_.  Second edition. Prentice-Hall Inc., Englewood Cliffs, N.J.
  * Hwang _et al._ (2006) Hwang, W., Cho, Y.-R., Zhang, A. _et al._ 2006\.  A novel functional module detection algorithm for protein-protein interaction networks.  _Algorithms Mol Biol_ 1, 24.
  * Ito _et al._ (2001) Ito, T., Chiba, T., Ozawa, R. _et al._ 2001\.  A comprehensive two-hybrid analysis to explore the yeast protein interactome.  _Proc Natl Acad Sci U S A_ 98, 4569–4574.
  * Ito _et al._ (2000) Ito, T., Tashiro, K., Muta, S. _et al._ 2000\.  Toward a protein-protein interaction map of the budding yeast: A comprehensive system to examine two-hybrid interactions in all possible combinations between the yeast proteins.  _Proc Natl Acad Sci U S A_ 97, 1143–1147.
  * Jones _et al._ (2001–) Jones, E., Oliphant, T., Peterson, P. _et al._ 2001–.  SciPy: Open source scientific tools for Python.  URL http://www.scipy.org/.
  * Kassir _et al._ (2003) Kassir, Y., Adir, N., Boger-Nadjar, E. _et al._ 2003\.  Transcriptional regulation of meiosis in budding yeast.  _Int Rev Cytol_ 224, 111–71.
  * Kelley _et al._ (2003) Kelley, B. P., Sharan, R., Karp, R. M. _et al._ 2003\.  Conserved pathways within bacteria and yeast as revealed by global protein network alignment.  _Proc Natl Acad Sci U S A_ 100, 11394–11399.
  * Kemeny and Snell (1976) Kemeny, J. G. and Snell, J. L. 1976.  _Finite Markov chains_.  Springer-Verlag, New York.  Reprinting of the 1960 original, Undergraduate Texts in Mathematics.
  * Krogan _et al._ (2006) Krogan, N. J., Cagney, G., Yu, H. _et al._ 2006\.  Global landscape of protein complexes in the yeast Saccharomyces cerevisiae.  _Nature_ 440, 637–643.
  * Krogan and Greenblatt (2001) Krogan, N. J. and Greenblatt, J. F. 2001.  Characterization of a six-subunit holo-elongator complex required for the regulated expression of a group of genes in Saccharomyces cerevisiae.  _Mol Cell Biol_ 21, 8203–8212.
  * Li _et al._ (2004) Li, S., Armstrong, C. M., Bertin, N. _et al._ 2004\.  A map of the interactome network of the metazoan C. elegans.  _Science_ 303, 540–543.
  * Lorch _et al._ (2000) Lorch, Y., Beve, J., Gustafsson, C. _et al._ 2000\.  Mediator-nucleosome interaction.  _Mol Cell_ 6, 197–201.
  * Ma and Zeng (2003) Ma, H. and Zeng, A.-P. 2003.  Reconstruction of metabolic networks from genome data and analysis of their global structure for various organisms.  _Bioinformatics_ 19, 270–277.
  * Martinez (2002) Martinez, E. 2002.  Multi-protein complexes in eukaryotic gene transcription.  _Plant Mol Biol_ 50, 925–947.
  * Mead _et al._ (2002) Mead, J., Bruning, A. R., Gill, M. K. _et al._ 2002\.  Interactions of the Mcm1 MADS box protein with cofactors that regulate mating in yeast.  _Mol Cell Biol_ 22, 4607–4621.
  * Messenguy and Dubois (1993) Messenguy, F. and Dubois, E. 1993.  Genetic evidence for a role for MCM1 in the regulation of arginine metabolism in Saccharomyces cerevisiae.  _Mol Cell Biol_ 13, 2586–2592.
  * Miklos and Maleszka (2004) Miklos, G. L. G. and Maleszka, R. 2004.  Microarray reality checks in the context of a complex disease.  _Nat Biotechnol_ 22, 615–621.
  * Mizzen _et al._ (1996) Mizzen, C. A., Yang, X. J., Kokubo, T. _et al._ 1996\.  The TAF(II)250 subunit of TFIID has histone acetyltransferase activity.  _Cell_ 87, 1261–1270.
  * Myer and Young (1998) Myer, V. E. and Young, R. A. 1998.  RNA polymerase II holoenzymes and subcomplexes.  _J Biol Chem_ 273, 27757–27760.
  * Nabieva _et al._ (2005) Nabieva, E., Jim, K., Agarwal, A. _et al._ 2005\.  Whole-proteome prediction of protein function via graph-theoretic analysis of interaction maps.  _Bioinformatics_ 21 Suppl 1, 302–310.
  * Neuwald and Landsman (1997) Neuwald, A. F. and Landsman, D. 1997.  GCN5-related histone N-acetyltransferases belong to a diverse superfamily that includes the yeast SPT10 protein.  _Trends Biochem Sci_ 22, 154–155.
  * Newman (2004) Newman, M. E. J. 2004.  Coauthorship networks and patterns of scientific collaboration.  _Proc Natl Acad Sci U S A_ 101 Suppl 1, 5200–5205.
  * Osada _et al._ (2001) Osada, S., Sutton, A., Muster, N. _et al._ 2001\.  The yeast SAS (something about silencing) protein complex contains a MYST-type putative acetyltransferase and functions with chromatin assembly factor ASF1.  _Genes Dev_ 15, 3155–3168.
  * Papamichos-Chronakis _et al._ (2002) Papamichos-Chronakis, M., Petrakis, T., Ktistaki, E. _et al._ 2002\.  Cti6, a PHD domain protein, bridges the Cyc8-Tup1 corepressor and the SAGA coactivator to overcome repression at GAL1.  _Mol Cell_ 9, 1297–1305.
  * Pellegrini _et al._ (2004) Pellegrini, M., Haynor, D., and Johnson, J. M. 2004.  Protein interaction networks.  _Expert Rev Proteomics_ 1, 239–249.
  * Pereira-Leal _et al._ (2004) Pereira-Leal, J. B., Enright, A. J., and Ouzounis, C. A. 2004.  Detection of functional modules from protein interaction networks.  _Proteins_ 54, 49–57.  Evaluation Studies.
  * Reggiori _et al._ (2003) Reggiori, F., Wang, C.-W., Stromhaug, P. E. _et al._ 2003\.  Vps51 is part of the yeast Vps fifty-three tethering complex essential for retrograde traffic from the early endosome and Cvt vesicle completion.  _J Biol Chem_ 278, 5009–5020.
  * Ricci _et al._ (2002) Ricci, A., Genereaux, J., and Brandl, C. 2002.  Components of the SAGA histone acetyltransferase complex are required for repressed transcription of ARG1 in rich medium.  _Mol Cell Biol_ 22, 4033–42.
  * Rives and Galitski (2003) Rives, A. and Galitski, T. 2003.  Modular organization of cellular networks.  _Proc Natl Acad Sci U S A_ 100, 1128–33.
  * Rual _et al._ (2005) Rual, J.-F., Venkatesan, K., Hao, T. _et al._ 2005\.  Towards a proteome-scale map of the human protein-protein interaction network.  _Nature_ 437, 1173–1178.
  * Salwinski _et al._ (2004) Salwinski, L., Miller, C. S., Smith, A. J. _et al._ 2004\.  The Database of Interacting Proteins: 2004 update.  _Nucleic Acids Res_ 32, 449–451.
  * Sanders and Weil (2000) Sanders, S. L. and Weil, P. A. 2000.  Identification of two novel TAF subunits of the yeast Saccharomyces cerevisiae TFIID complex.  _J Biol Chem_ 275, 13895–13900.
  * Sharan _et al._ (2005a) Sharan, R., Ideker, T., Kelley, B. _et al._ 2005a.  Identification of protein complexes by comparative analysis of yeast and bacterial protein interaction data.  _J Comput Biol_ 12, 835–846.
  * Sharan _et al._ (2005b) Sharan, R., Suthram, S., Kelley, R. M. _et al._ 2005b.  Conserved patterns of protein interaction in multiple species.  _Proc Natl Acad Sci U S A_ 102, 1974–1979.
  * Shore and Sharrocks (1995) Shore, P. and Sharrocks, A. D. 1995.  The MADS-box family of transcription factors.  _Eur J Biochem_ 229, 1–13.
  * Smith _et al._ (1998) Smith, E. R., Eisen, A., Gu, W. _et al._ 1998\.  ESA1 is a histone acetyltransferase that is essential for growth in yeast.  _Proc Natl Acad Sci U S A_ 95, 3561–3565.
  * Spirin and Mirny (2003) Spirin, V. and Mirny, L. 2003.  Protein complexes and functional modules in molecular networks.  _Proc Natl Acad Sci U S A_ 100, 12123–8.
  * Sprinzak _et al._ (2003) Sprinzak, E., Sattath, S., and Margalit, H. 2003.  How reliable are experimental protein-protein interaction data?  _J Mol Biol_ 327, 919–923.
  * Steffen _et al._ (2002) Steffen, M., Petti, A., Aach, J. _et al._ 2002\.  Automated modelling of signal transduction networks.  _BMC Bioinformatics_ 3, 34.
  * Stelzl _et al._ (2005) Stelzl, U., Worm, U., Lalowski, M. _et al._ 2005\.  A human protein-protein interaction network: a resource for annotating the proteome.  _Cell_ 122, 957–968.
  * Sterner and Berger (2000) Sterner, D. E. and Berger, S. L. 2000.  Acetylation of histones and transcription-related factors.  _Microbiol Mol Biol Rev_ 64, 435–459.
  * Sternglanz and Schindelin (1999) Sternglanz, R. and Schindelin, H. 1999.  Structure and mechanism of action of the histone acetyltransferase Gcn5 and similarity to other N-acetyltransferases.  _Proc Natl Acad Sci U S A_ 96, 8807–8808.  Comment.
  * Struhl (1998) Struhl, K. 1998.  Histone acetylation and transcriptional regulatory mechanisms.  _Genes Dev_ 12, 599–606.
  * Suthram _et al._ (2006) Suthram, S., Shlomi, T., Ruppin, E. _et al._ 2006\.  A direct comparison of protein interaction confidence assignment schemes.  _BMC Bioinformatics_ 7, 360.
  * Thouless (1974) Thouless, D. J. 1974.  Electrons in disordered systems and the theory of localization.  _Phys Rep C_ 13, 93–142.
  * Timmers and Tora (2005) Timmers, H. T. M. and Tora, L. 2005.  SAGA unveiled.  _Trends Biochem Sci_ 30, 7–10.
  * Traven _et al._ (2006) Traven, A., Jelicic, B., and Sopta, M. 2006.  Yeast Gal4: a transcriptional paradigm revisited.  _EMBO Rep_ 7, 496–9.
  * Uetz _et al._ (2000) Uetz, P., Giot, L., Cagney, G. _et al._ 2000\.  A comprehensive analysis of protein-protein interactions in Saccharomyces cerevisiae.  _Nature_ 403, 623–627.
  * Wasserman and Faust (1994) Wasserman, S. and Faust, K. 1994.  _Social network analysis_.  Cambridge University Press, Cambridge.
  * Whaley and Petitet (2005) Whaley, R. C. and Petitet, A. 2005.  Minimizing development and maintenance costs in supporting persistently optimized BLAS.  _Software: Practice and Experience_ 35, 101–121.
  * Wittschieben _et al._ (1999) Wittschieben, B. O., Otero, G., de Bizemont, T. _et al._ 1999\.  A novel histone acetyltransferase is an integral subunit of elongating RNA polymerase II holoenzyme.  _Mol Cell_ 4, 123–128.
  * Wolffe (1992) Wolffe, A. P. 1992.  _Chromatin Structure and Function_.  Academic Press, London, England.
  * Wolffe (2001) Wolffe, A. P. 2001.  Transcriptional regulation in the context of chromatin structure.  _Essays Biochem_ 37, 45–57.
  * Workman and Kingston (1998) Workman, J. L. and Kingston, R. E. 1998.  Alteration of nucleosome structure as a mechanism of transcriptional regulation.  _Annu Rev Biochem_ 67, 545–579.
  * Xenarios _et al._ (2002) Xenarios, I., Salwinski, L., Duan, X. J. _et al._ 2002\.  DIP, the Database of Interacting Proteins: a research tool for studying cellular networks of protein interactions.  _Nucleic Acids Res_ 30, 303–305.
  * Xiong _et al._ (2005) Xiong, H., He, X., Ding, C. _et al._ 2005\.  Identification of functional modules in protein complexes via hyperclique pattern discovery.  _Pac Symp Biocomput_ 221–232.
  * Zhang _et al._ (2007) Zhang, Y.-C., Blattner, M., and Yu, Y.-K. 2007.  Heat conduction process on community networks as a recommendation model.  Submitted to Phys. Rev. Lett.