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

% initial conditions
E = [10 10 0];
% model parameters
VR = 0.1; VB = 0.1; Kb = 2; Kr = 2;
dt = 0.5;
c = 1;
residual = 1;
% loop until convergence
while abs(residual) > 0.000001
   NE = E;
   % update concentrations
   NE = NE + dt*[ -VR*E(1)/(Kr+E(1))+VB*E(3)/(Kb+E(3)) VR*E(1)/(Kr+E(1)) -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;
   c = c+1;
end
stableactivity = E(3);    

end

