%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%            Sim_ODE_modelv1.m                                       %
%                                                                         %
%       Valencia_UPV iGEM 2015                            %
%
%    Switch deterministic model                                           %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% Model:
%
% SPECIES IN THE FIRST CONSTITUTIVE LEVEL
% x1= A (E-PIF6)
% dot{A} = alpha_A - dA*A
     % We denote by alphaA = k_{mA}*c_{gA}*k_{pA}/d_{mA} the A protein expression
     % net formation rate.
% x2= B  (PhyB-VP16)
% dot{B} = alphaB - dB*B - kRed*B*Red
    % When there is red light B is activated, transforming into Bstar
% x3= C  (Gal4-KDronp)
% dot{C} = alphaC - dC*C
% x4= D   (KDronpa)
% dot{D} = alphaD - dD*D - kBlue*D*Blue
% x5= Bstar   (Red light activated PhyB-VP16)
% dot{Bstar} = kRed*B*Red - dBstar*Bstar
% x6= Dstar   (Blue light activated KDronpa)
% dot{Dstar} = kBlue*D*Blue - dDstar*Dstar

% List of states:
%
% x1= A (E-PIF6)
% x2= B  (PhyB-VP16)
% x3= C  (Gal4-KDronp)
% x4= D   (KDronpa)
% x5= Bstar   (Red light activated PhyB-VP16)
% x6= Dstar   (Blue light activated NDronpa-VP16)
% x7= gE+ 
% x8= E (BD3-PIF6)
% x9) gF+ 
% x10= F (BD4-KDronpa)
% x11= gG+
% x12= G (Bxb1)
% x13= gH+
% x14= H (BD5-PIF6)
% x15= gI+
% x16= I (BD6-KDronpa)
% x17= gJ+
% x18= J (phiC31)
% x19= alpha
% x20= beta
% x21= gamma
% x22= Omega
% Inputs:
%
% u(1): Red(t)  Red light. Amplitude represents intensity
% u(2): Blue(t)  Blue light. Amplitude represents intensity
%
% Parameters:
%
%  alphaA, alphaB, alphaC, alphaD,... = k_m*c_g*k_p/d_m:  protein expression
     % net formation rate                         
%  dA,dB,dC,dD,dDstar,...: protein degradation rates 
%  kRed:  PhyB-VP16 activation rate 
%  kBlue:  KDronpa activation rate 


clear all
close all

% Number of Circuits
param.N = 1;

param.alphaA=5;  %Parameters iGEM Zurich 2014
param.alphaB=5; 
param.alphaC=5; 
param.alphaD=5; 
param.alphaE=5;
param.alphaF=5; 
param.alphaG=5; 
param.alphaH=5; 
param.alphaI=5; 
param.alphaJ=5;
param.alphaalp=5;
param.alphabe=5;
param.alphaga=5;
param.alphaOm=5;

basal=0.01;
param.alphaEb= basal ;
param.alphaFb= basal ;
param.alphaGb= basal ;
param.alphaHb= basal;
param.alphaIb= basal ;
param.alphaJb= basal ;
param.alphaalpb= basal ;
param.alphabeb= basal ;
param.alphagab= basal ;
param.alphaOmb= basal ;

param.dA=0.02; % min^{-1};
param.dB=0.2; 
param.dC=0.02; 
param.dD=0.2; 
param.dBstar=0.2;%TAGGED
param.dDstar=0.2;
param.dE=0.02; % min^{-1};
param.dF=0.02; 
param.dG=0.01; 
param.dH=0.02; 
param.dI=0.02;
param.dJ=0.01;
param.dalp=0.002; % min^{-1};
param.dbe=0.002; 
param.dga=0.002; 
param.dOm=0.002; 

param.kRed=1; % min^{-1} nM^{-1}
param.kBlue=1; % min^{-1} nM^{-1}

param.kphiC31=1; 
param.kbxb1=1;
N=100;
M=100;
param.me= M;
param.ne= N;
param.mf=M; 
param.nf= N;
param.mg= M;
param.ng= N;
param.mh=M; 
param.nh= N;
param.mi=M; 
param.ni= N;
param.mj=M; 
param.nj= N;
param.malp= M;
param.nalp= N;
param.mbe=M; 
param.nbe= N;
param.mga= M;
param.nga= N;
param.mOm=M; 
param.nOm= N;

param.t4=0; %TETRAMERIZATION IN %

% Options for the ode45/ode15s
opt = odeset('AbsTol',1e-3,'RelTol',1e-2,'MaxStep',1);

%%%%%%%%%%%%%%%%%%%%%%%%  Simulate %%%%%%%%%%%%%%%%
t_final=3100; % with input_t=(0:Ts:t_final)' (minutes)
t_recomb=280; %time needed for recombinases work
time_span = [0,t_final];
%We set the initial value for states: %genes no son cero condiciones inic.
gin=1;
ic = [0,0,0,0,0,0,gin,0,gin,0,gin,0,gin,0,gin,0,gin,0,0,0,0,0]; 

