Community detection for networks with unipartite and bipartite structure

Because our model is an extension and generalization of BKN, we briefly review this generative model for a unipartite undirected network in this subsection.

Let ${{G}_{U}}=(V,E)$ denote a unipartite undirected network, where V is the vertex set and E is the edge set. Here, A represents the adjacency matrix, where ${{A}_{ij}}=1$ if there is a link between vertex i and vertex j, and ${{A}_{ii}}=2$ if there is a self-loop of vertex i. Further, suppose the number of modules K is given. BKN is parameterized by a set of parameters θ, such that ${{\theta }_{iz}}$ represents the propensity of vertex i to belong to module z. Specifically, the physical meaning of ${{\theta }_{iz}}{{\theta }_{jz}}$ is the expected number of links between vertices i and j belonging to module z, with the exact number being Poisson distributed. By summing over all K modules, $\sum z {{\theta }_{iz}}{{\theta }_{jz}}$ represents the expected number between vertices i and j, the exact number of which is also Poisson distributed, because the sum of independent Poisson variables is still Poisson distributed. Given that the actual number of edges between i and j is known, the likelihood function of the unipartite undirected network ${{G}_{U}}$ with the adjacency matrix A is

$$\begin{eqnarray}\begin{array}{lllllllllllllll} P\left( {{G}_{U}}|\theta \right) & = & \mathop \prod \limits_{i\lt j} \frac{{{\left( \mathop \sum \limits_{z} {{\theta }_{iz}}{{\theta }_{jz}} \right)}^{{{A}_{ij}}}}}{{{A}_{ij}}!}{\rm exp} \left( -\mathop \sum \limits_{z} {{\theta }_{iz}}{{\theta }_{jz}} \right) \\ {} & {} & \times \mathop \prod \limits_{i} \frac{{{\left( \frac{1}{2}\mathop \sum \limits_{z} {{\theta }_{iz}}{{\theta }_{iz}} \right)}^{{{A}_{ii}}/2}}}{\left( {{A}_{ii}}/2 \right)!}{\rm exp} \left( -\frac{1}{2}\mathop \sum \limits_{z} {{\theta }_{iz}}{{\theta }_{iz}} \right). \\ \end{array}\end{eqnarray} \tag{ 1 }$$

BKN can be fitted to an observed network according to the maximum likelihood principle, with respect to parameters ${{\theta }_{iz}}$. By taking the logarithm of equation (1), and discarding constant terms, the log likelihood function of the model can be written as

$$\begin{eqnarray}&&{\rm ln} P\left( {{G}_{U}}|\theta \right)=\mathop \sum \limits_{ij} {{A}_{ij}}{\rm ln} \left( \mathop \sum \limits_{z} {{\theta }_{iz}}{{\theta }_{jz}} \right)-\mathop \sum \limits_{ijz} {{\theta }_{iz}}{{\theta }_{jz}}.\end{eqnarray} \tag{ 2 }$$

Applying Jensen's inequality, we obtain

$$\begin{eqnarray}&&{\mkern 1mu} {\rm ln} P({{G}_{U}}|\theta ){\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} \geqslant \mathop \sum \limits_{ijz} \left[ {{A}_{ij}}{{q}_{ij}}(z){\rm ln} \left( \frac{{{\theta }_{iz}}{{\theta }_{jz}}}{{{q}_{ij}}(z)} \right)-{{\theta }_{iz}}{{\theta }_{jz}} \right]\end{eqnarray} \tag{ 3 }$$

by introducing an arbitrary variable ${{q}_{ij}}(z)$, which satisfies $\sum z {{q}_{ij}}(z)=1$. When

$$\begin{eqnarray}&&{{q}_{ij}}(z)=\frac{{{\theta }_{iz}}{{\theta }_{jz}}}{\mathop \sum \limits_{z} {{\theta }_{iz}}{{\theta }_{jz}}},\end{eqnarray} \tag{ 4 }$$

the inequality becomes an exact equality. By differentiating equation (3), with respect to ${{\theta }_{iz}}$, we obtain

$$\begin{eqnarray}&&{{\theta }_{iz}}=\frac{\mathop \sum \limits_{j} {{A}_{ij}}{{q}_{ij}}(z)}{\mathop \sum \limits_{j} {{\theta }_{jz}}},\end{eqnarray} \tag{ 5 }$$

which gives the optimal value of ${{\theta }_{iz}}$. In [6], Ball et al further change the form of equation (5) to

$$\begin{eqnarray}&&{{\theta }_{iz}}=\frac{\mathop \sum \limits_{j} {{A}_{ij}}{{q}_{ij}}(z)}{\sqrt{\mathop \sum \limits_{ij} {{A}_{ij}}{{q}_{ij}}(z)}},\end{eqnarray} \tag{ 6 }$$

