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    