%Phase diagram for the mutual repression switch

%You can run this first and then run "SwitchIntegration" if you want to
%display the simulated trajectory on top of the phase diagram.

clear all

%Parameters:
r=30;      %Production at maximum mRNA number in 1/min
gamma=1/30; %Degradation rate in 1/min
Kd=200;     %Dissociation constant

%Create vectors with numbers equally spaced between 0 and 1000. These are
%the values we'll use for R1 and R2.
R1=linspace(0,1000);
R2=linspace(0,1000);

%Plot the phase diagram
figure(1)
plot(R1,r/gamma./(1+(R1/Kd).^2),'-r')
hold on
plot(r/gamma./(1+(R2/Kd).^2),R2,'-k')
hold off

%In order to show the vectors we need to create a grid.
[R1M,R2M]=meshgrid(R1(1:2:end),R2(1:2:end));

%Calculate the x and y components fo rh derivative at each R1 and R2 point
%in phase space.
dR1M=-gamma*R1M+r./(1+(R2M/Kd).^2);
dR2M=-gamma*R2M+r./(1+(R1M/Kd).^2);

hold on
quiver(R1M,R2M,dR1M,dR2M,1.5,'LineWidth',1)
hold off