using the fact that $\sum i {{\theta }_{iz}}=\sum j {{\theta }_{jz}}$.

BKN is an expectation–maximization (EM) algorithm, which can be used to find the maximum likelihood in an iterative manner. The log likelihood increases monotonically by updating the parameters iteratively using equations (4) and (5) or equations (4) and (6), which can be viewed as two variants of EM—the former is known as the expectation conditional maximization (ECM) [29], and the latter as an incremental variant of the EM [30]. Multiple random initializations are required to escape from local maxima and the division with the highest log likelihood is chosen as the result.

Because the physical meaning of ${{q}_{ij}}(z)$ is the weight of the edge between vertices i and j belonging to module z, this quantity is used to infer link communities in the network. Mathematically, vertex i belongs to module z if $\sum j {{A}_{ij}}{{q}_{ij}}(z)\geqslant 1$ ³ , and the physical meaning of the left-hand side of this inequality is the average number of links belonging to module z that link to vertex i.

In this subsection, we redefine BKN in the context of a bipartite network. In a unipartite network, all links exist within the same group of vertices. Consequently, one set of parameters, θ, is enough to describe the underlying structure of a unipartite network. On the other hand, there are two sets of vertices in a bipartite network; thus, two sets of parameters should be used for the parameterization. In a bipartite network ${{G}_{B}}=(U,V,E)$, U and V denote two disjointed vertex sets and E denotes the edge set. We use two sets of parameters ${{\theta }^{(U)}}$ and ${{\theta }^{(V)}}$ to describe the underlying structure of a bipartite network, such that $\theta _{iz}^{(U)}$ represents the propensity of vertex i in the vertices of set U belonging to module z, and $\theta _{jz}^{(V)}$ represents the propensity of vertex j in the vertices of set V belonging to module z. Similar to BKN, the physical meaning of $\theta _{iz}^{(U)}\theta _{jz}^{(V)}$ is the expected number of links between vertex i in vertex set U and vertex j in vertex set V belonging to module z, with the exact number being Poisson distributed. By summing over all K modules, $\sum z \theta _{iz}^{(U)}\theta _{jz}^{(V)}$ represents the expected number of links between vertex i in vertex set U and vertex j in vertex set V, the exact number of which remains Poisson distributed. The likelihood function of generating a bipartite network ${{G}_{B}}$ is

$$\begin{eqnarray}&&P\left( {{G}_{B}}|{{\theta }^{(U)}},{{\theta }^{(V)}} \right)=\mathop \prod \limits_{ij} \frac{{{\left( \mathop \sum \limits_{z} \theta _{iz}^{(U)}\theta _{jz}^{(V)} \right)}^{{{A}_{ij}}}}}{{{A}_{ij}}!}{\rm exp} \left( -\mathop \sum \limits_{z} \theta _{iz}^{(U)}\theta _{jz}^{(V)} \right).\end{eqnarray} \tag{ 7 }$$

By taking the logarithm of the likelihood function, and applying Jensen's inequality with an arbitrary variable ${{q}_{ij}}(z)$ that satisfies $\sum z {{q}_{ij}}(z)=1$, we obtain

$$\begin{eqnarray}&&{\rm ln} P\left( {{G}_{B}}|{{\theta }^{(U)}},{{\theta }^{(V)}} \right){\mkern 1mu} {\mkern 1mu} {\mkern 1mu} \geqslant \mathop \sum \limits_{ijz} \left[ {{A}_{ij}}{{q}_{ij}}(z){\rm ln} \left( \frac{\theta _{iz}^{(U)}\theta _{jz}^{(V)}}{{{q}_{ij}}(z)} \right)-\theta _{iz}^{(U)}\theta _{jz}^{(V)} \right],\end{eqnarray} \tag{ 8 }$$

where the optimal value of ${{q}_{ij}}(z)$ is given by

$$\begin{eqnarray}&&{{q}_{ij}}(z)=\frac{\theta _{iz}^{(U)}\theta _{jz}^{(V)}}{\mathop \sum \limits_{z} \theta _{iz}^{(U)}\theta _{jz}^{(V)}},\end{eqnarray} \tag{ 9 }$$

with the physical meaning unchanged compared to that of BKN. The optimal values of $\theta _{iz}^{(U)}$ and $\theta _{jz}^{(V)}$ are given by

