function dydt = StochasticSwitch(t,y)
%-----------------------------------------------------------------------
% Specify the ODE
%-----------------------------------------------------------------------

% Work in a more comprehensible set of variables
%%  Initialisation
mRNA_LacI  = y(1);
LacI= y(2);
IPTG_LacI=y(3);
mRNA_XylR=y(4);
XylR=y(5);
Xylose_XylR=y(6);
mRNA_Rfp=y(7);
Rfp=y(8);
mRNA_Gfp=y(9);
Gfp=y(10);
mRNA_Hin_RightToLeft=y(11);
Hin=y(12);
mRNA_Hin_LeftToRight=y(13);
mRNA_AraR=y(14);
AraR=y(15);
Arabinose_AraR=y(16);
mRNA_SspB=y(17);
SspB=y(18);
SspB_Hin=y(19);
%%
global IPTG;
global Xylose;
global Arabinose;
global flip;
global Hin_ToFlip;
global iteration;

%%  Parameters
gamma_mRNA = 0.0058;
gamma_protein = 0.0012;
k_ClpXP_degradation=0.2;
%Generic kf and kb parameters
kf=0.001;
kb=0.5;
kf_DNA=0.0001;


Km_LacI_pspac=100;
Km_XylR_xylA=400;
Km_AraR_araR=800;
Km_AraR_araE=100;

n=2;

k_transcription_LacI=0.07;
k_translation_LacI= 0.1;

k_transcription_XylR=0.075;
k_translation_XylR= 0.11;

k_transcription_HinLeftToRight=0.13;
k_transcription_HinRightToLeft=0.05;
k_translation_Hin= 0.2;

k_transcription_Gfp=0.1;
k_translation_Gfp= 0.2;

k_transcription_Rfp=0.1;
k_translation_Rfp= 0.2;

k_transcription_AraR=0.06;
k_translation_AraR= 0.1;

k_transcription_SspB=0.18;
k_translation_SspB= 0.26;


%%    IPTG, Xylose, Arabinose Reactions
J_IPTG_LacI_Formation=IPTG*LacI*kf-kb*IPTG_LacI;
J_Xylose_XylR_Formation=Xylose*XylR*kf-kb*Xylose_XylR;
J_Arabinose_AraR_Formation=Arabinose*AraR*kf-kb*Arabinose_AraR;
J_ClpXP_Degradation=k_ClpXP_degradation * SspB_Hin;
J_SspB_Hin_Formation=Hin*SspB*kf_DNA-kb*SspB_Hin;


%%    Dynamic equations
dydt(1)=  k_transcription_LacI * 1.4 - gamma_mRNA * mRNA_LacI; % dmRNA_LacI/dt
dydt(2)=  k_translation_LacI * mRNA_LacI - gamma_protein * LacI - J_IPTG_LacI_Formation; % dLacI/dt
dydt(3)= J_IPTG_LacI_Formation-gamma_protein*IPTG_LacI;%dIPTG_LacI/dt
dydt(4)=  k_transcription_XylR * 1.4 - gamma_mRNA * mRNA_XylR; % dmRNA_XylR/dt
dydt(5)=  k_translation_XylR * mRNA_XylR - gamma_protein * XylR - J_Xylose_XylR_Formation; % dXylR/dt
dydt(6)= J_Xylose_XylR_Formation-gamma_protein*Xylose_XylR;%dXylose_XylR/dt
dydt(7)=  0;% dmRNA_Rfp/dt
dydt(8)=  0;% dRfp/dt
dydt(9)=  0;% dmRNA_Gfp/dt
dydt(10)= 0;% dGfp/dt
dydt(11)= 0;%mRNA_Hin_RightToLeft/dt
dydt(12)= 0;%dHin/dt
dydt(13)= 0;%mRNA_Hin_LeftToRight/dt
dydt(14)=  k_transcription_AraR * 1.4 * power(Km_AraR_araR, n)/(power(Km_AraR_araR, n)+power(AraR, n))- gamma_mRNA * mRNA_AraR; % dmRNA_AraR/dt
dydt(15)=  k_translation_AraR * mRNA_AraR - gamma_protein * AraR - J_Arabinose_AraR_Formation; % dAraR/dt
dydt(16)= J_Arabinose_AraR_Formation-gamma_protein*Arabinose_AraR;%dArabinose_AraR/dt
dydt(17)=  k_transcription_SspB * 1.4 * power(Km_AraR_araE, n)/(power(Km_AraR_araE, n)+power(AraR, n))- gamma_mRNA * mRNA_SspB; % dmRNA_SspB/dt
dydt(18)=  k_translation_SspB * mRNA_SspB - gamma_protein * SspB - J_SspB_Hin_Formation + J_ClpXP_Degradation; % dSspB/dt
dydt(19)= J_SspB_Hin_Formation-gamma_protein* SspB_Hin - J_ClpXP_Degradation; %dSspB_Hin/dt



%% Flipping Decision
numberOfIterationsToPass=50;
iteration=iteration+1;
% Hin is being expressed from left to right
% If Hin level reaches to Hin_ToFlip and not flipped, flip it
% with %40 of chance
if (Hin>=Hin_ToFlip && iteration>=numberOfIterationsToPass && flip==0 && randint(1,1,[1,10])>6);
    flip=1;
    iteration=0;
    %if still not flipped try it again later when Hin level reaches to
    %the new Hin_ToFlip level. Do not flip this time
