function r_init = Gene_Basal()
%Chemical partition function for basal gene expression

%Conversion factors & constants
molec2uM = 0.00166057788110;
RT = 1.987e-3*310;

%Plasmid / chromosome copy number
copynumber = 100;

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

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

%Constant Concentrations -- do not change
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

%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

%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;

%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;

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

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