function [t,y] = RedLightSensor()
%% Parameters

tf= 1/60; % in seconds
%% guessed params
global K1           % response param OmpR-P,tRNA
K1=5;               % nM


%% lac Z
global n
global alpha_M
global alpha_B
global alpha_A
global gamma_M
global gamma_B

n=2;
alpha_M=0.997*tf;
alpha_B=1.66*1e-5*tf;
alpha_A=20*tf;
gamma_M=0.411*tf;
gamma_B=8.33*1e-4*tf;


%conversion rates

global muMpernM
global molecpernM
muMpernM=1e3;
molecpernM=0.60221;

%% Initial values
x0=zeros(1,10);
%#OmpR - bionumbers
x0(6)=3500/molecpernM; % molecules/cell;
%#EnvZ - bionumbers
x0(1)=100/molecpernM; % molecules/cell;




%% ODE
timeframe=[0 40000];

options=odeset('MaxStep','5');
[t,y]=ode23s(@odeeq,timeframe,x0,options);

%figure(7)
%title('Complete System')
%clf
%plot(t,y)
%legend('EnvZ','EnvZ-P','EnvZ-P.OmpR','EnvZ.OmpR-P','OmpR-P','OmpR','EnvZ.OmpR','lacZ mRNA','\beta-Gal','dye')

end


function [dxdt]=odeeq(t,x)
global K1

global n
global alpha_M
global alpha_B
global alpha_A
global gamma_M
global gamma_B



dxdt=[    OmpReq(t,x(1:7));
          % 8 lacZ mRNA
          alpha_M*Hill(x(5)/K1,n,1) - gamma_M*x(8);
          % 9 \beta-Galactosidase
          alpha_B*x(8) - gamma_B*x(9);
          % 10 result
          alpha_A*x(9);
          ];
 
end

function H=Hill(x,n,tresh)
H=(x.^n)/(tresh.^n+x.^n);
end
