function [t,x,x0,v,Y0,Finv] = daes_linearchain_dae(loadfunctions);

% linear chain, to study the enzyme synthesis rates

if nargin==0
  loadfunctions=0;
end

format compact
n = 7; % number of metabolites
nenz = 6; % number of reactions/enzymes
% parameters
A = [
3 1 3 3 1 2 3 1 3 2; 
5 1 3 3 2 3 2 1 3 4 ; 
4 6 5 1 3 4 1 3 3 1;
3 2 1 4 1 2 3 2 1 2];
BCs = [10;1]; %externally fixed concentrations
sensors = [2,4]; % note that sensors(1) == 1 <=> sensor = x0, niet x1!

% the stoichiometry matrix
N = [0  0  0  0  0  0;
	 1 -1  0  0  0  0;
	 0  1 -1  0  0  0;
	 0  0  1 -1  0  0;
	 0  0  0  1 -1  0;
	 0  0  0  0  1 -1;
	 0  0  0  0  0  0];

e = [0.23 0.2 0.1 0.15 0.05 0.17];
es = sum(e); e=e/es; % normalise to sum=1
params{1} = A;
params{2} = e;
params{3} = n;
params{4} = BCs;
params{5} = sensors;
params{6} = N;
params{7} = length(e);

if loadfunctions
  disp('... loading Symbolic functions ...')
  load('Symbolic_functions/daes_linearchain_dae_functions.mat');
else
  %define the object function symbolically, and compute its symbolic Jacobian, 
% which yields the optimum equations.
  syms Objsym Obj_parts dObjdx dObjdxx f1 f2 f3 f4 f5 f6 x1 x2 x3 x4 x5 x6 x7 
  vars=[x1 x2 x3 x4 x5 x6 x7];
  f1=symfun((x1-x2)/(A(1,1) + A(2,1)*x1 + A(3,1)*x2), vars);
  f2=symfun((x2-x3)/(A(1,2) + A(2,2)*x2 + A(3,2)*x3), vars);
  f3=symfun((x3-x4)/(A(1,3) + A(2,3)*x3 + A(3,3)*x4), vars);
  f4=symfun((x4-x5)/(A(1,4) + A(2,4)*x4 + A(3,4)*x5), vars);
  f5=symfun((x5-x6)/(A(1,5) + A(2,5)*x5 + A(3,5)*x6), vars);
  f6=symfun((x6-x7)/(A(1,6) + A(2,6)*x6 + A(3,6)*x7), vars);
  prod = symfun(((x1-x2)*(x2-x3)*(x3-x4)*(x4-x5)*(x5-x6)*(x6-x7)),vars);

  disp('... computing the Jacobians....')
  FTNS      = symfun([f1;f2;f3;f4;f5;f6],vars);
  OBJ_PARTS = symfun([1/f1; 1/f2; 1/f3; 1/f4; 1/f5; 1/f6],vars);
  OBJSYM    = symfun(1/f1 + 1/f2 + 1/f3 + 1/f4 + 1/f5 + 1/f6, vars);
  DOBJDX    = symfun(jacobian(OBJSYM,[x2 x3 x4 x5 x6]),vars);
  DOBJDXPROD= symfun(prod*DOBJDX,vars);
%  DOBJDXPROD= symfun(DOBJDX,vars);
  DOBJDXX   = symfun(jacobian(prod*DOBJDX,vars),vars);
%  ENZ       = symfun((x1-x2)*[1/f1; 1/f2; 1/f3; 1/f4; 1/f5; 1/f6]/((x1-x2)*OBJSYM),vars);
  ENZ       = symfun([1/f1; 1/f2; 1/f3; 1/f4; 1/f5; 1/f6]/(OBJSYM),vars);

  % now make these normal matlab functions: much faster evaluation
  ftns      = matlabFunction(FTNS);
  Obj_parts = matlabFunction(OBJ_PARTS);
  Objsym    = matlabFunction(OBJSYM);
  dObjdx    = matlabFunction(DOBJDX);
  dObjdxprod= matlabFunction(DOBJDXPROD);
  dObjdxx   = matlabFunction(DOBJDXX);
  Enz	  	= matlabFunction(ENZ);
  save('Symbolic_functions/daes_linearchain_dae_functions.mat', 'vars','ftns', 'Obj_parts','Objsym','dObjdx','dObjdxprod','dObjdxx','Enz');
end

