%
% sim_geneosc.m
%
% simulate the genetic oscillator model.
%

% model parameters
a = 1.9;    % additional level of mRNA
a0 = 0;      % baseline level of mRNA
n = 4;       % number of proteins required to bind and repress
b = 100;      % decay rate of protein 

T = 50;    % number of time units to run

p = [a, a0, n, b];  % put all the parameters into one vector

ic = rand(6,1); % initial conditions for voltage and three gating variables

% run ode45, which will solve the system defined by the other file 'geneosc.m'
% input:
% [0 T] gives the timespan to run the system
% ic is the initial conditions
% p is the vector of parameters
%
% output:
% t will be the time vector
% y will be a matrix, each column is one of five variables in the system
opt = odeset('abstol',1e-6);
[t,y] = ode45(@geneosc,[0 T],ic,opt,p);

% pull out all six variables
m1 = y(:,1);     % lac mRNA
m2 = y(:,2);     % tet mRNA
m3 = y(:,3);     % cL mRNA
p1 = y(:,4);     % lac protein
p2 = y(:,5);     % tet protein
p3 = y(:,6);     % cL protein

% plot voltage as a function of time
figure(1), hold on, plot(t,p1,'linewidth',4);
figure(1), hold on, plot(t,p2,'r','linewidth',4);
figure(1), hold on, plot(t,p3,'y','linewidth',4);
figure(1), hold on, plot(t,(2/(n-2))^(1/n),'k','linewidth',2);
set(gca,'fontsize',20);
xlabel('time (min)','fontsize',30,'interpreter','latex');
ylabel('proteins','fontsize',30,'interpreter','latex');
