Mathematical Modeling of Bistable Toggle Switch

Nns=5*10^6; %Number of non-specific sites
Kb=1.3806*10^-23; %(JK-1) %Boltzman's constant
T=298; %(K) %Temperature in Kelvin
P=3000; %Number of RNA polymerase (2000-4000)

Kspd=100*10^-9; %dissociation constant for specific binding of polymerase on DNA
Knspd=10000; %dissociation constant for non-specific binding of polymerase on DNA
delEpd=Kb*T*log(Kspd/Knspd); %Binding energy between polymerase and DNA

Ksrd=(10)*10^-9; %TALE %dissociation constant for specific binding of repressor on DNA (TALE: 1/150nM~, paper
example: 0.02)
Knsrd=10000; %dissociation constant for non-specific binding of repressor on DNA
delErd=Kb*T*log(Ksrd/Knsrd); %Binding energy between repressor and DNA

Ksrd_strong=(0.1)*10^-9; %TALE %dissociation constant for specific binding of repressor on DNA (TALE: 1/150nM~,
paper example: 0.02) --> weaker binding (x0.1)
delErd_strong=Kb*T*log(Ksrd_strong/Knsrd); %Binding energy between repressor and DNA

R=logspace(-6.5, -3,500); %Number of repressors
% R=linspace(0,100);
Freg1=1./(1+(R./Nns).*exp(-delErd/(Kb*T))); %Regulation factor 1x (<1 for repression)
Freg3=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^3; %Regulation factor 3x (<1 for repression)
Freg5=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^5; %Regulation factor 5x (<1 for repression)

Freg5_0=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^5.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^0; %Regulation
factor 5x (<1 for repression)
Freg4_1=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^4.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^1;  %Regulation
factor 5x (<1 for repression)
Freg3_2=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^3.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^2;  %Regulation
factor 5x (<1 for repression)
Freg2_3=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^2.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^3;  %Regulation
factor 5x (<1 for repression)
Freg1_4=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^1.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^4;  %Regulation
factor 5x (<1 for repression)
Freg0_5=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^0.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^5;  %Regulation
factor 5x (<1 for repression)

p_bound_base=1./(1+(Nns/P)*exp(delEpd/(Kb*T))); %basal level
p_bound_1=1./(1+(Nns./(P.*Freg1)).*exp(delEpd/(Kb*T))); %level with repressor bound (1x)
p_bound_3=(1./(1+(Nns./(P.*Freg3)).*exp(delEpd/(Kb*T)))); %level with repressor bound (3x)
p_bound_5=(1./(1+(Nns./(P.*Freg5)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x)

p_bound_5_0=(1./(1+(Nns./(P.*Freg5_0)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x)
p_bound_4_1=(1./(1+(Nns./(P.*Freg4_1)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x)
p_bound_3_2=(1./(1+(Nns./(P.*Freg3_2)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x)
p_bound_2_3=(1./(1+(Nns./(P.*Freg2_3)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x)
p_bound_1_4=(1./(1+(Nns./(P.*Freg1_4)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x)
p_bound_0_5=(1./(1+(Nns./(P.*Freg0_5)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x)



Fold_Change1=p_bound_1./p_bound_base; %Fold Change 1x
Fold_Change3=p_bound_3./p_bound_base; %Fold Change 3x
Fold_Change5=p_bound_5./p_bound_base; %Fold Change 5x

Fold_Change5_0=p_bound_5_0./p_bound_base; %Fold Change 5x
Fold_Change4_1=p_bound_4_1./p_bound_base; %Fold Change 5x
Fold_Change3_2=p_bound_3_2./p_bound_base; %Fold Change 5x
Fold_Change2_3=p_bound_2_3./p_bound_base; %Fold Change 5x
Fold_Change1_4=p_bound_1_4./p_bound_base; %Fold Change 5x
Fold_Change0_5=p_bound_0_5./p_bound_base; %Fold Change 5x

fit1=1./(0.0111.*(R).^1+1);
fit3=1./(4.7215e9.*(R).^2.9929+1);
fit5=1./(6.6001e20.*(R).^4.8635+1);


figure1=figure(1);
subplot(1,1,1),semilogx(R,Fold_Change0_5,'k.',
R,Fold_Change1_4,'b.',R,Fold_Change2_3,'g.',R,Fold_Change3_2,'y.',R,Fold_Change4_1,'r.',R,Fold_Change5_0,'m.')
axis([10^-6.5 10^-3 0 1])
title(char(' Thermodynamic Model and Apparent Hill Coefficient','Unbalanced Binding Strengths : Fold Change vs Repressor Level (Lin - Log)'))
xlabel('R (number of repressor molecules)')
ylabel('Fold-Change of P_b_o_u_n_d')
%ylabel(char(' Fold-Change of P_b_o_u_n_d','(probability of Polymerase binding to promoter)'))
legend('5xStrong, 0xWeak','4xStrong, 1xWeak', '3xStrong, 2xWeak', '2xStrong, 3xWeak','1xStrong, 4xWeak','0xStrong, 5xWeak','Location','NorthEastOutside')


% subplot(2,1,2),loglog(R,Fold_Change5_0,'k.',
R,Fold_Change4_1,'b.',R,Fold_Change3_2,'g.',R,Fold_Change2_3,'y.',R,Fold_Change1_4,'r.',R,Fold_Change0_5,'m.')
% axis([10^-6 10^-3 10^-3 10^0])
% title(char(' Thermodynamic Model and Sensitivity','Fold Change vs Repressor Level (Log - Log)'))
% legend('1x Binding Site', '3x Binding Site', '5x Binding Site','Location','NorthEastOutside')
% xlabel('R (number of repressor molecules)')
% ylabel('Fold-Change of P_b_o_u_n_d')
% % ylabel(char(' Fold-Change of P_b_o_u_n_d','(probability of Polymerase binding to promoter)'))
%