%% This code produces one-parameter bifurcation diagram for the 5D Lactose Operon Model
function Code_for_Ex_2_5
clc
%Define a range of allolactose concentration as,
ssA            = [10^(-8):10^(-3):1e0];
% Choose an initial steady state concentration for mRNA , beta-gal,
% internal lactose, external lactose and permease,
init           = [1 1 1 1 1]*0.006;
%Create a foor loop and compute steady state for each allolactose concentration
% above with Matlab's fsolve command
for i=1:length(ssA)
    As              = ssA(i);
    SSValues(i,:)   = fsolve(@bifur,init,[],As);    
    ssM             = SSValues(:,1);
    ssB             = SSValues(:,2);
    ssL             = SSValues(:,3);
    ssLe            = SSValues(:,4);
    ssP             = SSValues(:,5);
    init            = SSValues(i,:) ;
end
 figure(1),semilogy(ssLe,ssA,'linewidth',2), hold on
% Now plot steady state curve for external lactose on x-axis and
% allolactose on y-axis. Choosey-axis in log scale.
............
axis([0.001 0.07 0.001 1])




%% Function for the Steady State Equations
function dy = bifur(y,As)
% Input the parameter values below from the paper published in Biophysical Journal by Yildirim et al in 2003 
..........................


% Input the steady state equations below from the same paper, dA is the equation for allolactose, dM is for mRNA etc 
dA = ............;
dM = ............;
dB = ............;
dL = ............;
dP = ............;
dy = [dA dM dB dL dP]';

