function [ stableactivity ] = Barkai( alpha )
% uses a very crude ODE solver to identify the steady state of a
% methylation-activation model with integral control
% takes alpha, the activation probability of methylated protein
% returns steady-state activity

% initial conditions
E = [100 10 0];
% model parameters
VR = 0.1; VB = 1; Kb = 2; 
dt = 0.5;
rec = zeros(length(0:dt:100), 4);
c = 1;
residual = 1;
% loop until convergence
while abs(residual) > 0.000001
   NE = E;
   % update concentrations
   NE = NE + dt*[ -VR+VB*E(3)/(Kb+E(3)) VR -VB*E(3)/(Kb+E(3)) ];
   temp = NE(2)+NE(3);
   NE(2) = (1-alpha)*temp;
   NE(3) = alpha*temp;
   residual = NE(3)-E(3);
   E = NE;
   rec(c,:) = [c*dt E(1) E(2) E(3)];
   c = c+1;
end
stableactivity = E(3);    

end

