function [output] = sigb_model_liebal(varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sigb_full-reference_0902
% Generated: 24-Aug-2009 16:55:24
% 
% [output] = sigb_model_liebal() => output = initial conditions in column vector
% [output] = sigb_model_liebal('states') => output = state names in cell-array
% [output] = sigb_model_liebal('algebraic') => output = algebraic variable names in cell-array
% [output] = sigb_model_liebal('parameters') => output = parameter names in cell-array
% [output] = sigb_model_liebal('parametervalues') => output = parameter values in column vector
% [output] = sigb_model_liebal(time,statevector) => output = time derivatives in column vector
% 
% State names and ordering:
% 
% statevector(1): X    (Biomass)
% statevector(2): Glc  (Glucose)
% statevector(3): B    (sigB,BSU04730)
% statevector(4): W2   (RsbW,BSU04720)
% statevector(5): W2B
% statevector(6): V    (RsbV,BSU04710)
% statevector(7): W2V
% statevector(8): W2V2
% statevector(9): VP
% statevector(10): Lacz(beta-galactosidase, recombinant reporter gene)
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global time
parameterValuesNew = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HANDLE VARIABLE INPUT ARGUMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin == 0,
	% Return initial conditions of the state variables (and possibly algebraic variables)
	output = [0.0225593, 4.56095, 0.0219636, 0.00029401, 0.825264, 5.5762e-06, 1.64632e-05, 0.285385, 1.56677, 6.86523];
	output = output(:);
	return
elseif nargin == 1,
	if strcmp(varargin{1},'states'),
		% Return state names in cell-array
		output = {'X', 'Glc', 'B', 'W2', 'W2B', 'V', 'W2V', 'W2V2', 'VP', 'Lacz'};
	elseif strcmp(varargin{1},'algebraic'),
		% Return algebraic variable names in cell-array
		output = {};
	elseif strcmp(varargin{1},'parameters'),
		% Return parameter names in cell-array
		output = {'k_turn', 'Kmdth', 'Kidth', 'Vmdth', 'Vmglc', 'Kmglc', 'nhglc', 'k_awbx', 'k_dwbx', 'k_awv1', ...
			'k_dwv1', 'n_awv2', 'n_dwv2', 'k_pho1', 'n_pho2', 'Vmdep', 'Kmdep', 'Kidep', 'nhdep', 'Kddpg', ...
			'k_bsyn', 'k_back', 'n_degl', 'n_degb', 'k_degs', 'k_xprZ', 'k_dgrZ', 'compart'};
	elseif strcmp(varargin{1},'parametervalues'),
		% Return parameter values in column vector
		output = [0.00396846, 31.7746, 0.00779687, 0.0626591, 0.0309228, 2.81366e-06, 1.74713, 0.999019, 0.00272426, 1.47964, ...
			0.00529827, 32.9994, 0.00911102, 0.00643309, 2.56261, 0.0212604, 0.306731, 1.77878e-12, 2.92192, 9.85965e-10, ...
			0.00108488, 3.32338e-05, 0.400092, 9.99997, 0.000762165, 356.023, 0.314823, 1];
	else
		error('Wrong input arguments! Please read the help text to the ODE file.');
	end
	output = output(:);
	return
elseif nargin == 2,
	time = varargin{1};
	statevector = varargin{2};
elseif nargin == 3,
	time = varargin{1};
	statevector = varargin{2};
	parameterValuesNew = varargin{3};
	if length(parameterValuesNew) ~= 28,
		parameterValuesNew = [];
	end
elseif nargin == 4,
	time = varargin{1};
	statevector = varargin{2};
	parameterValuesNew = varargin{4};
else
	error('Wrong input arguments! Please read the help text to the ODE file.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STATES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = statevector(1);
Glc = statevector(2);
B = statevector(3);
W2 = statevector(4);
W2B = statevector(5);
V = statevector(6);
W2V = statevector(7);
W2V2 = statevector(8);
VP = statevector(9);
Lacz = statevector(10);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(parameterValuesNew),
	k_turn = 0.00396846;
	Kmdth = 31.7746;
	Kidth = 0.00779687;
	Vmdth = 0.0626591;
	Vmglc = 0.0309228;
	Kmglc = 2.81366e-06;
	nhglc = 1.74713;
	k_awbx = 0.999019;
	k_dwbx = 0.00272426;
	k_awv1 = 1.47964;
	k_dwv1 = 0.00529827;
	n_awv2 = 32.9994;
	n_dwv2 = 0.00911102;
	k_pho1 = 0.00643309;
	n_pho2 = 2.56261;
	Vmdep = 0.0212604;
	Kmdep = 0.306731;
	Kidep = 1.77878e-12;
	nhdep = 2.92192;
	Kddpg = 9.85965e-10;
	k_bsyn = 0.00108488;
	k_back = 3.32338e-05;
	n_degl = 0.400092;
	n_degb = 9.99997;
	k_degs = 0.000762165;
	k_xprZ = 356.023;
	k_dgrZ = 0.314823;
	compart = 1;
else
	k_turn = parameterValuesNew(1);
	Kmdth = parameterValuesNew(2);
	Kidth = parameterValuesNew(3);
	Vmdth = parameterValuesNew(4);
	Vmglc = parameterValuesNew(5);
	Kmglc = parameterValuesNew(6);
	nhglc = parameterValuesNew(7);
	k_awbx = parameterValuesNew(8);
	k_dwbx = parameterValuesNew(9);
	k_awv1 = parameterValuesNew(10);
	k_dwv1 = parameterValuesNew(11);
	n_awv2 = parameterValuesNew(12);
	n_dwv2 = parameterValuesNew(13);
	k_pho1 = parameterValuesNew(14);
	n_pho2 = parameterValuesNew(15);
	Vmdep = parameterValuesNew(16);
	Kmdep = parameterValuesNew(17);
	Kidep = parameterValuesNew(18);
	nhdep = parameterValuesNew(19);
	Kddpg = parameterValuesNew(20);
	k_bsyn = parameterValuesNew(21);
	k_back = parameterValuesNew(22);
	n_degl = parameterValuesNew(23);
	n_degb = parameterValuesNew(24);
	k_degs = parameterValuesNew(25);
	k_xprZ = parameterValuesNew(26);
	k_dgrZ = parameterValuesNew(27);
	compart = parameterValuesNew(28);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VARIABLES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vfree = VP+V;
Vtot = VP+V+W2V+2*W2V2;
Wtot = 2*(W2+W2B+W2V+W2V2);
Btot = B+W2B;
glcinh = nhdep*Glc^nhdep/(Kddpg+Glc^nhdep);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REACTION KINETICS 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R_bioequ = k_turn*X*Glc-Vmdth*X/(Kmdth*(1+Glc/Kidth+X/Kmdth*(1+Glc/Kidth)));
R_glccon = sign(Glc)*Vmglc*Glc^nhglc/(Kmglc+Glc^nhglc);
R_w2b1eq = 2*k_awbx*B*W2-k_dwbx*W2B;
R_w2v1eq = 2*k_awv1*W2*V-k_dwv1*W2V;
R_w2v2eq = n_awv2*k_awv1*W2V*V-2*n_dwv2*k_dwv1*W2V2;
R_dephos = Vmdep*VP/(Kmdep*(1+glcinh/Kidep+VP/Kmdep*(1+glcinh/Kidep)));
R_w2v1ph = k_pho1*W2V;
R_w2v2ph = 2*n_pho2*k_pho1*W2V2;
R_prosyn = k_back+k_bsyn*B;
R_lacZ = k_xprZ*B-k_dgrZ*Lacz;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DIFFERENTIAL EQUATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X_dot = R_bioequ;
Glc_dot = -X*R_glccon;
B_dot = R_prosyn-R_w2b1eq-n_degb*k_degs*B;
W2_dot = R_w2v1ph-R_w2b1eq-R_w2v1eq+1/2*R_prosyn-n_degl*k_degs*W2;
W2B_dot = R_w2b1eq-n_degl*k_degs*W2B;
V_dot = R_dephos-R_w2v1eq-R_w2v2eq+R_prosyn-n_degl*k_degs*V;
W2V_dot = R_w2v1eq-R_w2v2eq-R_w2v1ph+R_w2v2ph-n_degl*k_degs*W2V;
W2V2_dot = R_w2v2eq-R_w2v2ph-k_degs*W2V2;
VP_dot = R_w2v1ph+R_w2v2ph-R_dephos-n_degl*k_degs*VP;
Lacz_dot = R_lacZ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RETURN VALUES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STATE ODEs
output(1) = X_dot;
output(2) = Glc_dot;
output(3) = B_dot;
output(4) = W2_dot;
output(5) = W2B_dot;
output(6) = V_dot;
output(7) = W2V_dot;
output(8) = W2V2_dot;
output(9) = VP_dot;
output(10) = Lacz_dot;
% return a column vector 
output = output(:);
return