% We build the red light signal
Ts=1; % Sampling time to generate input signal. 
input_t=(0:Ts:t_final)';
red_intensity=1; % light intensity
blue_intensity=1; % light intensity

input_red =[zeros(1626,1);red_intensity*ones(1375,1); zeros(100,1)];
input_blue=[zeros(1626,1);blue_intensity*ones(1375,1); zeros(100,1)];

[t,x] = ode23s(@model_switch_v2, time_span, ic, opt,input_t,input_red,input_blue,param);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Plot the results.
%CONSTITUTIVE LEVEL
figure(1)
subplot(321)
plot(t,[x(:,1)],'b-') 
set(gcf, 'color', 'white');
title('Switching Network');
xlabel('Time (minutes)');
ylabel('A (nMolar)');
grid on
subplot(323)
plot(t,[x(:,2)],'b-') 
hold on
plot(t,[x(:,5)],'r-')
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('B,B* (nMolar)');
legend('B','B*')
grid on
subplot(325)
plot(input_t,input_red,'r-') 
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('Red light');
grid on

subplot(322)
plot(t,[x(:,3)],'b-') 
set(gcf, 'color', 'white');
title('Switching Network');
xlabel('Time (minutes)');
ylabel('C (nMolar)');
grid on
subplot(324)
plot(t,[x(:,4)],'b-') 
hold on
plot(t,[x(:,6)],'r-')
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('D,D* (nMolar)');
legend('D','D*')
grid on
subplot(326)
plot(input_t,input_blue,'b-') 
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('Blue light');
grid on

%Figure 2: Graphs  gE, E, gF, F, gG, G, gH, H, gI, I,gJ, J
figure(2)
%%TOGGLE RED
subplot(321)
plot(input_t,input_red,'r-') 
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('Red light');
grid on

subplot(323)
plot(t,[x(:,7)],'k-') 
hold on
plot(t,[x(:,8)],'r-')
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('gE+,E (nMolar)');
legend('gE+','E')
grid on

subplot(325)
plot(t,[x(:,11)],'k-') 
hold on
plot(t,[x(:,12)],'r-')
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('gG+,G (nMolar)');
legend('gG+','G')
grid on

%%TOGGLE BLUE
subplot(322)
plot(input_t,input_blue,'b-') 
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('Blue light');
grid on

subplot(324)
plot(t,[x(:,13)],'k-') 
hold on
plot(t,[x(:,14)],'b-')
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('gH+,H (nMolar)');
legend('gH+','H')
grid on

subplot(326)
plot(t,[x(:,17)],'k-') 
hold on
plot(t,[x(:,18)],'b-')
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('gJ+,J (nMolar)');
legend('gJ+','J')
grid on
% 
%Figure 3: Graphs  alpha, beta, gamma and omega.
figure(3)
subplot(321)
plot(input_t,input_red,'r-') 
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('Red light');
grid on

subplot(322)
plot(input_t,input_blue,'b-') 
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('Blue light');
grid on
subplot(323)
plot(t,[x(:,19)],'g-') 
set(gcf, 'color', 'white');
title('Switching Network');
xlabel('Time (minutes)');
ylabel('Alpha (nMolar)');
grid on
subplot(324)
plot(t,[x(:,20)],'g-') 
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('Beta (nMolar)');
grid on
subplot(325)
plot(t,[x(:,21)],'g-') 
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('Gamma (nMolar)');
grid on
subplot(326)
plot(t,[x(:,22)],'g-') 
set(gcf, 'color', 'white');
xlabel('Time (minutes)');
ylabel('Omega (nMolar)');
grid on
% 
figure(4)
plot(t,[x(:,19)],'g')
hold on
plot(t,[x(:,20)],'m')
plot(t,[x(:,21)],'y')
plot(t,[x(:,22)],'c')
hold off
legend ('alfa','beta','gamma','Omega')
set(gcf, 'color', 'white');
xlabel('Time(min)');
ylabel('[Outputs](nM)');
grid on

figure (5)
subplot(211)
plot([x(:,8)],[x(:,13)],'r')
set(gcf, 'color', 'white');
xlabel('Red toggle[···]');
ylabel('Blue Toggle genes[···]');
grid on
subplot(212)
plot([x(:,8)],[x(:,7)],'b');
set(gcf, 'color', 'white');
xlabel(' Blue toggle[···]');
ylabel(' Red toggle genes[···]');
grid on


OA=[x(:,22)]./[x(:,19)];
OB=[x(:,22)]./[x(:,20)];
OG=[x(:,22)]./[x(:,21)];
figure(9)
subplot(221);plot(t,OA,'g');xlabel('Time(min)');ylabel('Omega/alpha');grid on;set(gcf, 'color', 'white');
subplot(222);plot(t,OB,'m');xlabel('Time(min)');ylabel('Omega/beta');grid on;set(gcf, 'color', 'white');
subplot(223);plot(t,OG,'y');xlabel('Time(min)');ylabel('Omega/gamma');grid on;set(gcf, 'color', 'white');


