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=(3)*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
 
 R=logspace(-6, 5,100); %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)
 
 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) 
 
 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
 
 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(2,1,1),semilogx(R,Fold_Change1,'k.', R,fit1,'k-', R,Fold_Change3,'b.', R,fit3,'b-',     
 R,Fold_Change5,'r.',R,fit5,'r-')
 semilogx(R,Fold_Change1,'k.', R,fit1,'k-', R,Fold_Change3,'b.', R,fit3,'b-', R,Fold_Change5,
 'r.',R,fit5,'r-')
 axis([10^-6 10^5 0 1])
 title(char('Thermodynamic Model and Apparent Hill Coefficient','    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('1x Binding Site Model','1x Binding Site Hill Function Fit', '3x Binding Site Model', '3x Binding Site Hill   
 Function Fit',
 '5x Binding Site Model','5x Binding Site Hill Function Fit','Location','NorthEastOutside')
 
 
 
 % subplot(2,1,2),loglog(R,Fold_Change1,'k--', R,Fold_Change3,'b--', R,Fold_Change5,'r--')
 figure(2)
 loglog(R,Fold_Change1,'k--', R,Fold_Change3,'b--', R,Fold_Change5,'r--')
 axis([10^-6 10^5 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)'))