%% Erythrocyte Metabolism
%
% This example shows how to analyze various features of a metabolic network.

% Copyright 2015 Domenico Gatti.


%% Aims
%
% * Learn about the stoichiometries in the network.
%
% * Calculate fluxes, sensitivities, and control coefficients.
% 
% * Determine the effects of changing the boundary conditions.
%

%% Background

%% Load the Project
% Begin by loading the project at the command-line using |sbioloadproject|.

sbioloadproject Erythrocyte_Metabolism_Practice

% Alternatively we can directly import the SBML file:

% m1 = sbmlimport('Erythrocyte_Metabolism_Practice.xml');

%% Add some external programs
% The calculations and the plots will require some external programs:
% will add them to our path.

addpath(genpath(Ô../../GENERAL_SCRIPTS_FUNCTIONS'));
addpath(genpath(Ô../../TOOLBOXES/MetaTool'));
addpath(genpath(Ô../../TOOLBOXES/FastICA_25'));
addpath(genpath(Ô../../TOOLBOXES/RADICAL_1_2'));

%%
% The project contains a model, |m1|, of the wild-type erythrocyte
% metabolic network displayed below. (You can explore the network
% interactively by starting the SimBiology desktop with |simbiology| and
% opening the project file Erythrocyte_Metabolism_2.sbproj.

%%
% <<../Erythrocyte_Metabolism_Practice.png>>

%% 

m1
m1.Reactions
m1.Species
% m1.Parameters

n_species = size(m1.Species,1);
n_reactions = size(m1.Reactions,1);

%% Reorder reactions
% To simplify some of the following operations we will reorder the reactions 
% in such a way that all those internal to the erythrocyte appear first
% followed by all the exchanges with the plasma.
%
% for i = 1:n_reactions
%     if regexp(m1.Reactions(i).name,'exchange')
%         ind(i) = 1;
%     elseif regexp(m1.Reactions(i).name,'transport')
%         ind(i) = 1;
%     end
% end
% 
% internal_ind = find(~ind);
% exchange_ind = find(ind);
% 
% reorder(m1, m1.Reactions([internal_ind exchange_ind]));
% 
% sbmlexport(m1,'Erythrocyte_Metabolism_Reordered.xml');

%% Compute the stoichiometric and adjacency matrices.

% We retrieve the stoichiometric matrix for the erythrocyte model by
% passing the model object as an input argument to the 'getstoichmatrix'
% function. 

[S,objSpecies,objReactions] = getstoichmatrix(m1);
spy(S)

objSpecies 
objReactions 

% The stoichiometric matrix is typically stored as a sparse matrix. We
% convert it to a full matrix to see more easily the relationships between
% species and reactions:

SM = full(S);

% Notice that in SimBiology the stoichiometric matrix contains only
% unidirectional reactions; thus if we want a representation of all
% reactions we need to duplicate the columns corresponding to a reversible
% reaction.

rev_reactions = cell2mat(get(m1.Reactions,'Reversible'));
irrev_reactions = + (rev_reactions == 0);

rSM = -SM(:,rev_reactions);
aSM = [SM rSM];
spy(aSM);

% An alternative view of a metabolic network can be obtained by taking the
% negative transpose of the stoichiometric matrix. In this case the rows
% represent reactions and the columns are metabolites. Metabolites that
% enter a reaction are shown with a positive sign, those that leave the
% reaction with a negative sign.

RSM = -SM';

% The fundamental topological features of the metabolic network are
% determined by the non-zero elements of SM. A binary form of SM can be
% constructed by replacing any non-zero element (positive or negative) of
% SM with a 1:

SM_bin = +(SM ~= 0);

% Post-multiplication of this binary matrix by its transpose leads to a
% m x m symmetric matrix, the 'compound adjacency matrix', cAM:

cAM = SM_bin*SM_bin';

% The diagonal elements of cAM are the number of reactions in which each
% compound (species) participates (the 'reaction participation number').
% The off-diagonal elements of cAM show the number of reactions in which
% two compounds participate.

% Pre-multiplication of this binary matrix by its transpose leads to a
% symmetric matrix, the 'reaction adjacency matrix', rAM:

rAM = SM_bin'*SM_bin;

% The diagonal elements of rAM are the number of compounds (species) that
% participate in each reaction (the 'compound participation number'). The
% off-diagonal elements of rAM count how many compounds are shared by two
% reactions.

% In addition to the 'compound' and 'reaction' adjacency matrix the more
% general term of 'adjacency matrix' is reserved for a n x n matrix, where
% n equals the total number of species + reactions in a model. Each row
% corresponds to a species or reaction, and each column corresponds to a
% species or reaction. The matrix indicates which species and reactions are
% involved as reactants and products:
% 
% Reactants are represented in the matrix with a 1 at the appropriate
% location (row of species, column of reaction). Reactants appear above the
% diagonal. Products are represented in the matrix with a 1 at the
% appropriate location (row of reaction, column of species). Products
% appear below the diagonal. All other locations in the matrix contain a 0.
% For example, if a model object contains one reaction equal to A + B -> C
% and the Name property of the reaction is R1, the adjacency matrix is:
% 
%        A    B    C   R1
%   A    0    0    0   1
%   B    0    0    0   1
%   C    0    0    0   0
%   R1   0    0    1   0

% Get the adjacency matrix for m1:

[AM, Headings] = getadjacencymatrix(m1);
Headings 

% Convert the adjacency matrix from sparse to full:

AM = full(AM);

%% SVD analysis of the stoichiometric matrix.

% We calculate the 'economy size' svd. This means that if SM is m x n with
% m >= n, svd computes only the first n columns of U and S is n x n. For m
% < n, only the first m columns of V are computed and S is m-by-m.

[U,S,V] = svd(SM,'econ');

% Here we plot a spectral analysis of the stoichiometric matrix.

figure;plot(diag(S))
vline([5,32] , {'g-',});
box('on');grid off

% Create labels
xlabel('EigenReaction');
ylabel('EigenReaction Ranking');

% Percentage contribution
cumS = cumsum(diag(S));
figure;plot(cumS/sum(diag(S)))

%%
% Here we square all the values in U so we can rank the species in each
% eigenreaction by magnitude.

sU = U.^2;

% Here we look at all the species that either appear or disappear in terms
% of descending contributions.

[sorted_sU,ind_sU] = sort(sU,'descend');
% objSpecies(ind_sU(:,1))

% After sorting the elements of the U matrix in descending order we create
% a second matrix that conserves the sign of the elements of the original U
% matrix. This will be useful in plotting the eigenreactions and
% eigenconnectivity with the correct original sign.

sign_U = sign(U);
s_sorted_sU = sorted_sU.*sign_U(ind_sU);

sign_cell = ones(n_species,n_reactions);
sign_cell = num2str(sign_cell);
sign_cell(sign_U == 1) = '+';
sign_cell(sign_U == -1) = '-';
sign_cell = sign_cell(:,1:38);
% sign_cell = sign_cell(ind_sU);

%%
% We can do a similar treatment also for the eigenconnectivities.
sV = V.^2;

[sorted_sV,ind_sV] = sort(sV,'descend');
% objReactions(ind_sU(:,1))

sign_V = sign(V);
s_sorted_sV = sorted_sV.*sign_V(ind_sV);

sign_cell_c = ones(n_reactions,n_reactions);
sign_cell_c = num2str(sign_cell_c);
sign_cell_c(sign_V == 1) = '+';
sign_cell_c(sign_V == -1) = '-';
sign_cell_c = sign_cell_c(:,1:38);

%% Plot the eigenmodes
% We will plot the top 10 species in terms of contribution, in the top 4
% 'EigenReactions' and 'EigenConnectivities of the erythrocyte metabolism.

% EigenReaction 1
FIG_1 = figure;
    set(FIG_1,'Units','normalized','Position',[0.5 0.55 0.5 0.45],...
        'Name','EigenMode 1');clf

axes1 = axes('Parent',FIG_1,'YTick',[1 2 3 4 5 6 7 8 9 10],...
    'YDir','reverse','Position',[0.075 0.1 0.4 0.83]);
hold(axes1,'all');

scatter(log(sorted_sU(10:-1:1,1)),(10:-1:1),30,'MarkerEdgeColor',[1,0,0],'MarkerFaceColor',[1,0,0]); % red
scatter(log(sorted_sU(10:-1:1,2)),(10:-1:1),30,'MarkerEdgeColor',[0,1,0]); % green
scatter(log(sorted_sU(10:-1:1,3)),(10:-1:1),30,'MarkerEdgeColor',[0,0,1]); % blue
scatter(log(sorted_sU(10:-1:1,4)),(10:-1:1),30,'MarkerEdgeColor',[0,1,1]); % cyan
legend('EigenReaction 1','EigenReaction 2','EigenReaction 3',...
    'EigenReaction 4','Location','NorthWest');
box(axes1,'on');

s_sign = sign_cell(ind_sU(10:-1:1,1));
string = objSpecies(ind_sU(10:-1:1,1));
for i = 1:10
string{i,1} = [s_sign(i,1) '  ' string{i,1}(1,10:end)];
end
lscatter(log(sorted_sU(10:-1:1,1)),[9.6:-1:0.6],string,...
    'TextColor','k','FontSize',10,'FontWeight','bold');
set(gca,'XLim',[-6 -0.],'YLim',[0 11]);

% Create labels
xlabel('Natural Log of Species Contribution to the EigenReaction');
ylabel('Ranking of Metabolic Species by Contribution');

% EigenConnectivity 1
axes2 = axes('Parent',FIG_1,'YTick',[1 2 3 4 5 6 7 8 9 10],'YDir',...
    'reverse','Position',[0.55 0.1 0.4 0.83]);
hold(axes2,'all');

scatter(log(sorted_sV(10:-1:1,1)),(10:-1:1),30,'MarkerEdgeColor',[1,0,0],'MarkerFaceColor',[1,0,0]); % red
scatter(log(sorted_sV(10:-1:1,2)),(10:-1:1),30,'MarkerEdgeColor',[0,1,0]); % green
scatter(log(sorted_sV(10:-1:1,3)),(10:-1:1),30,'MarkerEdgeColor',[0,0,1]); % blue
scatter(log(sorted_sV(10:-1:1,4)),(10:-1:1),30,'MarkerEdgeColor',[0,1,1]); % cyan
legend('EigenConnectivity 1','EigenConnectivity 2','EigenConnectivity 3',...
    'EigenConnectivity 4');
set(legend,'Location','NorthWest');
box(axes2,'on')

s_sign_c = sign_cell_c(ind_sV(10:-1:1,1));
string = objReactions(ind_sV(10:-1:1,1));
for i = 1:10
string{i,1} = [s_sign_c(i,1) '  ' string{i,1}];
end
lscatter(log(sorted_sV(10:-1:1,1)),[9.6:-1:0.6],string,...
    'TextColor','k','FontSize',10,'FontWeight','bold');
set(gca,'XLim',[-6 -0.5],'YLim',[0 11]);

% Create labels
xlabel('Natural Log of Reaction Contribution to the EigenConnectivity');
ylabel('Ranking of Reactions by Contribution');

%%
% EigenReaction 2
FIG_2 = figure;
    set(FIG_2,'Units','normalized','Position',[0.5 0.55 0.5 0.45],...
        'Name','EigenMode 1');clf

axes1 = axes('Parent',FIG_2,'YTick',[1 2 3 4 5 6 7 8 9 10],...
    'YDir','reverse','Position',[0.075 0.1 0.4 0.83]);
hold(axes1,'all');

scatter(log(sorted_sU(10:-1:1,1)),(10:-1:1),30,'MarkerEdgeColor',[1,0,0]); % red
scatter(log(sorted_sU(10:-1:1,2)),(10:-1:1),30,'MarkerEdgeColor',[0,1,0],'MarkerFaceColor',[0,1,0]); % green
scatter(log(sorted_sU(10:-1:1,3)),(10:-1:1),30,'MarkerEdgeColor',[0,0,1]); % blue
scatter(log(sorted_sU(10:-1:1,4)),(10:-1:1),30,'MarkerEdgeColor',[0,1,1]); % cyan

legend('EigenReaction 1','EigenReaction 2','EigenReaction 3',...
    'EigenReaction 4');
set(legend,'Location','NorthWest');
box(axes1,'on');

s_sign = sign_cell(ind_sU(10:-1:1,2));
string = objSpecies(ind_sU(10:-1:1,2));
for i = 1:10
string{i,1} = [s_sign(i,1) '  ' string{i,1}(1,10:end)];
end
lscatter(log(sorted_sU(10:-1:1,2)),[9.6:-1:0.6],string,...
    'TextColor','k','FontSize',10,'FontWeight','bold');
set(gca,'XLim',[-6.5 0],'YLim',[0 11]);

% Create labels
xlabel('Natural Log of Species Contribution to the EigenReaction');
ylabel('Ranking of Metabolic Species by Contribution');

% EigenConnectivity 2
axes2 = axes('Parent',FIG_2,'YTick',[1 2 3 4 5 6 7 8 9 10],'YDir',...
    'reverse','Position',[0.55 0.1 0.4 0.83]);
hold(axes2,'all');

scatter(log(sorted_sV(10:-1:1,1)),(10:-1:1),30,'MarkerEdgeColor',[1,0,0]); % red
scatter(log(sorted_sV(10:-1:1,2)),(10:-1:1),30,'MarkerEdgeColor',[0,1,0],'MarkerFaceColor',[0,1,0]); % green
scatter(log(sorted_sV(10:-1:1,3)),(10:-1:1),30,'MarkerEdgeColor',[0,0,1]); % blue
scatter(log(sorted_sV(10:-1:1,4)),(10:-1:1),30,'MarkerEdgeColor',[0,1,1]); % cyan

legend('EigenConnectivity 1','EigenConnectivity 2','EigenConnectivity 3',...
    'EigenConnectivity 4');
set(legend,'Location','NorthWest');
box(axes2,'on')

s_sign_c = sign_cell_c(ind_sV(10:-1:1,2));
string = objReactions(ind_sV(10:-1:1,2));
for i = 1:10
string{i,1} = [s_sign_c(i,1) '  ' string{i,1}];
end
lscatter(log(sorted_sV(10:-1:1,2)),[9.6:-1:0.6],string,...
    'TextColor','k','FontSize',10,'FontWeight','bold');
set(gca,'XLim',[-6.5 0.5],'YLim',[0 11]);

% Create labels
xlabel('Natural Log of Reaction Contribution to the EigenConnectivity');
ylabel('Ranking of Reactions by Contribution');

%%
% EigenReaction 3
FIG_3 = figure;
    set(FIG_3,'Units','normalized','Position',[0.5 0.55 0.5 0.45],...
        'Name','EigenMode 1');clf

axes1 = axes('Parent',FIG_3,'YTick',[1 2 3 4 5 6 7 8 9 10],...
    'YDir','reverse','Position',[0.075 0.1 0.4 0.83]);
hold(axes1,'all');

scatter(log(sorted_sU(10:-1:1,1)),(10:-1:1),30,'MarkerEdgeColor',[1,0,0]); % red
scatter(log(sorted_sU(10:-1:1,2)),(10:-1:1),30,'MarkerEdgeColor',[0,1,0]); % green
scatter(log(sorted_sU(10:-1:1,3)),(10:-1:1),30,'MarkerEdgeColor',[0,0,1],'MarkerFaceColor',[0,0,1]); % blue
scatter(log(sorted_sU(10:-1:1,4)),(10:-1:1),30,'MarkerEdgeColor',[0,1,1]); % cyan
legend('EigenReaction 1','EigenReaction 2','EigenReaction 3',...
    'EigenReaction 4');
set(legend,'Location','NorthWest');
box(axes1,'on');

s_sign = sign_cell(ind_sU(10:-1:1,3));
string = objSpecies(ind_sU(10:-1:1,3));
for i = 1:10
string{i,1} = [s_sign(i,1) '  ' string{i,1}(1,10:end)];
end
lscatter(log(sorted_sU(10:-1:1,3)),[9.6:-1:0.6],string,...
    'TextColor','k','FontSize',10,'FontWeight','bold');
set(gca,'XLim',[-5.5 0.5],'YLim',[0 11]);

% Create labels
xlabel('Natural Log of Species Contribution to the EigenReaction');
ylabel('Ranking of Metabolic Species by Contribution');

% EigenConnectivity 3
axes2 = axes('Parent',FIG_3,'YTick',[1 2 3 4 5 6 7 8 9 10],'YDir',...
    'reverse','Position',[0.55 0.1 0.4 0.83]);
hold(axes2,'all');

scatter(log(sorted_sV(10:-1:1,1)),(10:-1:1),30,'MarkerEdgeColor',[1,0,0]); % red
scatter(log(sorted_sV(10:-1:1,2)),(10:-1:1),30,'MarkerEdgeColor',[0,1,0]); % green
scatter(log(sorted_sV(10:-1:1,3)),(10:-1:1),30,'MarkerEdgeColor',[0,0,1],'MarkerFaceColor',[0,0,1]); % blue
scatter(log(sorted_sV(10:-1:1,4)),(10:-1:1),30,'MarkerEdgeColor',[0,1,1]); % cyan
legend('EigenConnectivity 1','EigenConnectivity 2','EigenConnectivity 3',...
    'EigenConnectivity 4');
set(legend,'Location','NorthWest');
box(axes2,'on')

s_sign_c = sign_cell_c(ind_sV(10:-1:1,3));
string = objReactions(ind_sV(10:-1:1,3));
for i = 1:10
string{i,1} = [s_sign_c(i,1) '  ' string{i,1}];
end
lscatter(log(sorted_sV(10:-1:1,3)),[9.6:-1:0.6],string,...
    'TextColor','k','FontSize',10,'FontWeight','bold');
set(gca,'XLim',[-6.5 0.5],'YLim',[0 11]);

% Create labels
xlabel('Natural Log of Reaction Contribution to the EigenConnectivity');
ylabel('Ranking of Reactions by Contribution');

%%
% EigenReaction 4
FIG_4 = figure;
    set(FIG_4,'Units','normalized','Position',[0.5 0.55 0.5 0.45],...
        'Name','EigenMode 1');clf

axes1 = axes('Parent',FIG_4,'YTick',[1 2 3 4 5 6 7 8 9 10],...
    'YDir','reverse','Position',[0.075 0.1 0.4 0.83]);
hold(axes1,'all');

scatter(log(sorted_sU(10:-1:1,1)),(10:-1:1),30,'MarkerEdgeColor',[1,0,0]); % red
scatter(log(sorted_sU(10:-1:1,2)),(10:-1:1),30,'MarkerEdgeColor',[0,1,0]); % green
scatter(log(sorted_sU(10:-1:1,3)),(10:-1:1),30,'MarkerEdgeColor',[0,0,1]); % blue
scatter(log(sorted_sU(10:-1:1,4)),(10:-1:1),30,'MarkerEdgeColor',[0,1,1],'MarkerFaceColor',[0,1,1]); % cyan
legend('EigenReaction 1','EigenReaction 2','EigenReaction 3',...
    'EigenReaction 4');
set(legend,'Location','NorthWest');
box(axes1,'on');

s_sign = sign_cell(ind_sU(10:-1:1,4));
string = objSpecies(ind_sU(10:-1:1,4));
for i = 1:10
string{i,1} = [s_sign(i,1) '  ' string{i,1}(1,10:end)];
end
lscatter(log(sorted_sU(10:-1:1,4)),[9.6:-1:0.6],string,...
    'TextColor','k','FontSize',10,'FontWeight','bold');
set(gca,'XLim',[-5.5 0],'YLim',[0 11]);

% Create labels
xlabel('Natural Log of Species Contribution to the EigenReaction');
ylabel('Ranking of Metabolic Species by Contribution');

% EigenConnectivity 4
axes2 = axes('Parent',FIG_4,'YTick',[1 2 3 4 5 6 7 8 9 10],'YDir',...
    'reverse','Position',[0.55 0.1 0.4 0.83]);
hold(axes2,'all');

scatter(log(sorted_sV(10:-1:1,1)),(10:-1:1),30,'MarkerEdgeColor',[1,0,0]); % red
scatter(log(sorted_sV(10:-1:1,2)),(10:-1:1),30,'MarkerEdgeColor',[0,1,0]); % green
scatter(log(sorted_sV(10:-1:1,3)),(10:-1:1),30,'MarkerEdgeColor',[0,0,1]); % blue
scatter(log(sorted_sV(10:-1:1,4)),(10:-1:1),30,'MarkerEdgeColor',[0,1,1],'MarkerFaceColor',[0,1,1]); % cyan
legend('EigenConnectivity 1','EigenConnectivity 2','EigenConnectivity 3',...
    'EigenConnectivity 4');
set(legend,'Location','NorthWest');
box(axes2,'on')

s_sign_c = sign_cell_c(ind_sV(10:-1:1,4));
string = objReactions(ind_sV(10:-1:1,4));
for i = 1:10
string{i,1} = [s_sign_c(i,1) '  ' string{i,1}];
end
lscatter(log(sorted_sV(10:-1:1,4)),[9.6:-1:0.6],string,...
    'TextColor','k','FontSize',10,'FontWeight','bold');
set(gca,'XLim',[-6.5 0],'YLim',[0 11]);

% Create labels
xlabel('Natural Log of Reaction Contribution to the EigenConnectivity');
ylabel('Ranking of Reactions by Contribution');

%% 
% % These are all the species that disappear in the 1st EigenReaction
% sU_minus_ind = U<0;
% objSpecies(sU_minus_ind(:,1))
% 
% % These are the species that appear in the 1st EigenReaction
% sU_plus_ind = U>0;
% objSpecies(sU_plus_ind(:,1))
% 
% % These are all the species that disappear in the 1st EigenReaction ranked
% % according to their contribution
% sU_minus = sU;
% sU_minus(sU_plus_ind) = NaN;
% [sorted_sU_minus,ind_sU_minus] = sort(sU_minus,'descend');
% sU_minus_finite = isfinite(sorted_sU_minus);
% objSpecies(ind_sU_minus(sU_minus_finite(:,1),1))
% figure; plot(sU(ind_sU_minus(sU_minus_finite(:,1),1)))
% 
% % These are all the species that appear in the 1st EigenReaction ranked
% % according to their contribution
% sU_plus = sU;
% sU_plus(sU_minus_ind) = NaN;
% [sorted_sU_plus,ind_sU_plus] = sort(sU_plus,'descend');
% sU_plus_finite = isfinite(sorted_sU_plus);
% objSpecies(ind_sU_plus(sU_plus_finite(:,1),1))
% figure; plot(sU(ind_sU_plus(sU_plus_finite(:,1),1)))
% 
% % We can do a similar treatment also for the eigenconnectivities.
% sV = V.^2;
% [sorted_sV,ind_sV] = sort(sV,'descend');
% 
% objReactions(ind_sV(:,1))
% objReactions(ind_sV(:,2))
% objReactions(ind_sV(:,3))

%% Calculate Elementary Modes

% In order to use MetaTool we derive a vector in which every reversible
% reaction is represented by a 0 and each irreversible reaction is
% represented by a 1.

rev_reactions = zeros(1,n_reactions);
for i = 1:n_reactions
    rev_reactions(i) = get(m1.Reactions(i),'Reversible');
end
irrev_reactions = + (rev_reactions == 0);

% We will restrict the analysis to only the internal metabolites (these are
% the first 41 metabolites):

net_i.st = SM(1:41,:);
net_i.irrev_react = irrev_reactions; 

% The following line is flagged out because we have to make a keyboard
% input to execute the command.

net_i = metatool(net_i);

% We can quit here by pressing "q" or continue by pressing "enter" on the 
% keyboard.

% The 'net_i' structure contains many fields, the most important of which 
% are:
%
% st 	stoichiometric matrix (rows correspond to internal metabolites, 
%       columns to reactions)
% irrev_react 	row vector which contains 0 for a reversible and 1 for an 
%       irreversible reaction
% kn 	kernel (nullspace) of the stoichiometric matrix
% sub 	subset matrix (rows correspond to the subsets of reactions, columns  
%       to the reactions in st)
% rd 	reduced system stoichiometric matrix (rows are metabolites, columns
%       are subsets of reactions)
% irrev_rd 	reversibility of the subsets in the reduced system
% rd_ems 	elementary modes of the reduced system (rows correspond to the 
%       subsets, columns are elementary modes)
% irrev_ems 	row vector which contains 0 for a reversible and 1 for an 
%       irreversible elementary mode
% ext 	same structure as st, but rows correspond to external metabolites
% int_met 	names of the internal metabolites
% ext_met 	names of the external metabolites
% react_name 	names of the reactions

% MetaTool found 10 subsets of reactions and 28 elementary flux modes
% whose combinations lead to steady-state. We can find out which reactions
% are in each enzyme subsets:

net_i_subset_1 = m1.Reactions(find(net_i.sub(1,:)))
net_i_subset_7 = m1.Reactions(find(net_i.sub(7,:)))
net_i_subset_8 = m1.Reactions(find(net_i.sub(8,:)))

% Finally multiplying the transpose of the subset matrix * the matrix of
% the reduced elementary modes, we get the contribution of each reaction
% (the rows) to each elementary mode (the columns).

net_i.ems_part = net_i.sub' * net_i.rd_ems;

% After this expansion we can ask which reactions are involved in each
% elementary mode:

net_i_EM_1 = m1.Reactions(find(net_i.ems_part(:,1)));
net_i_EM_2 = m1.Reactions(find(net_i.ems_part(:,2)));
net_i_EM_3 = m1.Reactions(find(net_i.ems_part(:,3)));
net_i_EM_4 = m1.Reactions(find(net_i.ems_part(:,4)));
net_i_EM_26 = m1.Reactions(find(net_i.ems_part(:,26)));

% for all the modes:
% 
% n_EM = size(net_i.ems_part,2);
% net_i_EM_ind = false(n_EM,n_reactions);
% for i = 1:n_EM
%     net_i_EM_ind(i,:) = net_i.ems_part(:,i) ~= 0;
% end

% Therefore we could ask equivalently for the reactions in mode 15 or 20:

m1.Reactions(find(net_i.ems_part(:,15)))
m1.Reactions(find(net_i.ems_part(:,20)))

% or just for the enzymes involved in mode 15:

objReactions(find(net_i.ems_part(:,15)))
objReactions(find(net_i.ems_part(:,20)))

% We recall here that the columns of 'net_i.ems_part' are non-ortogonal
% vectors in the space containing all the combination of fluxes that would
% result in a steady state ( all concentration derivatives = 0 ).
% Therefore, although non-orthogonal to each other the columns of
% 'net.ems_part' should be orthogonal to the stoichiometric matrix
% 'net.st'. We can check that:

check = net_i.st*net_i.ems_part;

% We also notice that there is a very large number of possible modes.
% However, the effective dimensionality of the space is much smaller:

[Uem,Sem,Vem] = svd(net_i.ems_part);
plot(diag(Sem))

% or simply

rank(net_i.ems_part)

% eigen_net_i.ems_part = Uem(:,1:4)*Sem(1:4,1:4)*Vem(:,1:4)';
% check = net_i.st * eigen_net_i.ems_part;
% sum(check(:))

% Therefore, the first 4 columns of Uem represent an orthogonal basis for
% the space containing all the combination of fluxes that would result in a
% steady state (they are infact the 'column space' of net_i.ems_part). It's
% important to recall that there is not a UNIQUE basis for this space. For
% example, an alternative orthogonal basis would be provided by simply
% calculating the null space of 'net_i.st'.

alt_null_ortho_basis = null(net_i.st);

[~,eigenEM_ind] = sort(abs(Uem),'descend');
objReactions(eigenEM_ind(:,1))
objReactions(eigenEM_ind(:,2)) 
objReactions(eigenEM_ind(:,3))

% Furthermore, the space is NOT composed of linearly independent vectors:

check = null(net_i.ems_part);

%--------------------------------------------------------------------------
% Alternatively we could carry out the entire analysis using MetaTool as a
% stand-alone program: to achieve this we must import directly the SBML
% file into MetaTool. For this purpose we need to add to the path the
% Toolbox necessary to import the SBML file into MetaTool:

% addpath('/usr/local/lib');
% rehash
% addpath(genpath('SBMLToolbox-4.1.0'));
% 
% RCM_Model = TranslateSBML('Erythrocyte_Metabolism_2.xml');
% RCM_MT_Model = sbmlModel2metatool(RCM_Model);
% RCM = metatool(RCM_MT_Model);
% 
% % We can find out which reactions are in each enzyme subsets:
% subset_12 = m1.Reactions(find(RCM.sub(12,:)));
% subset_26 = m1.Reactions(find(RCM.sub(26,:)));
% 
% The '.rd_ems' matrix contains the contribution of each subset of
% reactions (the rows) to each elementary mode (the columns).

% Finally multiplying the transpose of the subset matrix * the matrix of
% the reduced elementary modes, we get the contribution of each reaction
% (the rows) to each elementary mode (the columns).

% RCM.ems_part = RCM.sub' * RCM.rd_ems;
%--------------------------------------------------------------------------

%%
% For brevity we will call EM the Elementary Modes matrix 'net_i.ems_part.

EM = net_i.ems_part;

% and by analogy with what we did with SM we can calculate its binary
% version;

EM_bin = +(EM ~= 0);

% Post-multiplication of this binary matrix by its transpose leads to a
% symmetric matrix, the 'reaction participation matrix', cEM:

cEM = EM_bin*EM_bin';

% The diagonal elements of cEM are the number of modes in which each
% reaction participates (the 'reaction participation number').
% The off-diagonal elements of cEM show the number of modes in which
% two reactions participate.

[sorted_diag_cEM,sorted_diag_cEM_ind] = sort(diag(cEM),'descend');
plot(sort(diag(cEM),'descend'))

% Pre-multiplication of this binary matrix by its transpose leads to a
% symmetric matrix, the 'reaction adjacency matrix', rAM:

rEM = EM_bin'*EM_bin;

% The diagonal elements of rAM are the number of reactions that participate
% in each elementary mode.The off-diagonal elements of rAM count how many
% reactions are shared by two modes. If we take the mean of the diagonal of
% rEM we will know how many reaction on average participate to each mode
% (this is also called the 'length' of the mode, and therefore the matrix
% is also called the 'pathway length matrix'):

mean_react_part = mean(diag(rEM));

% This turns out to be about 19 reactions per elementary mode. The
% distribution of the number of reactions can also be easily visualized as
% a hystogram plot. Experience with several metabolic networks has shown
% that this distribution is very often multimodal. This is indeed the case
% also for the Erythrocyte Metabolism:

hist(diag(rEM),15)

% Additional important information can be derived from an analysis of
% covariance and correlation in the EM and EM_bin matrices. First we
% calculate:

COV_m = cov(EM);
COV_r = cov(EM');
CORR_m = corr(EM);
CORR_r = corr(EM');

COV_bm = cov(EM_bin);
COV_br = cov(EM_bin');
CORR_bm = corr(EM_bin);
CORR_br = corr(EM_bin');

% The most useful ones are CORR_r which tells us which groups of reactions
% are similarly involved ( = working together) in the different modes and
% COVr, whose diagonal tells us which reactions have the widest range of
% fluxes in the different steady state elementary modes. For example:

figure;imagesc(CORR_r)
figure;plot(diag(COV_r))

objReactions([8 9 15 19:20 36])

%% Independent pathways
% However, one problem with the set of elementary flux modes is that they
% are not a true convex basis because that they are not linearly
% independent. Instead, we can derive a real 4-dimensional convex basis
% using Independent Component Analysis (ICA). We recall here that we had
% already found a 4-dimensional orthogonal basis as the null space of
% net_i.st.

iSM = net_i.st;
NS = null(iSM,'r')

% We can use fastICA or RADICAL
[ica_IC,ica_A,ica_W] = fastica(NS');
ica_IC = ica_IC';
check = ica_IC * ica_A';

% The independent component are in ica_IC as rows; we take their transpose.
% The mixing matrix is ica_A. The separating matrix is ica_W. 

[rad_IC,rad_W] = RADICAL(NS');
rad_IC = rad_IC';


% 'check should be essentially identical to EM. We can also check that the
% basis obtained in this way is really composed of linearly independent
% vectors:

null(ica_IC)
null(rad_IC)

% should give an empty matrix of dimensions 4x0. Finally we can confirm
% that the basis is convex (non-orthogonal):

NS = null(iSM)
NS'*NS
ica_IC'*ica_IC
rad_IC'*rad_IC

% Finally we make sure that ica_IC is really a convex basis for the null
% space of SM:

check = iSM * ica_IC
check = iSM * rad_IC

% Thus the columns of the IC matrix represent indeed a convex basis for the
% space containing all the flux combinations associated with a
% steady-state. It's important to remember that these bases are identified
% up to a sign, which means that any of them could be multiplied by -1.
% Finally, the ICA algorithm starts each time from a different set of
% random number, and thus typically, the solutions identified in different
% runs of the program are also slightly different. Nothwithstanding these
% limitations, the basis set identified by ICA represents the essential
% modes of the metabolic network associated with all conditions of
% steady-state.

%% Calculation of minimal pathways

% An alternative way of obtaining a biologically meaningful convex basis is 
% to use a MILP method to minimize the number of non-zero entries in the 
% basis.  

size(iSM)
[nrows,ncols] = size(iSM);

% 'rational' null space
NS = null(iSM,'r'); 
dxdt = zeros(nrows,1);

% Lower bound can be negative if there are reversible reactions
% lb = zeros(ncols,1);                 
% lb = ones(ncols,1)*min((NS(:)));
% Upper bound
% ub = ones(ncols,1)*max((NS(:)));      

lb = ones(ncols,1)*-10000; 
ub = ones(ncols,1)*10000;      
 
c = -eye(ncols);
v_opt = zeros(ncols,ncols);
fval_opt = zeros(ncols,1);
intcon = [1:ncols];
options = optimoptions('intlinprog','Display','off');
 
for i = 1:ncols
 [v_opt(:,i),fval_opt(i),exitflag,output] = ...
     intlinprog(c(i,:),intcon,[],[],iSM,dxdt,lb,ub,options);
end

v_opt_unique = (unique(round(v_opt'),'rows'))'
size(v_opt_unique,2)
sum(v_opt_unique(:) == 0)

check = iSM*v_opt_unique

% We can further notice, that some of the v_opt_unique solution
% are almost identical

v_opt_unique_trimmed = v_opt_unique(:,[2 3 5 6 8])

% Minimal Path 4 and 5 are also very similar. We can choose only pathway 4 
% bringing the total number of minimal pathways to the dimensions of the 
% null space of iSM. 

v_opt_unique_trimmed = v_opt_unique_trimmed(:,1:4)

% We can compare the number of 0's in the null space 'rational' solution NS
% and in the MILP solution:

sum(v_opt_unique_trimmed(:) == 0)
sum(NS(:) == 0)

% Finally we check that the MILP solution represents indeed steady state
% flux vectors.
check = iSM*NS
check = iSM*v_opt_unique_trimmed

% And that both solution are convex:
NS'*NS
v_opt_unique_trimmed'*v_opt_unique_trimmed

%% Compute the Conserved Quantities 
% The Simbiology function |sbioconsmoiety| examines the structure of a
% model's stoichiometry matrix to find linear combinations of species that
% are conserved. This analysis is structural in that it relies only on the
% stoichiometry and structure of the network and not on reaction kinetics.
% In fact, all the reaction rates in this model have been set to 0 because
% these rates are not important to our analysis. Here we call
% |sbioconsmoiety| with the algorithm specification |'semipos'|, so that
% all conserved quantities returned involve only positive sums of species.
% The third argument |'p'| asks for the output to be printed to a cell
% array of strings.

cons_wt = sbioconsmoiety(m1,'semipos','p');

% The output of |sbioconsmoiety| contains quantities whose time rate of
% change is zero, regardless of reaction kinetics. The species on the
% boundary of the system are present in the experimental model, and thus
% their amounts are conserved.

% However, we can be more 'hands-on' and calculate ourselves the conserved
% pools. The left null space of the stoichiometric matrix contains
% combinations of the time derivatives that add up to zero. This summation
% adds up to zero in both dynamic and steady states of the network; hence
% it contains time-invariant quantities. This condition is also called
% "node neutrality". As with the (right) null space, the choice of basis
% for the left null space is important in describing its contents in
% biochemically and biologically meaningful terms. A convex representation
% of the left null space has proven useful.
%
% The left null space of S is defined by:
%
% LS = 0 
%
% The dimension of the left null space of the stoichiometric matrix is m-r.
% Equation LS = 0 represents a linear combination of nodes whose links sum
% to zero, as can be seen pictorially from:

% (-)(|||) = 0 

% where the vertical lines represent column vectors and the horizontal
% lines represent row vectors. Thus the rows of L are a set of linearly
% independent vectors spanning the left null space, and they are orthogonal
% to the reaction vectors comprising S:
%
% l(i) · s(j) = 0
%
% The set of lines may represent mass conservation of atomic elements
% (e.g., hydrogen, carbon, oxygen), molecular subunits (e.g., carboxyl
% group, hydroxyl group), or chemical moieties (e.g., pentose, purine). The
% elemental vectors e(i), for instance, represent conservation
% relationships for the elements. Thus, the vectors that span the left null
% space represent the conservation relationships that are orthogonal to the
% reaction vectors s(j) that span the column space.
% 
% Multiplying both sides of l(i) · s(j) = 0 by the flux vector 'v' we
% obtain:
%
% l(i)  ·  (dx/)dt = 0 
%
% which is a conservation relationship, or a summation of concentrations,
% called a "pool", that is time invariant. All the time derivatives can be
% written simultaneously as:
%
% d/dt(Lx) = 0
%
% which represents all the conservation relationships and defines the
% pools. While there can be dynamic motion taking place in the column space
% along the reaction vectors, these motions do not change the total amount
% of mass in the time-invariant pools . Note that since the basis for the
% left null space is nonunique, there are many ways to represent these
% metabolic pools. The time derivatives can be positive or negative.
% Equation d/dt(Lx) = 0 can be integrated to give the mass conservation
% equations:
%
% Lx = a 
%
% where a is a vector that gives the sizes (the total concentration) of the
% pools.
% 
% The L matrix (left-null space of SM) is the same as the null space of the 
% transpose of SM:

LN = null(SM');

% To obtain an intuitive representation of the pools it is better to use
% the 'rational' representation of the null space.
LN = null(SM','r');

% We can prove that it is orthogonal to SM:

check = LN'*SM

% and also to the column space of SM which is a basis for the space of all
% the time derivatives of the concentrations:
[U,S,V] = svd(SM);

rankSM = rank(SM);
check = U(:,rank(SM)+1:end)'*U(:,1:rank(SM))
check = LN'*U(:,1:rank(SM))

% We can check the species in the various pools
m1.species(logical(LN(:,1)))
m1.species(logical(LN(:,2)))
m1.species(logical(LN(:,3)))
m1.species(logical(LN(:,4)))
m1.species(logical(LN(:,5)))
m1.species(logical(LN(:,6)))
m1.species(logical(LN(:,7)))
m1.species(logical(LN(:,8)))

% or in any combination of these pools
m1.species(logical(LN(:,8)+(LN(:,1:7)*ones(7,1))))  
% 
%% Dynamic Simulation of the Erythrocyte Metabolism
% We will simulate the model for 600s and store the results. To accomplish
% this we change the StopTime of the default configuration set object from
% to 600s.

 configsetObj = getconfigset(m1);
 set(configsetObj, 'StopTime', 600);
 set(configsetObj.SolverOptions, 'AbsoluteTolerance', 1.e-9);
 
% configsetObj.RuntimeOptions.StatesToLog = ...
%      sbioselect(m1, 'Type', 'species', 'Where', 'Name', '~=', 'example');

% Now we run a small simulation with the existing parameters.

RBC_metabolism = sbiosimulate(m1);
sbioplot(RBC_metabolism)

Plasma = sbioselect(m1,'Name','plasma');
Plasma_Glucose = sbioselect(m1,'Name','Plasma Glucose');
Plasma_Phosphate = sbioselect(m1,'Name','Plasma Phosphate');
Plasma_Pyruvate = sbioselect(m1,'Name','Plasma Pyruvate');
Plasma_Lactate = sbioselect(m1,'Name','Plasma Lactate');

Plasma_Init = Plasma;
Plasma_Glucose_Init = Plasma_Glucose;
Plasma_Phosphate_Init = Plasma_Phosphate;
Plasma_Pyruvate_Init = Plasma_Pyruvate;
Plasma_Lactate_Init = Plasma_Lactate;

%%
% In normal subject the ratio between Plasma Volume and Interstitial Fluid
% is PV/IF = 0.222
% Considering an hematochrit of 45% as being 98% provided by RBC we get:
% IF = PV/0.222
% PV = BV*0.55 
% RBCV = BV*0.45 -> BV = RBCV/0.45
% PV = (RBCV/0.45)*0.55
% IF = (RBCV/0.45)*0.55 * (1/0.222)
% nonRBCV = IF + PV = ((RBCV/0.45)*0.55) * (1/0.222) +  ((RBCV/0.45)*0.55)
% nonRBCV = IF + PV = ((RBCV/0.45)*0.55) * ((1/0.222) +1)
% nonRBCV = IF + PV = ((RBCV/0.45)*0.55) * 1.222/0.222
% nonRBCV = RBCV* (0.55*1.222)/(0.45*0.222) = RBCV * ~6.73
% RBCV + OCV = 2*IF = 2*[(RBCV/0.45)*0.55 * (1/0.222))]
% OCV = 2*[(RBCV/0.45)*0.55 * (1/0.222))] - RBCV
% OCV = RBCV*((0.55*1)/(0.45*0.222)) - RBCV = RBCV*(5.5055 - 1) 
% OCV = RBCV * ~4.51
% tot_nonRBCV = OCV + nonRBCV = RBCV * (6.73+4.51) = RBCV * 11.24
% Thus, the total amount of intra-cellular and extra-cellular fluid in the
% body is approximately 11.24 time the total volume of the erythrocytes.

% We can see what would happen to the amount of plasma lactate produced by 
% the RBCs if it was not used by other tissues (mostly the liver).

set(configsetObj, 'StopTime', 3600);
set(configsetObj.SolverOptions, 'AbsoluteTolerance', 1.e-7);

set(Plasma, 'Capacity', 11.24);
set(Plasma_Lactate, 'ConstantAmount', false);
set(Plasma_Lactate, 'BoundaryCondition', false);
set(Plasma_Pyruvate, 'ConstantAmount', false);
set(Plasma_Pyruvate, 'BoundaryCondition', false);

Plasma_Glucose = Plasma_Glucose_Init;

RBC_metabolism = sbiosimulate(m1);
sbioplot(RBC_metabolism)

% We can see how the concentrations of cellular fructose 1-6 biphosphate,
% and dihydroxyaceton phosphate, and of both cellular and plasma lactate
% increase dramatically.

% openvar RBC_metabolism.Time
% openvar RBC_metabolism.Data

%%
% Alternatively we can see what would happen under conditions in which the
% plasma glucose and lactate remain constant, but increase (e.g. Diabetes):

set(Plasma_Lactate, 'ConstantAmount', false);
set(Plasma_Lactate, 'BoundaryCondition', true);
set(Plasma_Pyruvate, 'ConstantAmount', false);
set(Plasma_Pyruvate, 'BoundaryCondition', true);
set(configsetObj, 'StopTime', 100);
set(configsetObj.SolverOptions, 'AbsoluteTolerance', 1.e-9);

% In this case it is better configure the solver so that states are logged
% only at specific times.

set(configsetObj.SolverOptions, 'OutputTimes',(0:1:100));

set(Plasma_Pyruvate, 'ConstantAmount', false);
set(Plasma_Pyruvate, 'BoundaryCondition', true);

set(Plasma_Glucose, 'InitialAmount', 5);
set(Plasma_Lactate, 'InitialAmount', 1.68);
set(Plasma_Pyruvate, 'InitialAmount', 0.084);

step0 = sbiosimulate(m1);

set(Plasma_Glucose, 'InitialAmount', 10);
set(Plasma_Lactate, 'InitialAmount', 3);
% set(Plasma_Pyruvate, 'InitialAmount', 0.16);

step1 = sbiosimulate(m1);

set(Plasma_Glucose, 'InitialAmount', 15);
set(Plasma_Lactate, 'InitialAmount', 5);
% set(Plasma_Pyruvate, 'InitialAmount', 0.32);

step2 = sbiosimulate(m1);

% Here we concatenate the simulations in order to visualize them in a
% single plot.
last_time = step0.Time(end);
Time_012 = [step0.Time;step1.Time(2:end) + 1*last_time;step2.Time(2:end) + 2*last_time];
Data_012_orig =[step0.Data;step1.Data(2:end,:);step2.Data(2:end,:)];
Time_Data_012_orig = [Time_012 Data_012_orig];
Data_012 = Time_Data_012_orig(:,2:42);
Data_012_zscore = zscore(Data_012);
Data_012_zscore_max = max(abs(Data_012_zscore))
[~,rank] = sort(Data_012_zscore_max,'descend');
%%
for i = 1:41
    string{i,1} = objSpecies{i,1}(1,10:end);
end

FIG_3 = figure;
    set(FIG_3,'Units','normalized','Position',[0.3 0.1 0.7 0.9],...
        'Name','');clf

axes1 = axes('Parent',FIG_3,...
    'Position',[0.075 0.53 0.4 0.45]);
hold(axes1,'all');
box(axes1,'on');
plot(Time_012,Data_012(:,1:10))
legend(string(1:10),'Location','Best')
% set(gca,'YLim',[0 1.6]);
ylabel('Species Concentration (mM)')

axes2 = axes('Parent',FIG_3,...
    'Position',[0.55 0.53 0.4 0.45]);
hold(axes2,'all');
box(axes2,'on');
plot(Time_012,Data_012(:,11:20))
legend(string(11:20),'Location','Best')
ylabel('Species Concentration (mM)')
    
axes3 = axes('Parent',FIG_3,...
    'Position',[0.075 0.05 0.4 0.45]);
hold(axes3,'all');
box(axes3,'on');
plot(Time_012,Data_012(:,21:30))
legend(string(21:30),'Location','Best')
ylabel('Species Concentration (mM)')
xlabel('Time (s)');
    
axes4 = axes('Parent',FIG_3,...
    'Position',[0.55 0.05 0.4 0.45]);
hold(axes4,'all');
box(axes4,'on');
plot(Time_012,Data_012(:,31:41),'LineStyle','-')
legend(string(31:41),'Location','Best')
ylabel('Species Concentration (mM)')
xlabel('Time (s)');

%%
s_string = string(rank);
Data_012_zscore_ranked = Data_012_zscore(:,rank);

FIG_4 = figure;
    set(FIG_4,'Units','normalized','Position',[0.3 0.1 0.7 0.9],...
        'Name','');clf

axes1 = axes('Parent',FIG_4,...
    'Position',[0.075 0.53 0.4 0.45]);
hold(axes1,'all');
box(axes1,'on');
plot(Time_012,Data_012_zscore_ranked(:,1:10))
legend(s_string(1:10),'Location','Best')
% set(gca,'YLim',[0 1.6]);
ylabel('Species Concentration (mM)')

axes2 = axes('Parent',FIG_4,...
    'Position',[0.55 0.53 0.4 0.45]);
hold(axes2,'all');
box(axes2,'on');
plot(Time_012,Data_012_zscore_ranked(:,11:20))
legend(s_string(11:20),'Location','Best')
ylabel('Species Concentration (mM)')
    
axes3 = axes('Parent',FIG_4,...
    'Position',[0.075 0.05 0.4 0.45]);
hold(axes3,'all');
box(axes3,'on');
plot(Time_012,Data_012_zscore_ranked(:,21:30))
legend(s_string(21:30),'Location','Best')
ylabel('Species Concentration (mM)')
xlabel('Time (s)');
    
axes4 = axes('Parent',FIG_4,...
    'Position',[0.55 0.05 0.4 0.45]);
hold(axes4,'all');
box(axes4,'on');
plot(Time_012,Data_012_zscore_ranked(:,31:41),'LineStyle','-')
legend(s_string(31:41),'Location','Best')
ylabel('Species Concentration (mM)')
xlabel('Time (s)');


