%template file for running BioNetGen simulations from Matlab

%In order to use the runbng function, you need to add to Matlab path
%the directories containing the BioNetGen installation.  Make sure you save
%the path before exiting Matlab so that you don't need to add the
%directories again in the future.

%Note, the first time you run the runbng function, a window will pop up and
%ask you to give the location of the file BNG2.pl

%Change the current directory to your RuleBender project directory

%define the same basic constants as in RuleBender
vol=1.1e-15;		%cell volume in  liters
avn=6.022e23*1e-9; %avogadros number
aPoPM=70*60/115; %average polymerase per minute per cell
aRiPM=.5*30; %average ribosomes per minute per cell
RNAd=.2; %average RNA degradation rate [per minute]

GT = 1/avn/vol;
RT = 1000;
nC = 1;

kPe = aPoPM/RNAd*aRiPM*.005;
kPd = .027;

iT = 10.^([2:.1:4]);
Pss = zeros(1,length(iT));

for i = 1:length(iT)
    changePar('NAR.bngl',{'iT'},[iT(i)]/vol/avn); %function for changing parameter values in rule bender
    runbng('','NAR.bngl');
    
    A = importdata('NAR.gdat',' ',1); %import the data
    time = A.data(:,1);
    
    %all other columns are in the order defined in RuleBender
    P = A.data(:,2);
    Pss(i) = P(end);

    figure(1), hold on
    
    h = plot(time/60,P*vol*avn);
   
    xlabel('time [hr]','fontsize',16);
    ylabel('P [#/cell]','fontsize',16);
    set(h,'linewidth',3)
    set(gca,'fontsize',16)
    axis([0 time(end)/60 0 max(P)])
   
end
hold off

figure(2)
h = semilogx(iT,Pss*vol*avn,'o');
xlabel('I [#/cell]','fontsize',16);
ylabel('P_{ss} [#/cell]','fontsize',16);
set(h,'linewidth',3)
set(gca,'fontsize',16)
%axis([iT(1) iT(end) 0 max(P)])