$$\begin{eqnarray}&&\theta _{iz}^{(U)}=\frac{\mathop \sum \limits_{j} {{A}_{ij}}{{q}_{ij}}(z)}{\mathop \sum \limits_{j} \theta _{jz}^{(V)}},\ \;\ \;\ \;\theta _{jz}^{(V)}=\frac{\mathop \sum \limits_{{\rm i}} {{A}_{ij}}{{q}_{ij}}(z)}{\mathop \sum \limits_{i} \theta _{iz}^{(U)}},\end{eqnarray} \tag{ 10 }$$

by differentiating equation (8) with respect to $\theta _{iz}^{(U)}$ and $\theta _{jz}^{(V)}$. In statistics, this is known as ECM [29]. The criterion that vertex i in U belongs to in module z is $\sum j {{A}_{ij}}{{q}_{ij}}(z)\geqslant 1$ and that of vertex j in V in module z is $\sum i {{A}_{ij}}{{q}_{ij}}(z)\geqslant 1$.

As noted, the mixture network can be considered as a general case, with the unipartite and bipartite networks as limiting cases. Consequently, we can define the generative model in a unified framework.

Given a graph $G=(U,V,E)$, where U and V represent two sets of vertices and E is the set of edges between vertices in U and V, the relationship between U and V determines the type of G, which can classify networks into three types. First is the unipartite network: when U and V are identical, all edges are within the same set of vertices, then G is unipartite. Second is the bipartite network: when U and V are disjointed, G is bipartite and each edge in E links nodes from two disjointed sets. Third is the mixture network: when U and V are neither identical nor disjointed, i.e., they have some common vertices but not all vertices are the same, it is a mixture network.

Vertex labels, or identifiers, are required to determine whether a vertex is shared by two vertex sets. Let $l_{i}^{(U)}$ denote the label of vertex i in U, and $l_{j}^{(V)}$ denote the label of vertex j in V. In each set, any label is unique. When the label for vertex i in U is the same as that of vertex j in V, i and j are the same vertex rather than two different vertices. For simplicity, let ${{O}_{U}}$ and ${{O}_{V}}$ denote two sets of nodes in U and V that share common labels (common vertices), while S_U and ${{S}_{V}}$ denote two sets of vertices in U and V that do not share any common labels (specific vertices).

An illustration of the relationship between the mixture network, the unipartite network, and the bipartite network is given in figure 1.

Zoom In Zoom Out Reset image size

The generalization of the bipartite network version model to that of the mixture network is based on a simple idea: if two vertices are actually the same, they should have only one underlying tendency to belong to each module, which is reached by constraining each pair of tendencies of these two nodes to each module to be the same. This idea is quite similar to that in the comment of Costa and Hensen [31], which suggests that a unipartite network can be transformed into a bipartite network by adding explicit constraints that specify that each pair of corresponding nodes must belong to the same community. Again, the tendencies of vertex i in U and vertex j in V belonging to module z are denoted as $\theta _{iz}^{(U)}$ and $\theta _{jz}^{(V)}$. Consequently, the likelihood function of a mixture network is

$$\begin{eqnarray}&&P\left( G|{{\theta }^{(U)}},{{\theta }^{(V)}} \right)=\mathop \prod \limits_{ij} \frac{{{\left( \mathop \sum \limits_{z} \theta _{iz}^{(U)}\theta _{jz}^{(V)} \right)}^{{{A}_{ij}}}}}{{{A}_{ij}}!}{\rm exp} \left( -\mathop \sum \limits_{z} \theta _{iz}^{(U)}\theta _{jz}^{(V)} \right),\end{eqnarray} \tag{ 11 }$$

subject to constraints that

$$\begin{eqnarray}&&{{\delta }_{l_{i}^{(U)},l_{j}^{(V)}}}\left( \theta _{iz}^{(U)}-\theta _{jz}^{(V)} \right)=0\ \;\ \;\ \;\ \;\ \;\forall i,j,z,\end{eqnarray} \tag{ 12 }$$

where δ is the Kronecker delta function. Fundamentally speaking, a random network generated by this generative model is a directed network, because its unipartite part does not necessarily need to be symmetric⁴ . For an undirected network, we adopt the view that an undirected edge can be viewed as two directed edges in opposite directions [32]. For self-loops, a slight difference exists compared to those in BKN. Unlike other undirected edges, a self-loop only has one corresponding element in the adjacency matrix. Thus, when vertex i in U and vertex j in V are actually one node and they have a link (i.e., self-loop), the value of ${{A}_{ij}}$ is 1 rather than 2, regardless of whether the network is directed or undirected.

Such an optimization problem with equality constraints can be solved using the Lagrange multiplier. (The proof is available in appendix A.)

We introduce an arbitrary variable ${{q}_{ij}}(z)$ that satisfies $\sum z {{q}_{ij}}(z)=1$. The optimal value of ${{q}_{ij}}(z)$ is given by

