% combined models
%% Repressilator Model Solver.
%% Call Equation funciton of Repressilator rep.m
%% 6 variables, 10 Parameters
%% Created by Kai, 2011/05/25
%% modified by Bruce,  2011/06/08

clc;
clear all;

dt = 0.1 ;
t = 0:dt:150 ;% Simulation time

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% opamp model
%set up output array for 3 voltages
v = zeros(length(t),3) ;
v(1,1:3) = [ 0.1 0 0];

% model parameters
tau = 2; %Time constant in Sec of RC
%k = tau/dt ;% constant from differential equation
g = -3*[1 1 1] ; % gain of opamp inverter
vMax = 10;
vMin = 0;
vThreshold = 1.5;
%solve by Euler iteration
for i=2:length(t)
    % output of opamps
    va = g .* ([v(i-1,3) v(i-1,1) v(i-1,2)] - vThreshold) ;
    va(va > vMax) = vMax;
    va(va < vMin) = vMin;
    % change of voltage on capacitor
    v(i,1:3) = v(i-1,1:3)+ dt/tau*(va - v(i-1,1:3));
    %v(i,1:3) = (va - k*v(i-1,1:3))/(1+k) ;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Full Hill model
%set up output array for 6 state variables
y = zeros(length(t),6) ;
y(1,1:6) = [0 0.1 0 0 0 0 ] ;% Initial Values
%y(:,1) = mRNA   :LacI
%y(:,2) = Protein:LacI
%y(:,3) = mRNA   :tetR
%y(:,4) = Protein:tetR
%y(:,5) = mRNA   :cI
%y(:,6) = Protein:cI

% reaction paramenters
alpha = 15*[1 1 1]; %% Regulated Growth Rate 15
beta  = 1*[1 1 1]; % Ratio of Protein Decay Rate to mRNA Decay Rate 4
alpha0= 0.1*[ 1 1 1]; % Basal Growth Rate 0.012
n = 2;% Hill coefficent

%solve by Euler iteration
for i=2:length(t)
   
%     dy(1) = -y(1) + alpha(1)/(1+y(6)^n) + alpha0(1); %y(1) = mRNA    LacI
%     dy(2) = -beta(1)*(y(2)-y(1));                    %y(2) = Protein LacI
%     dy(3) = -y(3) + alpha(2)/(1+y(2)^n) + alpha0(2); %y(3) = mRNA    tetR
%     dy(4) = -beta(2)*(y(4)-y(3));                    %y(4) = Protein tetR
%     dy(5) = -y(5) + alpha(3)/(1+y(4)^n) + alpha0(3); %y(5) = mRNA    cI
%     dy(6) = -beta(3)*(y(6)-y(5));                    %y(6) = Protein cI
    
    y(i,1) = y(i-1,1) + dt*(-y(i-1,1) + alpha(1)/(1+y(i-1,6)^n) + alpha0(1));
    y(i,3) = y(i-1,3) + dt*(-y(i-1,3) + alpha(2)/(1+y(i-1,2)^n) + alpha0(2));
    y(i,5) = y(i-1,5) + dt*(-y(i-1,5) + alpha(3)/(1+y(i-1,4)^n) + alpha0(3));
    y(i,2) = y(i-1,2) + dt*(-beta(1)*(y(i-1,2)-y(i-1,1))) ;
    y(i,4) = y(i-1,4) + dt*(-beta(2)*(y(i-1,4)-y(i-1,3))) ;
    y(i,6) = y(i-1,6) + dt*(-beta(3)*(y(i-1,6)-y(i-1,5))) ;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1); clf;
subplot(2,1,1)
plot(t,v(:,1),t,v(:,2),t,v(:,3));
legend('v1','v2','v3');
title('opamp Euler')

subplot(2,1,2)
plot(t,y(:,2), t,y(:,4), t,y(:,6));
legend('LacI','tetR','cI');
title('Hill full model Euler')

figure(2); clf;
plot(t,v(:,1),t,v(:,2),t,v(:,3));
hold on
plot(t,y(:,2), t,y(:,4), t,y(:,6), 'linewidth', 3);
title('Both -- heavy lines Hill -- light lines opamp')

figure(3); clf
plot3(v(:,1),v(:,2),v(:,3));
hold on
plot3(y(:,2),y(:,4),y(:,6),'r');

xlabel('x1')
ylabel('x2')
zlabel('x3')
box on
%set(gca, 'Projection', 'perspective')
%set(gcf, 'Renderer', 'opengl')
axis vis3d


