%% Model Gene Regulation Pathway
%% Create a model.

Mobj = sbiomodel('cell')

%
%% Add a compartment.

compObj = addcompartment(Mobj,'comp');
compObj.CapacityUnits = 'liter';

%
%% Add a reaction for transcription.

Robj1 = addreaction(Mobj,'DNA -> DNA + mRNA');

%
%% Display species.

Mobj.Species

%
%% Set initial amount for species DNA.

Mobj.Species(1).InitialAmount       = 50;
Mobj.Species(1).InitialAmountUnits  = 'molecule';
Mobj.Species(2).InitialAmountUnits  = 'molecule';

%
%% Add kinetic law.

Kobj1 = addkineticlaw(Robj1,'MassAction');

%
%% Display parameterVariables.

Kobj1.ParameterVariables

%
%% Display speciesVariables.

Kobj1.SpeciesVariables

%
%% Create a parameter k1.

Pobj1               = addparameter(Kobj1,'k1');
Pobj1.Value         = 0.2;
Pobj1.ValueUnits    = '1/second';
Kobj1.ParameterVariableNames = 'k1';

%
%% Display speciesVariableNames.

Kobj1.SpeciesVariableNames

%
%% Display reaction rate of Robj1.

Robj1.ReactionRate

%
%% Add a reaction for translation and set its kinetic law.

Robj2 = addreaction(Mobj,'mRNA -> mRNA + protein');
Mobj.Species(3).InitialAmountUnits = 'molecule';
Kobj2 = addkineticlaw(Robj2,'MassAction');

%
%% Define the reaction rate constant k2 for the reaction.

Pobj2               = addparameter(Kobj2,'k2');
Pobj2.Value         = 20;
Pobj2.ValueUnits    = '1/second';
Kobj2.ParameterVariableNames = 'k2';

%
%% Display the reaction rate of Robj2.

Robj2.ReactionRate

%
%% Add reactions for gene regulation.

Robj3 = addreaction(Mobj,'DNA + protein <-> DNAProteinComplex');
Mobj.Species(4).InitialAmountUnits = 'molecule';
Kobj3 = addkineticlaw(Robj3,'MassAction');
Pobj3 = addparameter(Kobj3,'k3','Value',0.2,'ValueUnits','1/(molecule*second)');
Pobj3r = addparameter(Kobj3,'k3r','Value',1.0,'ValueUnits','1/second');
Kobj3.ParameterVariableNames = {'k3','k3r'};

%
%% Display the reaction rate of Robj3.

Robj3.ReactionRate

%
%% Add reactions for mRNA degradation to nucleotides.

Robj4 = addreaction(Mobj,'mRNA -> null');
Kobj4 = addkineticlaw(Robj4, 'MassAction');
Pobj4 = addparameter(Kobj4,'k4','Value',1.5,'ValueUnits','1/second');
Kobj4.ParameterVariableNames = 'k4';

%
%% Display the reaction rate of Robj4.

Robj4.ReactionRate

%
%% Add reaction for protein degradation to amino acids.

Robj5 = addreaction(Mobj,'protein -> null');
Kobj5 = addkineticlaw(Robj5,'MassAction');
Pobj5 = addparameter(Kobj5,'k5','Value',1.0,'ValueUnits','1/second');
Kobj5.ParameterVariableNames = 'k5';

%
%% Display the reaction rate of Robj5.

Robj5.ReactionRate

%
%% Turn on unit conversion by using the configset object.

configset = getconfigset(Mobj);
configset.CompileOptions.UnitConversion = true;

%
%% Simulate the model.

[t, simdata, names] = sbiosimulate(Mobj);

%
%% Plot the results.

plot(t,simdata)
legend(names,'Location','NorthEastOutside')
title('Gene Regulation');
xlabel('Time');
ylabel('Species Amount');

%
% Copyright 2015 The MathWorks, Inc.