$$\begin{eqnarray}&&{{q}_{ij}}(z)=\frac{\theta _{iz}^{(U)}\theta _{jz}^{(V)}}{\mathop \sum \limits_{z} \theta _{iz}^{(U)}\theta _{jz}^{(V)}}.\end{eqnarray} \tag{ 13 }$$

For $\theta _{iz}^{(U)}$ and $\theta _{jz}^{(V)}$, two conditions will be considered separately:

(1) For any pair of common vertices ${{i}^{\prime}}$ in ${{O}_{U}}$ and ${{j}^{\prime}}$ in ${{O}_{V}}$ that ${{\delta }_{l_{{{i}^{\prime}}}^{(U)},l_{{{j}^{\prime}}}^{(V)}}}=1$, we obtain

$$\begin{eqnarray}\begin{array}{lllllllllllllll} \theta _{{{i}^{\prime}}z}^{(U)} & = & \frac{\mathop \sum \limits_{i} {{A}_{i{{j}^{\prime}}}}{{q}_{i{{j}^{\prime}}}}(z)+\mathop \sum \limits_{j} {{A}_{{{i}^{\prime}}j}}{{q}_{{{i}^{\prime}}j}}(z)}{\mathop \sum \limits_{i} \theta _{iz}^{(U)}+\mathop \sum \limits_{j} \theta _{jz}^{(V)}}, \\ \theta _{{{j}^{\prime}}z}^{(V)} & = & \frac{\mathop \sum \limits_{i} {{A}_{i{{j}^{\prime}}}}{{q}_{i{{j}^{\prime}}}}(z)+\mathop \sum \limits_{j} {{A}_{{{i}^{\prime}}j}}{{q}_{{{i}^{\prime}}j}}(z)}{\mathop \sum \limits_{i} \theta _{iz}^{(U)}+\mathop \sum \limits_{j} \theta _{jz}^{(V)}}; \\ \end{array}\end{eqnarray} \tag{ 14 }$$

(2) For the specific vertex i in ${{S}_{U}}$ or j in S_V, we obtain

$$\begin{eqnarray}&&\theta _{iz}^{(U)}=\frac{\mathop \sum \limits_{j} {{A}_{ij}}{{q}_{ij}}(z)}{\mathop \sum \limits_{j} \theta _{jz}^{(V)}},\ \;\ \;\ \;\theta _{jz}^{(V)}=\frac{\mathop \sum \limits_{i} {{A}_{ij}}{{q}_{ij}}(z)}{\mathop \sum \limits_{i} \theta _{iz}^{(U)}}.\end{eqnarray} \tag{ 15 }$$

Our model for a mixture network is a general solution, because both BKN and the bipartite network version described in section 2.2 can be described as limiting cases. When all vertices of U and V are common vertices (unipartite network), only equation (14) remains. Because the network that BKN studied is symmetric, we have $\sum {\rm i} {{A}_{i{{j}^{\prime}}}}{{q}_{i{{j}^{\prime}}}}(z)=\sum j {{A}_{{{i}^{\prime}}j}}{{q}_{{{i}^{\prime}}j}}(z)$ and $\sum i \theta _{iz}^{(U)}=\sum j \theta _{jz}^{(V)}$, and equation (14) can be reduced to equation (5). When all vertices of U and V are specific vertices (bipartite network), only equation (15) remains, which is the same as equation (10). Conversely, the form of equations for ${{q}_{ij}}(z)$ remains constant for all three versions of the models.

In the general model, the criteria for a node to belong to a module depend on whether this node is a common or specific vertex and whether the network is directed or undirected. For a specific vertex, the criterion is same as that of the bipartite network version model: vertex i in ${{S}_{U}}$ belongs to module z if $\sum j {{A}_{ij}}{{q}_{ij}}(z)\geqslant 1$, and vertex j in ${{S}_{V}}$ belongs to module z if $\sum i {{A}_{ij}}{{q}_{ij}}(z)\geqslant 1$. For a common vertex, the criterion depends on whether the network is directed or undirected. When the network is undirected, the unipartite part of the network is symmetric. For a pair of corresponding common vertices $i'$ and $j'$, we only need to count the links between the common vertices once, i.e., $\sum i\in {{{\rm S}}_{U}} {{A}_{i{{j}^{\prime}}}}{{q}_{i{{j}^{\prime}}}}(z)+\sum j {{A}_{{{i}^{\prime}}j}}{{q}_{{{i}^{\prime}}j}}(z)\geqslant 1$. When the network is directed, links between the common vertices should be counted twice rather than once, i.e., $\sum i {{A}_{i{{j}^{\prime}}}}{{q}_{i{{j}^{\prime}}}}(z)+\sum j {{A}_{{{i}^{\prime}}j}}{{q}_{{{i}^{\prime}}j}}(z)\geqslant 1$. According to the criteria, a node can belong to one module, several modules, or none of the modules.