% stitching parts together
fullxvector=@(x,BCs) [BCs(1) x(1:end)',BCs(2)];


disp('... computing the QSSA....')
% first we do a QSSA on x. No enzyme dynamics. The result is a well-ordered 
% concentration vector
x0 = 10*rand(n-length(BCs),1);
%X0 = fullxvector(x0,BCs);
%Metsys = @(t,x) metsys(t,x,ftns,params);
%options = odeset('RelTol',1e-6,'AbsTol',1e-8);
%[t,x0] = ode45(Metsys,[0,1000],X0,options);
%x0 = x0(end,[2:end-1])';

disp('... computing the initial optimum....')
%compute the initial optimum
%zopt = compute_initial_optimum([BCs(1) x0(:)' BCs(2)],A,n,sensors);
%zopt(zopt<0)=0
% in this case, we just take some arbitrary optimum, using fake BCs, which 
% are twice the original ones. they could be anything (realistic)
BCs0 = [3 9];
Obj = @(x0) -sum(Obj_parts(BCs0(1),x0(1),x0(2),x0(3),x0(4),x0(5),BCs0(2)));
z0 = fliplr([8 7 6 5 4]); % just some random IC that is in the admissable set.
opts = optimset('fminunc');
opts = optimset(opts,'DiffMaxChange',1);
opts = optimset(opts,'TolFun',eps);
opts = optimset(opts,'TolX',eps,'Display','off');
zopt = fminunc(Obj,z0,opts);
zopt = zopt(:);
zopt = [fullxvector(zopt,BCs0)];

% now take the sensors out of this, and replace the values in x0 with these values
x0(sensors-1) = zopt(sensors);
dObjdx(zopt(1),zopt(2),zopt(3),zopt(4),zopt(5),zopt(6),zopt(7))
%dz = dobjdx(zopt([2:end-1]),A,n,[zopt(1) zopt(3) zopt(end)]);

X0 = fullxvector(x0,BCs);
x0 = [X0(:); e(:); zopt(:)]; % length 3n-1;
x0 = x0(:)'

disp('... integrating the ODE ....')
Mass = zeros(2*n+nenz); % n x's, n-2 e's, n z's
Mass(1:n+nenz,1:n+nenz) = eye(n+nenz);
ODE = @(t,x) ode(t,x,ftns,Obj_parts,dObjdxprod,dObjdxx,Enz,params);
options = odeset('Mass',Mass,'MStateDependence','none','RelTol',1e-8,'AbsTol',1e-8);
[t,x] = ode15s(ODE,[0,300],x0,options);
ODE(t(end),x(end,:)')'

for i=1:length(t)
  XX = x(i,1:n);
  f = ftns(XX(1),XX(2),XX(3),XX(4),XX(5),XX(6),XX(7));
  v(i,1:nenz) =x(i,n+1:n+nenz).*f';
end

plot_optima_branched(t,x,v,ftns,params,1);
data = [t,x,v];
save('daes_linearchain_reversal_data.txt','-ASCII','data');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dx = metsys(t,x,ftns,params);

A = params{1};
e = params{2};
n = params{3};
BCs = params{4};
sensors = params{5};
N = params{6};

F = ftns(x(1),x(2),x(3),x(4),x(5),x(6),x(7));
V = e(:).*F(:);

dx = N*V(:); 
dx = dx(:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dxe = ode(t,X,ftns,Obj_parts,dObjdxprod,dObjdxx,Enz,params);


A = params{1};
n = params{3};
BCs = params{4};
sensors = params{5};
N = params{6};
nenz = params{7};
nsens = length(sensors);

x = X(1:n); % metabolites
e = X(n+1:n+nenz); % enzymzes
z = X(n+nenz+1:end); % raw sensor optimum

F = ftns(x(1),x(2),x(3),x(4),x(5),x(6),x(7)); 
V = e(:).*F(:);

dx = N*V(:);


% we make the DAE explicit, by converting it to an ODE
%dxsensors = [zeros(n-nsens,1); dx(sensors)];
%DZ = dObjdxx(z(1),z(2),z(3),z(4),z(5),z(6),z(7));
%[DZr,DZc] = size(DZ);
%F = [DZ;zeros(nsens,n)];
%for i=1:nsens,
%  F(DZr+i,sensors(i))=-1; 
%end
%condest(F);
%dz =  -F \ dxsensors;
% the optimum is tracked using two sensors

dxy(1) = x(sensors(1)) - z(sensors(1));
dxy(2) = x(sensors(2)) - z(sensors(2));
dz = dObjdxprod(z(1),z(2),z(3),z(4),z(5),z(6),z(7));

Eopt = Enz(z(1),z(2),z(3),z(4),z(5),z(6),z(7));
%som=sum(Eopt);
for i=1:n-1
  de(i) = 1/1000*(Eopt(i) - e(i));
end

[x(:)' e(:)' z(:)'];

dxe = [dx(:);de(:);dxy(:);dz(:)];
dxe(:)';
dz;

%%%%%%%%%%%%

function plot_optima_branched(t,y,v,ftns,params,doprint)

if nargin == 3
  doprint = 0;
end

A = params{1};
n = params{3};
BCs = [1,7];
sensors = params{5};
N = params{6};
nenz = params{7};

x = y(:,1:n);
e = y(:,n+1:n+nenz);
opt = y(:,n+nenz+1:end);
rest = setdiff([1:n],BCs);
sensor = [2,4];
rest = setdiff(rest,sensors);

figure(1)
  clf
  hold on
  plot([0,t(end)],[0,0],'k-')
  for i=1:size(v,2)
    q(i) = plot(t,v(:,i),'LineWidth',2);
  end
  set(gca,'Box','on')
  hold off
  legend([q],'v_1','v_2','v_3','v_4','v_5','v_6','v_7')
  %xlabel('time')
  %ylabel('reaction flux')

%  title ('enzymes')
figure(2)
  clf
  hold on
  p1=plot(t,x(:,rest),'r-','LineWidth',2);
  p2=plot(t,opt(:,1:end),'b-','LineWidth',2);
  plot(t,x(:,sensors),'r--','LineWidth',2);
  p4=plot(t,[x(:,BCs) ],'Color',[0 0.6 0],'LineWidth',2);
  set(gca,'Box','on')
  hold off
  grid off
  set(gca,'YLim',[0,max(x(:))*1.1]);
  %xlabel('time');
  %ylabel('metabolite/sensor concentration');
  legend([p1(1) p4(1) p2(1)],'Internal metabolites', 'External Metabolites', 'Sensor optimum')

if doprint
  print_fig(1, 'pics/through_therm_equil_fluxes.jpg',[0 0 5 3]); 
  print_fig(2, 'pics/through_therm_equil.jpg',[0 0 5 3]); 
end
