function [subnet, subnet_r, subnet_m, fin] = grow_modular_subnetwork(subnet, net, expr, class, beta, subnet_r, subnet_m)
%GROW_MODULAR_SUBNETWORK Grow a subnetwork by one node.
%
% [subnet, subnet_r, subnet_m, fin] = GROW_MODULAR_SUBNETWORK(subnet, net, expr, class, beta, subnet_r, subnet_m)
%
% INPUT:
% subnet the set of nodes in the current subnetwork (can be a seed gene)
% net the weighted adjacency matrix (weights have values between 0 and 1)
% expr genes are rows and samples are columns (values for each gene are z-score normalized across samples)
% class the class label for each sample (age of the sample)
% beta the modularity coefficient
% subnet_r the class relevance of subnet
% subnet_m the modularity of subnet
%
% OUTPUT:
% subnet the set of nodes in the current subnetwork
% subnet_r the class relevance of subnet
% subnet_m the modularity of subnet
% fin = 1 if there were no nodes found that increased both subnet_r and subnet_m
fin = 0;
%Identify network neighbors of current subnetwork
[r, c] = find(net(subnet,:));
nbors = setdiff(c,subnet);
%Calculate the relevance of candidate subnetworks
num = length(nbors); len = length(subnet)+1;
candidates_expr = (1/len)*(repmat(sum(expr(subnet,:), 1), num, 1) + expr(nbors, :));
r = abs(corr(candidates_expr',class','type','Spearman'));
%Calculate modularity of candidate subnetworks
col = find(any(net([subnet nbors],:)));
cnbors = setdiff(col,subnet);
in = sum(net(nbors,subnet),2) + sum(sum(net(subnet,subnet)))/2;
out = sum(net(nbors,cnbors),2) + sum(sum(net(subnet,cnbors))) - sum(net(nbors,subnet),2);
m = in./(((in+out).^2)+1);
%Calculate score for candidate subnetworks
keepers = find((m >= subnet_m).*(r >= subnet_r));
nbors = nbors(keepers); r = r(keepers); m = m(keepers);
score = r + beta*m;
%Add a node (randomly chosen from among those that are tied for the highest score)
tied = find(score == max(score));
if ~isempty(tied)
ind = tied(randperm(length(tied))); ind = ind(1);
subnet = union(subnet, nbors(ind)); subnet_r = r(ind); subnet_m = m(ind);
else
fin = 1;
end