Using the generative model described in section 2 to generate random networks can be viewed as the reverse process of finding network modules using the same model. When finding modules, we infer module structures from the adjacency matrix. When the model is used to generate random networks, community structures are known and adjacency matrices are generated based on the model's hidden variables.

The process of generating a random mixture network is described as follows: (1) assign module membership to vertices; (2) set expected number of links for nodes to each module; (3) calculate parameters; (4) calculate expected number of links between each pair of nodes; and (5) generate the adjacency matrix. In this study, we generate adjacency matrices with binary elements by forcing their elements to be 1 if they are larger than 1, despite that our model can deal with multi-edge problem technically.

Let $k_{iz}^{(U)}=\sum j {{A}_{ij}}{{q}_{ij}}(z)$ and $k_{jz}^{(V)}=\sum i {{A}_{ij}}{{q}_{ij}}(z)$, which are set by the user to denote the expected number of links belonging to module z that link to vertex i in U and vertex j in V, respectively. The values of $\sum i\in {{O}_{U}} \theta _{iz}^{(U)}$, $\sum i\in {{S}_{U}} \theta _{iz}^{(U)}$, $\sum j\in {{S}_{V}} \theta _{jz}^{(V)}$, and $\sum j\in {{O}_{V}} \theta _{jz}^{(V)}$ can be given by

$$\begin{eqnarray}&&\begin{array}{l} \sum\limits_{i \in {O_U}} {\theta _{iz}^{(U)}} = \sqrt {\frac{{{{\left( {\sum\nolimits_{i \in {O_U}} {k_{iz}^{(U)}} + \sum\nolimits_{j \in {O_V}} {k_{jz}^{(V)}} } \right)}^2} - {{\left( {\sum\nolimits_{i \in {S_U}} {k_{iz}^{(U)}} - \sum\nolimits_{j \in {S_V}} {k_{jz}^{(V)}} } \right)}^2}}}{{2\left( {\sum\nolimits_{i \in {S_U}} {k_{iz}^{(U)}} + \sum\nolimits_{j \in {S_V}} {k_{jz}^{(V)}} + \sum\nolimits_{i \in {O_U}} {k_{iz}^{(U)}} + \sum\nolimits_{j \in {O_V}} {k_{jz}^{(V)}} } \right)}},} \\ \sum\limits_{i \in {S_U}} {\theta _{iz}^{(U)}} = \frac{{2\sum\nolimits_{i \in {S_U}} {k_{iz}^{(U)}} }}{{\sum\nolimits_{j \in {S_V}} {k_{jz}^{(V)}} - \sum\nolimits_{i \in {S_U}} {k_{iz}^{(U)}} + \sum\nolimits_{i \in {O_U}} {k_{iz}^{(U)}} + \sum\nolimits_{j \in {O_V}} {k_{jz}^{(V)}} }}\sum\limits_{i \in {O_U}} {\theta _{iz}^{(U)},} \\ \sum\limits_{j \in {S_V}} {\theta _{jz}^{(V)}} = \frac{{2\sum\nolimits_{j \in {S_V}} {k_{jz}^{(V)}} }}{{\sum\nolimits_{i \in {S_U}} {k_{iz}^{(U)}} - \sum\nolimits_{j \in {S_V}} {k_{jz}^{(V)}} + \sum\nolimits_{i \in {O_U}} {k_{iz}^{(U)}} + \sum\nolimits_{j \in {O_V}} {k_{jz}^{(V)}} }}\sum\limits_{i \in {O_U}} {\theta _{iz}^{(U)},} \\ \sum\limits_{j \in {O_V}} {\theta _{jz}^{(V)}} = \sum\limits_{i \in {O_U}} {\theta _{iz}^{(U)}.} \end{array}\end{eqnarray} \tag{ 16 }$$

For common vertices i' in $O_{U}$ and j' in ${{O}_{V}}$ that ${{\delta }_{l_{{{i}^{\prime}}}^{(U)},l_{{{j}^{\prime}}}^{(V)}}}=1$,

