demo_small_network.m

％The Metabolic Networks Toolbox contains functions to create, ％modify, display, and simulate bioche

echo on;
clc
% We consider the following simple network
%
% S1             S5
%    \         /
%      S3 - S4
%    /         \
% S2             S6
%
% with Sn denoting the metabolites
% to illustrate how some of the functions in the
% Metabolic Network Toolbox can be used.
%
% Please hit any key to continue

pause
clc
% We define the network by its stoichiometric matrix N,
% the reversibility of reactions (bit string) 
% and the external metabolites (indices)
%
N            = [ -1 0 0 0 0; ...
		 0 -1 0 0 0; ...
		 1 1 -1 0 0; ...
		 0 0 1 -1 -1; ...
		 0 0 0 1 0; ...
		 0 0 0 0 1];

reversible   = [1 1 1 1 1]';
external_ind = [1 2 5 6 ]';
%
% Hit any key to continue
pause

% All information about the network is then stored in a matlab structure
network      = network_construct(N,reversible,external_ind)
% The fields 'kinetics' and 'graphics_par' contain information
% (set by default) about reaction kinetics and the graphical layout,
% respectively.

% Hit any key to continue
pause
clc
% The functions in the 'netgraph' directory change the layout
% of the network (defined by the field 'graphics_par'). 
% We can move the nodes in the graph by
%
network = netgraph_move(network);

% We can set/remove the reversible reactions ...
%
network = netgraph_reversible(network);

% ... and the external metabolites
%
network = netgraph_external(network);
clc
% Until here, we only considered the network structure.
% Now we define a kinetics.
% There are different ways to specify a kinetics (see 'kinetics_structure'). 
% Here we choose mass-action kinetics for all reactions.
%
network.kinetics = set_kinetics(network,'mass-action');
% This updates the field 'kinetics', which contains all kinetic data.
%
% Hit any key to continue
pause
clc
% Now we can run dynamical simulations.
% First, we compute a stationary state, starting from an initial guess
% of the concentrations ...

s = [1 1 0 0 0 0]';

[S, J] = network_steady_state(network,s);

% ... and plot it on the network graph. 
%

netgraph_draw(network,'metvalues',S,'arrowvalues',J,'arrowstyle','fluxes');

% The size (and color) of octagons 
% and squares represents the values of concentrations and fluxes.
%
% Hit any key to continue
pause
clc
% Likewise, we can compute properties from metabolic control analysis,
% for instance, elasticity coefficients and control coefficients
% (for more information: type 'help control_theory')
%

epsilon       = elasticities(network,S);
[ C_J, C_S ]  = control_coefficients(network.N,epsilon,network.external);

% Here we plot the control coefficients on the 4th metabolite 
%
netgraph_draw(network,'actvalues',C_S(4,:)');

% Hit any key to continue
pause
clc
% We compute a time course with initial conditions given by s:

[t,s_t,s_int] = network_integrate(network,s,5);
plot(t,s_t)

% Hit any key to continue
pause

% We can also show the time course as a movie
M = netgraph_movie(network,t,s_t);
movie(M,1);
% Hit any key to continue
pause

% Accordingly, we can compute time-dependent response coefficients

[t, s_t, RS, s_ind_t, RS_ind, parameters, parameter_names] = time_response_coefficients(network, s, 5);

plot(t,squeeze(RS(:,1,:)))
%
% Enjoy this toolbox!
echo off;