elseif (Hin>=Hin_ToFlip && iteration>=numberOfIterationsToPass && flip==0);
    iteration=0;
    % Hin is being expressed from right to left
    % If Hin level reaches to Hin_ToFlip and flipped already flip it
    % with %40 of chance
elseif (Hin>=Hin_ToFlip && iteration>=numberOfIterationsToPass && flip==1 && randint(1,1,[1,10])>6);
    flip=0;
    iteration=0;
    %if still not flipped try it again later when Hin level reaches to
    %the Hin_ToFlip level
elseif (Hin>=Hin_ToFlip && iteration>=numberOfIterationsToPass && flip==1 );
    iteration=0;
end;

%%      Hin mRNA transcription based on the flipping decision
if (flip==1);
    % It is flipped to the left. Express Hin mRNA from right to left
    %Calculate d[mRNA_Hin_RightToLeft]/dt
    mRNA_Hin_RightToLeft_Flux=k_transcription_HinRightToLeft*1.4* power(Km_XylR_xylA, n)/(power(Km_XylR_xylA, n)+power(XylR, n)) -gamma_mRNA*mRNA_Hin_RightToLeft;
    
    new_mRNA_Hin_RightToLeft=mRNA_Hin_RightToLeft + mRNA_Hin_RightToLeft_Flux;
    if (new_mRNA_Hin_RightToLeft>=0);
        dydt(11)= mRNA_Hin_RightToLeft_Flux;%dHin_mRNA_RightToLeft/dt
    else
        dydt(11)=0;
    end;

    %Calculate d[mRNA_Hin_LeftToRight]/dt. It only degrades. No expression!
    if ((mRNA_Hin_LeftToRight - gamma_mRNA*mRNA_Hin_LeftToRight)>=0);
        dydt(13)= -gamma_mRNA*mRNA_Hin_LeftToRight;
    else
        dydt(13)=0;
    end;


    %Calculate the expression of Gfp mRNA. It is flipped and it is expressed;
    dydt(9)=  k_transcription_Gfp * 1.4 * power(Km_XylR_xylA, n)/(power(Km_XylR_xylA, n)+power(XylR, n)) - gamma_mRNA * mRNA_Gfp; % dmRNA_Gfp/dt

    %Calculate mRNA_Rfp, It only degrades, no expression when flipped
    if (mRNA_Rfp - gamma_mRNA * mRNA_Rfp>0 );
        dydt(7)=- gamma_mRNA * mRNA_Rfp;
    else
        dydt(7)=0;
    end;

else
    % It is flipped to the right or not flipped at all. Express
    % Hin mRNA from left to right
    %Calculate d[mRNA_Hin_LeftToRight]/dt
    mRNA_Hin_LeftToRight_Flux=k_transcription_HinLeftToRight*1.4*power(Km_LacI_pspac, n)/(power(Km_LacI_pspac, n)+power(LacI, n))-gamma_mRNA*mRNA_Hin_LeftToRight;
    new_mRNA_Hin_LeftToRight=mRNA_Hin_LeftToRight + mRNA_Hin_LeftToRight_Flux;

    if (new_mRNA_Hin_LeftToRight>=0);
        dydt(13)= mRNA_Hin_LeftToRight_Flux;%dmRNA_Hin_LeftToRight/dt
    else
        dydt(13)=0;
    end;


    %             %Calculate d[mRNA_Hin_RightToLeft]/dt. It only degrades. No
    %             %expression!
    if ((mRNA_Hin_RightToLeft - gamma_mRNA*mRNA_Hin_RightToLeft)>=0);
        dydt(11)=-gamma_mRNA*mRNA_Hin_RightToLeft;
    else
        dydt(11)=0;
    end;
    %Calculate the expression of Rfp mRNA. It is expressed;
    dydt(7)=  k_transcription_Rfp * 1.4 * power(Km_LacI_pspac, n)/(power(Km_LacI_pspac, n)+power(LacI, n))- gamma_mRNA * mRNA_Rfp; % dmRNA_Rfp/dt

    %Calculate mRNA_Gfp, It only degrades, no expression when flipped to
    %the right
    if (mRNA_Gfp  - gamma_mRNA * mRNA_Gfp>0 );
        dydt(9)= -1 * gamma_mRNA * mRNA_Gfp;       
    else
        dydt(9)=0;
    end;


end;

%% Hin Translation
HinFlux= k_translation_Hin*(mRNA_Hin_RightToLeft + mRNA_Hin_LeftToRight)-gamma_protein*Hin - J_SspB_Hin_Formation;

newHin=Hin +  HinFlux;
if (newHin>=0)
    dydt(12)=HinFlux ;%dHin/dt
else
    dydt(12)=0;
end;

%% Rfp and Gfp Translation
dydt(8)=  k_translation_Rfp * mRNA_Rfp - gamma_protein * Rfp; % dRfp/dt
dydt(10)=  k_translation_Gfp * mRNA_Gfp - gamma_protein * Gfp; % dGfp/dt

%%
dydt = dydt';        %MATLAB wants a column vector