### Simulation study

Simulation utilized the sample expression data gal80R given in Cytoscape (http://cytoscape.org/). There are 331 genes and 361 interactions in this network. Within it, we randomly selected subnetworks at three different sizes n (n = 40, 60, 80), as condition-responsive. In each responsive subnetwork m% (80%, 90%, 100%) of genes are defined to be active. The significance values of active genes were assigned randomly with top n\times m\% significance values in gal80R, and that of the other genes were randomly sampled from the rest of the significance values. The phase locking index λ (see 2.3) of the interactions in the predefined responsive subnetwork were sampled from N\left(0.8,0.5\right), i.e. a normal distribution with μ = 0.8, σ = 0.5; while λ for the remaining edges were sampled from N\left(0.4,0.3\right). The choice of these values was based on the distribution of the λ values of gene pairs in protein complexes and of randomly selected gene pairs. For protein complexes we used the MIPS annotation (http://mips.helmholtz-muenchen.de/genre/proj/yeast) edited by Gerstein Lab (http://www.gersteinlab.org/proj/bottleneck/mips.txt).

A gene of the predefined responsive subnetworks that is in the TopoPL-identified subnetwork is considered a successful identification. This procedure was repeated 10 times and the true positive (TP, sensitivity) rate was defined to be the number of successful identifications divided by the size of the predefined network n. The false positive (FP, specificity or precision) rate was estimated as the number of false identifications divided by the size of the identified subnetwork. The F score is a measure of a test's accuracy. It considers both the precision and the sensitivity of the test:

F=\frac{specificity\mathsf{\text{*}}sensitivity}{\left(specificity+sensitivity\right)/2}

We used the average sensitivity, specificity and F score to measure the performance of TopoPL. The performance is also evaluated with Receiver Operating Characteristic (ROC) curve, a plot of the true positive rate against the false positive rate [11].

### Gene expression and protein-protein interaction data

Gene expression data was downloaded from EMBL's Huber group (http://www.ebi.ac.uk/huber-srv/scercycle/). It is a time course study of yeast cell cycle, where cells were arrested using alpha factor or cdc28. The alpha factor dataset contains 41 time points and the cdc28 dataset contains 44 time points, both at 5-minute resolution. These datasets provide strand-specific profiles of temporal expression during the mitotic cell cycle of S. cerevisiae, monitored for more than three complete cell divisions [14]. Yeast PPI data were downloaded from BioGRID (thebiogrid.org, version 3.1.69).

### Phase locking analysis

The details of definitions and steps of the phase locking analysis was described in our previous work [11] and briefly summarized here. Given a time series *s*(*t*), its Hilbert transformation is given by

{s}_{H}\left(t\right)=\frac{1}{\pi}\mathsf{\text{PV}}{\int}_{-\infty}^{\infty}\frac{s\left(t\right)}{t-\tau}d\tau

(1)

where PV stand for Cauchy Principal Value of integration. The corresponding analytical signal can then be constructed by:

s\left(t\right)+i{s}_{H}\left(t\right)=A\left(t\right){e}^{i\phi \left(t\right)}

(2)

where the instantaneous phase \phi \left(t\right) is thus uniquely determined. For two time series with instantaneous phase {\phi}_{i}\left(t\right) and {\phi}_{j}\left(t\right), their cyclic relative phase is determined by

\mathrm{\Psi}\left(t\right)=\left({\phi}_{i}\left(t\right)-{\phi}_{j}\left(t\right)\right)\mathsf{\text{mod}}\left(2\pi \right)

(3)

If two time series interact with each other, there will be rhythmic adjustment resulting in phase locking: \mathrm{\Psi}={\mathrm{\Psi}}_{0} is a constant. To evaluate the significance of phase locking, we utilize the circular mean of the phase difference

\lambda =\left|\mathsf{\text{exp}}\left(i\mathrm{\Psi}\left(t\right)\right)\right|=\left|\left(\left(\frac{1}{{t}_{N}}\right)\sum _{l=1}^{N}\mathsf{\text{exp}}\left(i\mathrm{\Psi}\left({t}_{l}\right)\right)\right)\right|

(4)

In a perfect locking \lambda =\left|\mathsf{\text{exp}}\left(i{\mathrm{\Psi}}_{0}\right)\right|=1, and \lambda \to 0 when \mathrm{\Psi}\left(t\right) is randomly distributed. *λ* offers a new measure to infer potential interaction between gene pairs [11].

### TopoPL

For each gene *i*, the EDGE software [15] was used to calculate {p}_{i}, the significance of its expression changes during the time course study. We convert {p}_{i} to a z-score through {z}_{i}={\varnothing}^{-1}\left(1-{p}_{i}\right), where {\varnothing}^{-1} is the inverse normal CDF. Let {A}^{\left(P\right)}=\left({{a}^{\left(p\right)}}_{ij}\right) be the adjacency matrix of genes in a PPI subnetwork and A=\left({a}_{ij}\right)=\left({{a}^{\left(p\right)}}_{ij}*{\lambda}_{ij}\right), TopoPL defines the overall activity of a subnetwork with:

{z}_{A}^{TopoPL}={\sum}_{i\in A}{\sum}_{j\in A}|{z}_{i}{|}^{0.5}*{a}_{ij}*|{z}_{j}{|}^{0.5}*sgn\left({z}_{i}+{z}_{j}\right)

(5)

{z}_{A}^{TopoPL} captures the dynamic topological property of the subnetwork, and hub genes (genes with high network degree) contribute more to this metric. |{z}_{i}{|}^{0.5}*{a}_{ij}*|{z}_{j}{|}^{0.5}*sgn\left({z}_{i}+{z}_{j}\right),i\ne j can be regarded as the "activity measurement" of the interaction. Gene pairs with significant and synchronized expression changes, and whose gene products interact, contribute more to the activity of the subnetwork.

This metric is an improved version over the PCI that we previously proposed to identify active pathways from gene expression data [13]: PCI={\sum}_{i,j}|{x}_{is}{|}^{0.5}*{a}_{ij}*|{x}_{js}{|}^{0.5}*sgn\left({x}_{is}+{x}_{js}\right), where {x}_{is} is normalized log expression measurement of gene *i* in sample *s*, and ({a}_{ij}) is the adjacency matrix of the PPI network of genes in the pathway. The merit of PCI has been demonstrated in previous works [13]. To reduce the potential impact on the network measure from residual inter-sample and inter-array biases after normalization, here we adopted the non-parametric measure {z}_{i} in place of {x}_{is}. A similar metric to Eq. (5) was developed recently by us to predict candidate disease genes for type 1 diabetes, where {z}_{i} is the z-score of disease relevance of gene *i*. There again we demonstrated the advantage of incorporating network structural information [16].

Obviously, {z}_{A}^{TopoPL} increase with the number of nodes and edges. To adjust for network size and density, we use the following equation

{z}_{A}^{TopoPL}\to {z}_{A}^{TopoPL}*\frac{1}{\left(\#nodes+\#edges\right)}

(6)

We implemented the searching procedure based on simulated annealing. The pseudocode of the algorithm is described below:

**Input:** the entire network {G}_{0}=\left(V,E\right); a set of parameters for running simulated annealing: start temperature {T}_{start} (= 1 in this study), end temperature {T}_{end} (= 1e-8 in this study), number of iterations N.

**Output:** the subnetwork with the highest score.

**Steps**: initialize each node with its expression significance score {z}_{i} and each edge with its phase locking index; select the largest connected component (subnetwork) {G}_{out} from top 10% significant nodes of {G}_{0}; calculate score of {G}_{out} and obtain its score {z}_{out}^{TopoPL}; then run the following:

For i = 1 to N, Do

Calculate the current temperature {T}_{i}={T}_{i}*0.{8}^{1/N}; {G}_{try}\leftarrow {G}_{ou{t}^{\prime}}

Exit loop if {T}_{i}<{T}_{end}

Randomly pick a node n\in V

IF (n\in {G}_{try}), remove n from {G}_{try};

ELSE add n to {G}_{try};

Calculate score {z}_{try}^{TopoPL}for the largest connected component of {G}_{try};

Calculate \mathrm{\Delta}={z}_{try}^{TopoPL}-{z}_{out}^{TopoPL};

IF Δ> 0, then{G}_{out}\leftarrow {G}_{try};

ELSE, accept {G}_{out}\leftarrow {G}_{try} with the probabilityp={e}^{\mathrm{\Delta}/{T}_{i}};

END

These steps can be iterated to identify subnetworks with the next highest scores and so on.

In this study we compared TopoPL with two other methods: (1) The commonly used network scoring method that sums significance levels of all genes in the network (hereafter referred to as the Additive scoring method):

{z}_{Additive}={\sum}_{i\in A}{z}_{i}

(7)

(2) A metric that we previously proposed in our TAPPA software package [13] (hereafter referred to as the TAPPA scoring method) that only utilize the topological characteristics of the PPI network:

{z}_{Topo}={\sum}_{i\in A}{\sum}_{j\in A}|{z}_{i}{|}^{0.5}*{{a}^{\left(p\right)}}_{ij}*|{z}_{j}{|}^{0.5}*sgn\left({z}_{i}+{z}_{j}\right)

(8)