$$\begin{eqnarray}\begin{array}{lllllllllllllll} \theta _{{{i}^{\prime}}z}^{(U)} & = & \mathop \sum \limits_{i\in {{O}_{U}}} \theta _{iz}^{(U)}\frac{k_{{{i}^{\prime}}z}^{(U)}+k_{{{j}^{\prime}}z}^{(V)}}{\mathop \sum \limits_{i\in {{O}_{U}}} k_{iz}^{(U)}+\mathop \sum \limits_{j\in {{O}_{V}}} k_{jz}^{(V)}}, \\ \theta _{{{j}^{\prime}}z}^{(V)} & = & \mathop \sum \limits_{j\in {{O}_{V}}} \theta _{jz}^{(V)}\frac{k_{{{i}^{\prime}}z}^{(U)}+k_{{{j}^{\prime}}z}^{(V)}}{\mathop \sum \limits_{i\in {{O}_{U}}} k_{iz}^{(U)}+\mathop \sum \limits_{j\in {{O}_{V}}} k_{jz}^{(V)}}. \\ \end{array}\end{eqnarray} \tag{ 17 }$$

For specific vertices i in ${{S}_{U}}$ and j in ${{S}_{V}}$,

$$\begin{eqnarray}&&\theta _{iz}^{(U)}=\mathop \sum \limits_{i\in {{S}_{U}}} \theta _{iz}^{(U)}\frac{k_{iz}^{(U)}}{\mathop \sum \limits_{i\in {{S}_{U}}} k_{iz}^{(U)}},\;\theta _{jz}^{(V)}=\mathop \sum \limits_{j\in {{S}_{V}}} \theta _{jz}^{(V)}\frac{k_{jz}^{(V)}}{\mathop \sum \limits_{j\in {{S}_{V}}} k_{jz}^{(V)}}.\end{eqnarray} \tag{ 18 }$$

The process for generating a mixture network, as described here, can also be used to generate a unipartite directed network. However, for our bipartite network version, the model identifiability problem should be addressed before the model is used to generate random bipartite networks (the solution to parameters is not unique, which is treated in detail in appendix B and discussed in section 8). This problem can be easily solved by introducing an extra constraint $\sum i \theta _{iz}^{(U)}=\sum j \theta _{jz}^{(V)}$. Consequently, we have

$$\begin{eqnarray}&&\theta _{iz}^{(U)}=\frac{k_{iz}^{(U)}}{\sqrt{\mathop \sum \limits_{i} k_{iz}^{(U)}}},\;\theta _{jz}^{(V)}=\frac{k_{jz}^{(V)}}{\sqrt{\mathop \sum \limits_{j} k_{jz}^{(V)}}}.\end{eqnarray} \tag{ 19 }$$

The details are available in appendix B.

We use three indices to measure the performance of the algorithm in this subsection: (1) the fraction of correctly classified vertices, (2) the Jaccard index for overlapping vertices (vertices belonging to more than one module), and (3) variant normalized mutual information (NMI) [33]. Because our algorithm allows one vertex to belong to multiple modules, a node is correctly classified if the predicted memberships for all modules are correct. The definition of the Jaccard index is as follows: let ${{S}_{1}}$ denote the known set of overlapping vertices and ${{S}_{2}}$ is the predicted set of overlapping vertices—the Jaccard index is defined as $|{{S}_{1}}\cap {{S}_{2}}|/|{{S}_{1}}\cup {{S}_{2}}|$. NMI measures the similarity between the actual community structures and the predictions from the perspective of information theory: it is one if the prediction is the same as the actual condition and zero when they are independent. If a generated random network contains any isolated nodes, they would be removed from the network before detecting the modules, and therefore cannot be used to measure the performance of the algorithm.

The total number of vertices is fixed to 2000 for all networks in our consistency tests. For random unipartite networks, adjacency matrices of 2000 × 2000 are generated. For random bipartite networks, adjacency matrices of 1000 × 1000 are generated so that the total number of vertices is unchanged. For mixture networks, adjacency matrices of 1500 × 1500 are generated with 1000 common vertices (50%); the total number of unique vertices is still 2000.

First, we perform three consistency tests on random networks with two modules. In the first sets of consistency testing, 55% of vertices belong to module 1 and 55% to module 2 (10% to both modules). The total expected number of links of each vertex is set to be k. That is, if a node is a specific vertex that only belongs to one module, its expected number of linkages to this module is k; on the other hand, if a specific vertex belongs to two modules, its expected number of linkages to both modules is $k/2$. For a common vertex, the number of linkages in both vertex sets is $k/2$ if it belongs to only one module and $k/4$ if it belongs to two modules.

We test the performance of the algorithm by varying k. When k is small, the actual modular structures of the networks are quite hard to find. When k is large, it is quite easy to reveal the modules of networks. As shown in the left panels of figure 2(a), the performance of the algorithm on all types of networks increases with k as expected. When $k\leqslant 3$, the fractions of correctly classified vertices are far away from zero, but NMIs are quite close to zero, indicating that predictions for networks with low degrees are random with respect to the actual community structures. The performance at $k=4$ increases rapidly. When $k\geqslant 10$, the algorithm can determine the modular structures for all three types of networks with high performance.

