function ss_5genes();
alpha=[5 10 10 8 10 0  0 0 ];
beta=[10 10 10  10 10 0 0 0];
g=zeros(8,8);
h=zeros(8,8);
g(1,3)=1;
g(1,5)=-1;
g(2,1)=2;
g(3,2)=-1;
g(4,3)=2;
g(4,5)=-1;
g(5,4)=2;
g(1,6)=1;
g(2,7)=1;
g(3,8)=1;
g(4,6)=1;
g(5,7)=1;
h(1,1)=2;
h(2,2)=2;
h(3,2)=-1;
h(3,3)=2;
h(4,4)=2;
h(5,5)=2;
p0=[0.10 0.70 0.10 0.70 0.10 1.25 1.00 0.75];
tspan = [0:(0.5/40):0.5];
[t,p] = ode15s(@ode_ssys, tspan, p0,[],alpha,beta,g,h);
figure(1);
plot(t,p(:,1),'rs-');; 
hold on;
plot(t,p(:,2),'gd-');
plot(t,p(:,3),'ko-');
plot(t,p(:,4),'mx-');
plot(t,p(:,5),'b*-');
legend('X1','X2','X3','X4','X5');
xlabel('Time')
xlabel('Concentration')
title('A genetic network');


function dp = ode_ssys(t,X,alpha,beta,g,h);
dp = zeros(8,1);   
for i=1:5   
  vplus=alpha(i); 
  for j=1:8
     vplus=vplus*X(j)^g(i,j);
  end
  vminus=beta(i);  
  for j=1:8
     vminus=vminus*X(j)^h(i,j);
  end
  dp(i)=vplus-vminus;
end
dp(6)=0;
dp(7)=0;
dp(8)=0;
