function r_init = Gene_Repression(Total_Repressor1,Total_Repressor2,Inducer1,Inducer2)
%Chemical partition function for repressing gene expression

%Conversion factors & constants
molec2M = 1.660577881102623e-009;
RT = 1.987e-3*310;

%Competing random DNA for RNAP/sigma -- do not change
dG_RNAP_NonSpecific = -4.5; %Estimated
dG_TF_NonSpecific = -7.5; %Estimated

%Plasmid / chromosome copy number
copynumber = 10;

%Thermodynamic Parameters [kcal/mol]
%Target Promoter
dG_Sigma70 =  -9.3 - dG_RNAP_NonSpecific; %A medium strength Sigma70 promoter
dG_SigmaS = -5 - dG_RNAP_NonSpecific; %Estimated

dG_Repressor1 = -14.5 - dG_TF_NonSpecific; %A repressor similar to the Lac repressor
dG_Repressor2 = -14.5 - dG_TF_NonSpecific; %A weaker binding repressor
dG_R12 = 0; %A protein-protein interaction between repressors 1 and 2 (for cooperative binding)

Ka_Inducer1 = 1.2e6; %M^-1
Ka_Inducer2 = 1e11; %M^-1

%Competing random DNA for RNAP/sigma -- do not change
dG_NonSpecific = -4.5; %Estimated

%All other competing promoters (average) -- do not change
dG_Sigma_OtherPromoters = -8 - dG_NonSpecific; %Estimated

%Kinetic Parameters
k_conf = 0.01 * 60; %min^-1

%Numbers of molecules of important species
Sigma70 = 4700;  %Measured experimentally
SigmaS = 100; %Measured experimentally
RNAP = 2600; %Measured experimentally

%All promoters in the genome -- do not change
Promoters = copynumber;
Other_Sigma70_Promoters = 10500; %Estimated from E. coli genome
Other_SigmaS_Promoters = 5000; %A completely arbitrary estimate
Random_DNA = 1e6; %E. coli genome size in base pairs

%Equilibrium Statements -- do not change

%1) Competition between RNAP and different sigma factors. Assumes each sigma
%factor binds to RNAP with equal affinity.
TotalSigma = Sigma70 + SigmaS;
Total_RNAP_Sigma70 = RNAP * Sigma70 / TotalSigma;
Total_RNAP_SigmaS = RNAP * SigmaS / TotalSigma;

%2) Competition between RNAP/sigma and every promoter in the genome. We assume
%that there is an average dG of binding for the non-target (other)
%promoters in the genome.
Sigma70_PromoterCompetition = Promoters * exp(-dG_Sigma70/RT) / (Promoters * exp(-dG_Sigma70/RT) + Other_Sigma70_Promoters*exp(-dG_Sigma_OtherPromoters/RT));
SigmaS_PromoterCompetition = Promoters * exp(-dG_SigmaS/RT) / (Promoters * exp(-dG_SigmaS/RT) + Other_SigmaS_Promoters*exp(-dG_Sigma_OtherPromoters/RT));

RNAP_Sigma70 = Total_RNAP_Sigma70 * Sigma70_PromoterCompetition;
RNAP_SigmaS = Total_RNAP_SigmaS * SigmaS_PromoterCompetition;

%3) Equilibrium binding of each repressor to its inducer (if
%present/existing) and to random DNA

Repressor1 = Total_Repressor1 / (1 +  Inducer1 * Ka_Inducer1 * molec2M) / Random_DNA;
Repressor2 = Total_Repressor2 / (1 +  Inducer2 * Ka_Inducer2 * molec2M) / Random_DNA;

%Energies & Microstates of Each State

dG(1) = 0;                  %Reference State
h(1) = 1;

dG(2) = dG_Sigma70;         %RNAP/Sigma70 bound to Promoter
h(2) = RNAP_Sigma70/Random_DNA;

dG(3) = dG_SigmaS;          %RNAP/SigmaS bound to Promoter
h(3) = RNAP_SigmaS/Random_DNA;

dG(4) = dG_Repressor1;
h(4) = Repressor1;

dG(5) = dG_Repressor2;
h(5) = Repressor2;

dG(6) = dG_Repressor1 + dG_Repressor2 + dG_R12;
h(6) = Repressor1 * Repressor2; 

P = exp(-dG./RT).*h;      
Q = sum(P);
Pinit = sum(P([2:3])) / Q;

r_init = Promoters * Pinit * k_conf; %mRNAs transcribed / min