Zoom In Zoom Out Reset image size

In the second set of consistency testing, k is fixed at 10, the fraction of overlapping vertices remains at 10%, and the fraction of vertices belonging to module 1 increases. The condition when the fraction of vertices belonging to module 1 is smaller than the fraction belonging to module 2 does not need to be discussed because it is equivalent to the condition considered here. As suggested by the middle panels of figure 2(a), all indices decrease with further unbalanced fractions of two modules.

In the third set of testing, k is fixed at 10, the fractions of modules 1 and 2 are set to be equal, and the fraction of overlapping vertices increases. The right panels of figure 2(a) indicate that all indices decrease with an increase of the fraction of overlapping vertices $\left( \leqslant 0.6 \right)$. When the fraction of overlapping vertices is large, the fractions of correctly classified vertices and the Jaccard indices increase, but the NMIs are close to zero, indicating that the predictions for networks with large fractions of overlapping vertices are also random with respect to the actual community structures.

We also perform three consistency tests on random networks of more than two modules. In these tests, the total number of vertices is still fixed to 2000, the fraction of overlapping nodes is set to 10%, and the fractions of nodes specifically belonging to each module are equal but k varies. For the overlapping nodes, the number of modules they belong to is restricted by an integer interval given by the user. Once the interval is given, we randomly generate the number of modules for each overlapping node and randomly assign it to a number of modules.

First, we set the number of modules to be five and each overlapping node belongs to two modules. Because the fraction of correctly classified vertices is difficult to calculate for random networks with more than two modules, and the Jaccard index for overlapping vertices cannot reflect whether they are assigned to the correct modules, these two indices are not calculated in the following tests. As shown in the left panel of figure 2(b), NMIs increase with increasing k. Similar results are observed for two sets of random networks with 10 modules, the first in which each overlapping vertex belongs to two modules (figure 2(b), middle panel) and the second of which belongs to 2–10 modules (figure 2(b), right panel). For the latter, the curves for all three types of networks increase more slowly than for the former.

In all sets of testing, the curves of the indices for all three types of networks overlap with each other very well at each point (figure 2). This suggests that our algorithm performs equally for all three types of networks generated by the same parameters.

In summary, the consistency test results suggest that our algorithm can perform well for a large range of parameters for all three types of networks. All indices increase with increasing k. For random networks with two modules, all indices approach one when the difference between the fractions of two modules or the fraction of overlapping vertices is small.

Next, we test the performance of the algorithm based on synthetic networks generated by sampling, the pseudocode of which is given in figure 3. In this process, r random unipartite networks are generated for each measured parameter of network quality. And for each random unipartite network, we generate s networks by sampling. Note that in sampling, some information is lost, because many elements in the adjacency matrix of the unipartite network become unknown in sampled networks.

Zoom In Zoom Out Reset image size

Random unipartite networks used in this section are generated using two well-known benchmarks that are of nonoverlapping community: the Girvan–Newman (GN) benchmark [2, 34] and the Lancichinetti–Fortunato–Radicchi (LFR) benchmark [35–37]. Using algorithms that allow an overlapping community to detect a nonoverlapping community is straightforward by assigning each node to its most probable module. For the GN benchmark, there are 128 vertices and 4 modules in each generated random network, with an expected degree of 16 for each node. Compared to the GN benchmark, it is more difficult to examine community structure in random networks generated by the LFR benchmark. For the LFR benchmark, we generate random unipartite networks with 1000 nodes, with community size ranging from 20 to 100 (labeled 'B'), and all other parameters are the same as those used in [37]. Both directed and undirected random networks of the two benchmarks are generated using the software described in [35–37]. The difficulty of finding network modules is controlled by the mixing parameter μ, which determines the ratio of the external degree linking to other modules to the total degree of a node. When μ increases, the community structure is more mixed. The community structure vanishes when μ is large enough.

We define two sampling rules: symmetric and asymmetric sampling. For symmetric sampling, the fractions of specific vertices in U and V are the same (denoted as u and v), with the fraction of common vertices being changed (denoted as c). For both benchmarks, we change c from 0 to 1, which determines the type of synthetic network generated by sampling: (1) $c=0$, meaning the sampled networks are bipartite; (2) $c=1$, meaning the sampled networks are unipartite and identical to the original unipartite network; and (3) $0\lt c\lt 1$, meaning the sampled networks are mixture networks. For asymmetric sampling, we fix c and then change u and v accordingly. These three parameters are not independent; their sum is always one.

