%ODE SOLVER 
%Indicate initial values for model state variables, indicate [arabinose],
%[aTc], [theophylline], Kd_tbs, units: (M)

%THEOPHYLLINE: example for 1% arabinose, 25 ng/mL aTc, [theo] = 0 - 10 mM
X0 = [0 0 0 0 0 0 .0666 5.6251*10^-8 0 1.3*10^-9]';
X02 = [0 0 0 0 0 0 .0666 5.6251*10^-8 1*10^-3 1.3*10^-9]';
X03 = [0 0 0 0 0 0 .0666 5.6251*10^-8 2*10^-3 1.3*10^-9]';
X04 = [0 0 0 0 0 0 .0666 5.6251*10^-8 5*10^-3 1.3*10^-9]';
X05 = [0 0 0 0 0 0 .0666 5.6251*10^-8 10*10^-3 1.3*10^-9]';
%call ode45 
[t,X] = ode45('model2alt10_17', 36000, X0);
[t2,X2] = ode45('model2alt10_17', 36000, X02);
[t3,X3] = ode45('model2alt10_17', 36000, X03);
[t4,X4] = ode45('model2alt10_17', 36000, X04);
[t5,X5] = ode45('model2alt10_17', 36000, X05);
%plot GFP output vs time for conditions of [theo] = 1 mM, 2 mM, 5, mM, 10 mM
figure
plot(t2/3600, X2(:,6), t3/3600, X3(:,6), t4/3600, X4(:,6), t5/3600, X5(:,6)), ylabel('[GFP] (uM)'), xlabel('Time (hours)')
legend('1mM theophylline', '2mM theophylline', '5mM theophylline', '10mM theophylline'), title('Riboswitch Tunability [RiboTALe8, 1% arabinose, 25 ng/mL aTc]')
%calculate fold-change in fluorescence for [theophylline] = 0 - 10 mM
[pk_0mM] = findpeaks(X(:,6));
[pk_1mM] = findpeaks(X2(:,6));
[pk_2mM] = findpeaks(X3(:,6));
[pk_5mM] = findpeaks(X4(:,6));
[pk_10mM] = findpeaks(X5(:,6)); 
fold_0mM = pk_0mM / pk_10mM;
fold_1mM = pk_1mM / pk_10mM;
fold_2mM = pk_2mM / pk_10mM;
fold_5mM = pk_5mM / pk_10mM;
fold_10mM = pk_10mM / pk_10mM;
%plot expected fold-change in expression
theo_induction = [1 2 5 10];
fold_change_theo = [fold_1mM fold_2mM fold_5mM fold_10mM];
figure
plot(theo_induction, fold_change_theo), xlabel('[Theophylline] (mM)'), ylabel('Expected Relative GFP Expression ([GFP]/[GFP])')
title('Expected Effect of Theophylline on GFP Expression')


%KD_TBS: example for 1% arabinose, 25 ng/mL aTc, 1 mM theophylline, Kd = 1.3 nM or
%240 nM
X0kd = [0 0 0 0 0 0 .0666 5.6251*10^-8 1*10^-3 1.3*10^-9]';
X0kd2 = [0 0 0 0 0 0 .0666 5.6251*10^-8 1*10^-3 240*10^-9]';
%call ode45
[tkd, Xkd] = ode45('model2alt10_17', 36000, X0kd);
[tkd2, Xkd2] = ode45('model2alt10_17', 36000, X0kd2);
%plot GFP output vs time
figure
plot(tkd/3600, Xkd(:,6), tkd2/3600, Xkd2(:,6)), ylabel('[GFP] (uM)'), xlabel('Time (hours)')
legend ('RiboTALe8, Kd = 1.3 nM', 'RiboTALe1, Kd = 240 nM'), title('RiboTALe Tunability [1% arabinose, 25 ng/mL aTc, 1mM theophylline]')
%calculate fold-change in fluorescence for Kd_tbs = 1.3 nM, 240 nM
[pk_13] = findpeaks(Xkd(:,6))
[pk_240] = findpeaks(Xkd2(:,6))
fold_13 = pk_13/pk_13
fold_240 = pk_240/pk_13
%plot
kd = [1.3 240];
fold_change_kd = [fold_13, fold_240];
figure
plot(kd, fold_change_kd), xlabel('Kd_tbs (nM)'), ylabel('Relative GFP expression'), title('Expected Effect of Kd_tbs on GFP Expression')


%ATC: example for 1% arabinose, 1 mM theophylline, Kd=1.3 nM, aTc = 0 - 100
%ng/mL
X0atc = [0 0 0 0 0 0 .0666 0 1*10^-3 1.3*10^-9]';
X0atc2 = [0 0 0 0 0 0 .0666 5.6251*10^-8 1*10^-3 1.3*10^-9]';
X0atc3 = [0 0 0 0 0 0 .0666 2.25*10^-7 1*10^-3 1.3*10^-9]';
%call ode45
[tatc, Xatc] = ode45('model2alt10_17', 36000, X0atc);
[tatc2, Xatc2] = ode45('model2alt10_17', 36000, X0atc2);
[tatc3, Xatc3] = ode45('model2alt10_17', 36000, X0atc3);
%plot GFP output vs time
figure
plot(tatc/3600, Xatc(:,6), tatc2/3600, Xatc2(:,6), tatc3/3600, Xatc3(:,6)), ylabel('[GFP] (uM)'), xlabel('Time (hours)')
legend('0 ng/mL aTc', '25 ng/mL aTc', '100 ng/mL aTc'), title('Transcript Tunability [RiboTALe8, 1% arabinose, 1 mM theophylline]')
%calculate fold-change in fluorescence for aTc = 0 - 100 ng/mL
[pk_0] = findpeaks(Xatc(:,6))
[pk_25] = findpeaks(Xatc2(:,6))
[pk_100] = findpeaks(Xatc3(:,6))
fold_0 = pk_0/pk_0
fold_25 = pk_25/pk_0
fold_100 = pk_100/pk_0
%plot
atc = [0 25 100];
fold_change_atc = [fold_0 fold_25 fold_100];
figure
plot(atc, fold_change_atc), xlabel('[aTc] (ng/mL)'), ylabel('Relative GFP expression')
title('Expected Effect of [aTc] on GFP Expression')




