%% Simulating a network
% Let's take a look at simulation.

%% Load the model
% First we load the model using the |loadModel| command.  Using the command
% in this format

m = loadModel('NFkB');

%% Defining an experiment
% In order to simulate the model, we must now let the simulator know how
% long to run a simulation, and how the input (if any) should vary during
% simulation.  This is accomplished by calling the function
% |constructExperiment| with the final time and input values.
% 
% In this case, the model measures time in seconds, but we want to
% simulate for a few hours.  

tF = 3600*5;  % simulate for 5 hrs
u = [1e-5     % TNFa quantity
        1     % IkBa+/+
        0     % IkBb-/-
        0];   % IkBe-/-

expt = constructExperiment(m, tF, u);

%% Simulate the model
% To simulate the model, we have a simple function |simulate|.  If called
% without output arguments, this function simulates the model and plots the
% output species.

simulate(m, expt);
% input('Press any key to continue.')


%% Manually plotting output
% Simulate can also be called with various output arguments, in which case
% it does not automatically plot the output.  We can do that ourselves
[t y] = simulate(m, expt);
plot(t, y(t))

%%
% What if we want to plot the third species, even if it isn't an output?
% We can access the full state trajectories x(t), and evaluate x(t) for any
% sets of species.
[t y x] = simulate(m, expt);
plot(t, x(t, [5 6]))