First, we apply our algorithm to a set of symmetric sampling undirected networks generated from the GN benchmark. As shown in figure 4 (left upper panel), when μ increases, NMIs decrease for random networks with any c, as expected. When $\mu \leqslant 0.3$, there is no large difference between the NMIs of networks with different c. When $\mu =0.4$, the performance for our model on networks with small c begins to decrease. The smaller c is, the worse the performance. It is not hard to understand, because the fraction of unknown edges increases with the decrease of c. When $\mu \geqslant 0.5$, the performance vanishes. The tendency for the performance of symmetric sampling for directed networks is generally similar to those of undirected ones, despite the difference in the point of decrease (figure 4, left lower panel).

Zoom In Zoom Out Reset image size

For the GN benchmark, we also apply the algorithm to four sets of asymmetric sampling networks, c of which are fixed to 0.2 and 0.8 (two tests for undirected networks and two for directed networks). For the former tests, the performance for networks of $u=0.5$ decreases the slowest, because they have more remaining known edges than those of other fractions (figure 4, middle panels). For the latter tests, performances for all u are quite similar because $c=0.8$, and no matter how u and v change, the information loss is quite similar (figure 4, right panels).

For the LFR benchmark, we also apply our algorithm to two sets of symmetric sampling networks and four sets of asymmetric sampling networks. Similar results are obtained from tests of this more difficult benchmark. For the symmetric sampling undirected networks (figure 5, left upper panel), when $\mu \lt =0.5$, the algorithm easily finds the modules, despite that the performance for networks with c close to 0 begins to decrease when μ is close to 0.5. When $\mu =0.6$, the performance decreases. The smaller c is, the worse the performance. When $\mu \geqslant 0.7$, the algorithm fails. For the remaining five tests, the performance tendency is generally similar to that of the GN benchmark, despite the difference in the point of decrease (figure 5, left bottom, middle, and right panels). The similar tendencies between the performance of our algorithm on the easier and more difficult tests suggest that our model is applicable not only to small and easy networks, but also to large and difficult ones.

Zoom In Zoom Out Reset image size

We apply our algorithm to a well-known real-world bipartite network: the southern women dataset [45], which consists of 18 women (W1–W18) and 14 social events (E1–E14). This dataset has been studied extensively to investigate the results of various methods [46].

To compare the result to those of other algorithms, we first display the result of our algorithm by assigning nodes to two nonoverlapping modules. The nonoverlapping partition divides the 18 women into W1–W9 and W10–W18, which agrees with the 'perfect' partition as suggested in [46] (figure 8(a)). Unlike many algorithms mentioned in [46], our algorithm assigns social events to modules along with women. For the social events, E1–E8 belong to the same module as W1–W9 and E9–E14 as W10–W18.

Zoom In Zoom Out Reset image size

When allowing overlapping modules, W8, W9, and E6–E9 belong to both modules (figure 8(b)). This is quite reasonable, because all links only exist between those overlapping nodes and the remaining specific nodes of two modules. These overlapping nodes maintain the connection between these two modules.

Next, we apply our algorithm to a recently released human transcriptional regulatory network⁵ obtained from [15], which is a mixture directed network. This dataset consists of 82 TFs and 3998 downstream target genes, which include 31 TFs. The number of modules is set at 30. One hundred random initializations are used for the algorithm and the run with the highest log likelihood is chosen as the result. For 30 modules, the average number of TFs is 5.3 and the average number of target genes is 242.2.

Modular structures identified in the network can find genes that are known to be related to one another or offer new biological insights. The underlying assumption is that genes in the same module tend to be involved in the same or similar biological pathways. Gene ontology (GO) [48] is a standard tool used to associate biological functions with gene sets. Here, we perform GO enrichment analysis [49] to determine the enriched GO terms for downstream target genes. We report two case studies on the modules found in this transcriptional regulatory network. Supplementary figure 1 displays a module consisting of 3 TFs and 469 downstream target genes (including 2 TFs). The enriched GO terms are mainly related to the cell cycle and the metabolic process (supplementary table 1). There are three TFs here: HDAC2, NFKB1, and NFYA. All three genes are known to be related to the cell cycle [50–53] and NFYA is also related to the metabolic process [54].

A second case study focuses on a module related to the immune system (supplementary figure 2). It consists of 130 downstream genes (including 2 TFs) and 10 TFs, 6 of which (IRF1, POLR3A, POU2F2, PRDM1, FOSL2, and IRF4) are related to the immune system [54–59]. For downstream genes, the enriched GO terms are shown in supplementary table 2, many of which are immune related. We predict that another 4 TFs (PBX3, POU5F1, PPARGC1A, and MAFK) may also be related to the immune system.