Team:Duke/Modeling/Codes/Thermo2

{{Template:Team:Duke}}

<div class="right_box">

= '''Mathematical Modeling of Bistable Toggle Switch''' =

<br>Nns=5*10^6; %Number of non-specific sites
<br>Kb=1.3806*10^-23; %(JK-1) %Boltzman's constant
<br>T=298; %(K) %Temperature in Kelvin
<br>P=3000; %Number of RNA polymerase (2000-4000)
<br>
<br>Kspd=100*10^-9; %dissociation constant for specific binding of polymerase on DNA
<br>Knspd=10000; %dissociation constant for non-specific binding of polymerase on DNA
<br>delEpd=Kb*T*log(Kspd/Knspd); %Binding energy between polymerase and DNA
<br>
<br>Ksrd=(10)*10^-9; %TALE %dissociation constant for specific binding of repressor on DNA (TALE: 1/150nM~, paper <br>example: 0.02)
<br>Knsrd=10000; %dissociation constant for non-specific binding of repressor on DNA
<br>delErd=Kb*T*log(Ksrd/Knsrd); %Binding energy between repressor and DNA
<br>
<br>Ksrd_strong=(0.1)*10^-9; %TALE %dissociation constant for specific binding of repressor on DNA (TALE: 1/150nM~, <br>paper example: 0.02) --> weaker binding (x0.1)
<br>delErd_strong=Kb*T*log(Ksrd_strong/Knsrd); %Binding energy between repressor and DNA
<br>
<br>R=logspace(-6.5, -3,500); %Number of repressors
<br>% R=linspace(0,100);
<br>Freg1=1./(1+(R./Nns).*exp(-delErd/(Kb*T))); %Regulation factor 1x (<1 for repression)
<br>Freg3=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^3; %Regulation factor 3x (<1 for repression)
<br>Freg5=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^5; %Regulation factor 5x (<1 for repression)
<br>
<br>Freg5_0=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^5.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^0; %Regulation <br>factor 5x (<1 for repression)
<br>Freg4_1=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^4.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^1;  %Regulation <br>factor 5x (<1 for repression)
<br>Freg3_2=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^3.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^2;  %Regulation <br>factor 5x (<1 for repression)
<br>Freg2_3=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^2.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^3;  %Regulation <br>factor 5x (<1 for repression)
<br>Freg1_4=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^1.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^4;  %Regulation <br>factor 5x (<1 for repression)
<br>Freg0_5=(1./(1+(R./Nns).*exp(-delErd/(Kb*T)))).^0.*(1./(1+(R./Nns).*exp(-delErd_strong/(Kb*T)))).^5;  %Regulation <br>factor 5x (<1 for repression)
<br>
<br>p_bound_base=1./(1+(Nns/P)*exp(delEpd/(Kb*T))); %basal level
<br>p_bound_1=1./(1+(Nns./(P.*Freg1)).*exp(delEpd/(Kb*T))); %level with repressor bound (1x)
<br>p_bound_3=(1./(1+(Nns./(P.*Freg3)).*exp(delEpd/(Kb*T)))); %level with repressor bound (3x) 
<br>p_bound_5=(1./(1+(Nns./(P.*Freg5)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x) 
<br>
<br>p_bound_5_0=(1./(1+(Nns./(P.*Freg5_0)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x) 
<br>p_bound_4_1=(1./(1+(Nns./(P.*Freg4_1)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x) 
<br>p_bound_3_2=(1./(1+(Nns./(P.*Freg3_2)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x) 
<br>p_bound_2_3=(1./(1+(Nns./(P.*Freg2_3)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x) 
<br>p_bound_1_4=(1./(1+(Nns./(P.*Freg1_4)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x) 
<br>p_bound_0_5=(1./(1+(Nns./(P.*Freg0_5)).*exp(delEpd/(Kb*T)))); %level with repressor bound (5x) 
<br>
<br>
<br>
<br>Fold_Change1=p_bound_1./p_bound_base; %Fold Change 1x
<br>Fold_Change3=p_bound_3./p_bound_base; %Fold Change 3x
<br>Fold_Change5=p_bound_5./p_bound_base; %Fold Change 5x
<br>
<br>Fold_Change5_0=p_bound_5_0./p_bound_base; %Fold Change 5x
<br>Fold_Change4_1=p_bound_4_1./p_bound_base; %Fold Change 5x
<br>Fold_Change3_2=p_bound_3_2./p_bound_base; %Fold Change 5x
<br>Fold_Change2_3=p_bound_2_3./p_bound_base; %Fold Change 5x
<br>Fold_Change1_4=p_bound_1_4./p_bound_base; %Fold Change 5x
<br>Fold_Change0_5=p_bound_0_5./p_bound_base; %Fold Change 5x
<br>
<br>fit1=1./(0.0111.*(R).^1+1);
<br>fit3=1./(4.7215e9.*(R).^2.9929+1);
<br>fit5=1./(6.6001e20.*(R).^4.8635+1);
<br>
<br>
<br>figure1=figure(1);
<br>subplot(1,1,1),semilogx(R,Fold_Change0_5,'k.', <br>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.')
<br>axis([10^-6.5 10^-3 0 1])
<br>title(char('            Thermodynamic Model and Apparent Hill Coefficient','Unbalanced Binding Strengths : Fold Change vs Repressor Level (Lin - Log)'))
<br>xlabel('R (number of repressor molecules)')
<br>ylabel('Fold-Change of P_b_o_u_n_d')
<br>%ylabel(char('             Fold-Change of P_b_o_u_n_d','(probability of Polymerase binding to promoter)'))
<br>legend('5xStrong, 0xWeak','4xStrong, 1xWeak', '3xStrong, 2xWeak', '2xStrong, 3xWeak','1xStrong, 4xWeak','0xStrong, 5xWeak','Location','NorthEastOutside')
<br>
<br>
<br>% 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.')
<br>% axis([10^-6 10^-3 10^-3 10^0])
<br>% title(char('    Thermodynamic Model and Sensitivity','Fold Change vs Repressor Level (Log - Log)'))
<br>% legend('1x Binding Site', '3x Binding Site', '5x Binding Site','Location','NorthEastOutside')
<br>% xlabel('R (number of repressor molecules)')
<br>% ylabel('Fold-Change of P_b_o_u_n_d')
<br>% % ylabel(char('             Fold-Change of P_b_o_u_n_d','(probability of Polymerase binding to promoter)'))
<br>% 




<br><br><br>
</div>