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)')) %