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')
   legend('5xStrong, 0xWeak','4xStrong, 1xWeak', '3xStrong, 2xWeak', '2xStrong, 3xWeak','1xStrong,
  4xWeak','0xStrong, 5xWeak','Location','NorthEastOutside')