clear all;
clc;
close all;

global kml dml kpl dpl kb1 kb2 kmo dmo kpo dpo Km n; %add in hill parameters

kml = 1;
dml = 1;
kpl = 1;
dpl = 1;
kb1 = 0.8;
kb2 = 0.2;
kmo = 1;
dmo = 1;
kpo = 1;
dpo = 1;
Km = 0.5;
n = 2;


%start at 1 1 2 because: 1 1: Steady state of Lac, since it is constitutive
%Vary the 2 value (initial amounts of IPTG) and see what is going to happen
%to the protein output: As you can see overshoot in production which
%decreases when you stop adding in IPTG. 
%Maybe add an external IPTG term to see how production varies and
%kickstarts accordingly

%Vary: amounts of IPTG
%Dissociation constants accordingly
%Promoter leakiness? Observations why do we get basal levels of
%encapsulation?

[T,Y]= ode45('m1_circuit', [0 100], [1 1 2 0 0 0]);

figure(1);
subplot(2,3,1);plot(T,Y(:,1));
subplot(2,3,2);plot(T,Y(:,2));
subplot(2,3,3);plot(T,Y(:,3));
subplot(2,3,4);plot(T,Y(:,4));
subplot(2,3,5);plot(T,Y(:,5));
subplot(2,3,6);plot(T,Y(:,6));

i=1;

for IPTG=0:1:10
[T1,Y]= ode45('m1_circuit', [0:0.1:100], [1 1 IPTG 0 0 0]);
A(:,i)=Y(:,1);
B(:,i)=Y(:,2);
C(:,i)=Y(:,3);
D(:,i)=Y(:,4);
E(:,i)=Y(:,5);
F(:,i)=Y(:,6); 
i=i+1;
end
%order: [dMlacdt; dPlacdt; dIPTGdt; dlacIPdt; dMoutdt; dPoutdt];
figure(2);
subplot(2,3,1);plot(T1,A);title('MLac');xlabel('time');
subplot(2,3,2);plot(T1,B);title('PLac');
subplot(2,3,3);plot(T1,C);title('IPTG');
subplot(2,3,4);plot(T1,D);title('IPTG-LacI complex');
subplot(2,3,5);plot(T1,E);title('Mout');
subplot(2,3,6);plot(T1,F);title('Pout');legend('I=0','I=1','I=2','I=3','I=4','I=5','I=6','I=7','I=8','I=9','I=10');