% Systems Biology		7.81/8.591/9.531		Fall 2004
% Problem set 1
%
% Problem 3: Positive feedback and multistability 
% file: ps1.m

clear; close;

% keep the following variables at their default values
Dtot=1;	% concentration of DNA
A=1;		% association constant
gamma=1;	% degradation rate

% change the values of the following variables
v= ;		% basal transcription rate
k= ;		% maximal transcription rate
n= ;		% Hill coefficient
			% n=1 for Scheme A, n=2 for Scheme B

% we later solve for roots, a procedure that requires 
% that n be an integer
if( rem(n,1)~=0 ) 
   error('n should be an integer');
end
            
% these variables define the range of x values;
% change them from the default to 'zoom in' on interesting values
xmin=0;	% lower limit
xmax=2;	% upper limit
nx=1000;	% number of points in mesh
x=linspace(xmin,xmax,nx);

% generation rate
gen=v+k*Dtot*(x.^n)./(A^n+x.^n);
% degradation rate
deg=gamma*x;

% plot curves
plot(x,deg,x,gen);
xlabel('concentration (x)');
ylabel('rate of change (dx/dt)');
title('Generation and Degradation of X');
legend('deg(x)','gen(x)',4);

% now solve for the roots of the system
if(n==1)
   c(1)=gamma;
   c(2)=-(v+k*Dtot)+gamma*A^n;
   c(3)=-v*A^n;
else
	c(1)=gamma;
	c(2)=-(v+k*Dtot);
	c(n+1)=gamma*A^n;
   c(n+2)=-v*A^n;
end
xsol=roots(c);

% count and plot roots
nsol=length(xsol);
nroots=0;
for i=1:nsol
   if( (imag(xsol(i))==0)&(real(xsol(i))>=0) )
      nroots=nroots+1;
      xbar=[xsol(i);xsol(i)];
      ybar=[0;gamma*xsol(i)];
      hl=line(xbar,ybar);
      set(hl,'Color',[1 0 0]);
   end
end
if( nroots==1 )
   etext=[''];
else
   etext=['s'];
end
nrtext=['The system has ' int2str(nroots) ' root' etext '.'];
text(xmin+0.1*(xmax-xmin),0.9*xmax,nrtext)