function dYdt = autocomplex2_001(t,Y)

%Y(1) = Ao
%Y(2) = ATP
%Y(3) = ADP
%Y(4) = Ap
%Y(5) = Tbound
%Y(6) = Dpbound_lsrR
%Y(7) = Dpbound_lsrA
%Y(8) = Dp
%Y(9) = D
%Y(10) = M
%Y(11) = mrnaR(mRNA|lsrR) 
%Y(12) = C1-cAMP-CRP
%Y(13) = mrnaA(mRNA|lsrA)
%Y(14) = cAMP-CRP
%Y(15) = Dbound_lsrA
%Y(16) = Dbound_lsrR
%Y(17) = cAMP
%Y(18) = CRP
%Y(19) = C2-cAMP-CRP
%Y(20) = A
%Y(21) = lsrR
%Y(22) = lsrA
%Y(23) = C1
%Y(24) = C2
%Y(25) = Bgal-reporter mRNA: lsrR
%Y(26) = B-gal reporter mRNA: lsrA
%Y(27) = Beta galactosidase: lsrR (will be reported in MU)
%Y(28) = Beta galactosidase: lsrA (will be reported in MU)
B = 1.225; %Transcription bias towards lsrA-side of lsr
k_lsrK = 1; %Rate constant in uM^-1*min^-1;
k_deg = 1; %Rate constant, ATP degradation; min^-1
k_degAp = 1;
k_degMRNA = 1;
k_degD = 1;
k_degM = 1;
kdeg_Bgal = 1;
k_bind = 1; %Rate constant, uM^-2*min^-1
k_regen = 1; %min^-1
k_cleave1 = 1; %min^-1
k_cleave2 = 1; %min^-1
k_degDp = 1; %min^-1
k_pol = 1; % uM^-1*min^-1
k_tcR = 2.1; %k_tc|lsrR (uM*min^-1)
k_tcA = B*k_tcR;
k_tl = 1; %min^-1
k_tlA = 1; %min^-1
k_tlM = k_tl;
k_cAMP = 1; %uM^-1*min^-1
k_c1 = 1; %uM^-1*min^-1
k_c2 = 1; %uM^-1*min^-1
%k_regenATP = 5;
nA = 1.4; nR = 1.4; %uM
Kbind_lsrR = 0.0001; %uM
Kbind_lsrA = 0.0001; %uM
Vmax_tet = 10; % uM^-1*min^-1
Vm_lsrA = 1; Vm_lsrR = 1; %uM/min
Kb_lsrA = 1; Kb_lsrR = 1; %uM
%DNAtot = 100; %Total number of lsrR operons available
%omega = 20; %Total concentration of CRP binding sites in genome

dAodt = -(((2*0.005)/110).*(t-240)).*exp(-(((t-240)/110).^2)); 
dATPdt = -k_lsrK*Y(1).*Y(2) - k_deg*Y(2)+k_deg*Y(2)+k_lsrK*Y(1).*Y(2);
dADPdt = k_lsrK*Y(1).*Y(2) - k_deg*Y(3);
dApdt = k_lsrK.*Y(1).*Y(2) - k_bind*(Y(4).^2)*Y(5) - k_degAp.*(Y(4));
dTdt = Vmax_tet.*(Y(15)).*(Y(16)) - k_bind.*Y(5).*Y(4).^2;
dDpbRdt = k_bind*Y(5).*(Y(4)).^2 - k_cleave1*Y(6);
dDpbAdt = k_bind*Y(5).*(Y(4)).^2 - k_cleave2*Y(7);
dDpdt = k_cleave1*Y(6) + k_cleave2*Y(7) - k_degDp*Y(8);
dDdt = k_pol*(Y(10).^2) - (Vm_lsrR.*Y(21).*((Y(9).^(nR))/(Kbind_lsrR^nR + (Y(9)).^nR))) - (Vm_lsrA.*Y(22).*((Y(9).^(nA))/(Kbind_lsrA^nA + (Y(9)).^nA))) ;
dMdt = k_tlM*Y(11) - k_pol*(Y(10).^2);
dmrnaRdt = k_tcR.*(Y(12)./(Y(12)+Y(23))).*Y(21) - (k_degMRNA).*Y(11);
dcAMPcrpbc1dt = k_c1*Y(14).*(Y(23)) - k_regen*Y(12);
dmrnaAdt = k_tcA.*(Y(19)./(Y(19)+Y(24))).*Y(22) - (k_degMRNA).*Y(13);
dcAMPCRPdt = k_cAMP.*Y(17).*Y(18) - Y(14).*(k_c1.*(Y(23)) + k_c2.*(Y(24)));
dDblsrAdt = Vm_lsrA.*(Y(22)).*(((Y(9)).^nA)./(Kbind_lsrA.^nA + (Y(9)).^nA)) - Vmax_tet*Y(15).*Y(16);
dDblsrRdt = Vm_lsrR.*(Y(21)).*(((Y(9)).^nR)./(Kbind_lsrR.^nR + (Y(9)).^nR)) - Vmax_tet*Y(15).*Y(16);
dcAMPdt = k_regen.*(Y(12)+Y(19)) - k_cAMP.*Y(17).*Y(18);
dCRPdt = k_regen.*(Y(12)+Y(19)) - k_cAMP.*Y(17).*Y(18);
dcAMPcrpbc2dt = k_c2.*Y(14).*(Y(24)) - k_regen*Y(19);
dAdt = k_tlA*Y(13) - k_deg*Y(20);
dlsrRdt = k_cleave1.*Y(6) - Vm_lsrR*Y(21).*(((Y(9)).^nR)./(Kbind_lsrR^nA + (Y(9)).^nR));
dlsrAdt = k_cleave2.*Y(7) - Vm_lsrA*Y(22).*(((Y(9)).^nA)./(Kbind_lsrA^nA + (Y(9)).^nA));
dC1dt = k_regen.*(Y(12)) - k_c1.*Y(14).*Y(23);
dC2dt = k_regen.*(Y(19)) - k_c2.*Y(14).*Y(24);
dmrnaXdt = k_tcR.*(Y(12)./(Y(12)+Y(23))).*Y(21) - (k_degMRNA).*Y(25);
dmrnaZdt = k_tcA.*(Y(19)./(Y(19)+Y(24))).*Y(22) - (k_degMRNA).*Y(26);
dBgalRdt = k_tl.*Y(25);
dBgalAdt = k_tlA.*Y(26);

dYdt = [dAodt; dATPdt; dADPdt; dApdt; dTdt; dDpbRdt; dDpbAdt; dDpdt; ...
    dDdt; dMdt; dmrnaRdt; dcAMPcrpbc1dt; dmrnaAdt; dcAMPCRPdt; dDblsrAdt; ...
    dDblsrRdt; dcAMPdt; dCRPdt; dcAMPcrpbc2dt; dAdt; dlsrRdt; dlsrAdt; ...
    dC1dt; dC2dt; dmrnaXdt; dmrnaZdt; dBgalRdt; dBgalAdt];
