function [t,x,x0,diffxopt] = daes_branched;

% as daes.m but for a branched EFM with two inputs

format compact
n = 12; % full vector, so n-4 internal metabolites. Roman numerals are the flux numbers
% x1->x2 \	/ x5 -> x6 -> x7 -> x8
%	I	  ==	 IV		V	 VI
% x3->x4 /III\ x9 -> x10 -> x11 -> x12
%   II 			  VII	 VIII   IX
% with x5 inhibiting x10->x11 and x10 inhibiting x6->x7

A = round(10*(rand(4,n)))+1; % just some 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;
0 0 1 4 5 2 3 1 2 1 
0 0 2 1 2 1 3 1 2 3];

loadfunctions = 0;
if loadfunctions == 0

  %define the object function symbolically, and compute its symbolic Jacobian, 
  % which yields the optimum equations.
  syms x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12;
  v=[x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12];
  f1=symfun((x1-x2)/(A(1,1) + A(2,1)*x1 + A(3,1)*x2), v);
  f2=symfun((x3-x4)/(A(1,2) + A(2,2)*x3 + A(3,2)*x4), v);
  f3=symfun((x2*x4-x5*x9)/(A(1,3) + A(2,3)*x2 + A(3,3)*x4 + A(4,3)*x5 + A(5,3)*x9), v);
  f4=symfun((x5-x6)/(A(1,4) + A(2,4)*x5 + A(3,4)*x6), v);
  f5=symfun((x6-x7)/(A(1,5) + A(2,5)*x6 + A(3,5)*x7 + A(4,5)*x10), v);
  f6=symfun((x7-x8)/(A(1,6) + A(2,6)*x7 + A(3,6)*x8), v);
  f7=symfun((x9-x10)/(A(1,7) + A(2,7)*x9 + A(3,7)*x10), v);
  f8=symfun((x10-x11)/(A(1,8) + A(2,8)*x10 + A(3,8)*x11 + A(4,8)*x5), v);
  f9=symfun((x11-x12)/(A(1,9) + A(2,9)*x11 + A(3,9)*x12), v);

  disp('... computing the Jacobians....')
  FTNS = symfun([f1;f2;f3;f4;f5;f6;f7;f8;f9],v);
  OBJ_PARTS = symfun([1/f1; 1/f2; 1/f3; 1/f4; 1/f5; 1/f6; 1/f7; 1/f8; 1/f9],v);
  OBJSYM = symfun(1/f1 + 1/f2 + 1/f3 + 1/f4 + 1/f5 + 1/f6 + 1/f7 + 1/f8 + 1/f9, v);
  DOBJDX = symfun(jacobian(OBJSYM,[x2 x4 x5 x6 x7 x9 x10 x11]),v); % all vars except BCs
  DOBJDXX = symfun(jacobian(DOBJDX,v),v);

  ftns      = matlabFunction(FTNS);
  Obj_parts = matlabFunction(OBJ_PARTS);
  Objsym    = matlabFunction(OBJSYM);
  dObjdx    = matlabFunction(DOBJDX);
  dObjdxx   = matlabFunction(DOBJDXX);
  
  save('double_branched.mat','v','ftns','Obj_parts','Objsym','dObjdx',...
	'dObjdxx','A');
else
  load('double_branched.mat');
end
%test the initial optimum guess first
BCs0 = [9 10 2 3]; % x1 x3 x8 x12
Obj = @(x) obj(x,Obj_parts,BCs0);
z0 = [8 8 5 4 3.5 6 5 4]; % [x2 x4 x5 x6 x7 x9 x10 x11]
disp('... making sure that the initial guess for the start optimum is admissable...')
disp('... i.e. the vector should be positive ...')
double(Obj_parts(BCs0(1),z0(1),BCs0(2),z0(2),z0(3),z0(4),z0(5),BCs0(3),...
	z0(6),z0(7),z0(8),BCs0(4)))								

e = [0.23 0.2 0.1 0.15 0.05 0.17 0.1 0.2 0.3];
es = sum(e); e=e/es;
BCs = [10;7;3;1]; %x(1), x(3), x(9)
sensors = [2,4,7,11]; 

x0 = 10*rand(n-length(BCs),1); % [x2 x4 x5 x6 x7 x9 x10 x11]
params{1} = A;
params{2} = e;
params{3} = n;
params{4} = BCs;
params{5} = sensors;

disp('... computing the QSSA....')
% first we do a QSSA on x. No enzyme dynamics. The result is a well-ordered 
% concentration vector
Metsys = @(t,x) metsys(t,x,ftns,params);
options = odeset('RelTol',1e-6,'AbsTol',1e-8);
[t,x0] = ode45(Metsys,[0,100],[BCs(1) x0(1) BCs(2) x0(2:5)' BCs(3) x0(6:8)' BCs(4)],options);
x=x0;
x0 = x0(end,[2 4: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)
% just some random IC that is in the admissable set.

opts = optimset('fminunc');
opts = optimset(opts,'DiffMaxChange',1);
zopt = fminunc(Obj,z0,opts);
dobjdx(zopt,Obj_parts,BCs0);
zopt = zopt(:)';
zopt = [BCs0(1) zopt(1) BCs0(2) zopt(2:5) BCs0(3) zopt(6:8) BCs0(4)];
% now take the sensors out of this, and replace the values in x0 with these values
x0(1) = zopt(2); % note that x(3) is a BC, so be careful with the numbering!
x0(2) = zopt(4);
x0(5) = zopt(7);
x0(8) = zopt(11);
double(dObjdx(zopt(1),zopt(2),zopt(3),zopt(4),zopt(5),zopt(6),zopt(7),zopt(8),zopt(9),zopt(10),zopt(11),zopt(12)))

%dz = dobjdx(zopt([2 4:end-1]),A,n,[zopt(1) zopt(3) zopt(end)]);

x0 = [BCs(1) ;x0(1) ;BCs(2); x0(2:5); BCs(3); x0(6:8); BCs(4); e(:); zopt(:)]; 
x0 = x0(:)'

disp('... integrating the ODE ....')
ODE = @(t,x) ode(t,x,ftns,Obj_parts,dObjdxx,params);
options = odeset('RelTol',1e-5,'AbsTol',1e-5);
[t,x] = ode15s(ODE,[0,1000],x0,options);
ODE(t(end),x(end,:)')';

BCs=[6,7,3,1];
params{4}=BCs;
x0=x(end,:);
x0(1)=BCs(1); 
%x0(n-1:n+nenz) = [0 0 0 0 0 0]; %knocking out all enzymes but one
[t1,x1] = ode15s(ODE,[t(end),t(end)+1000],x0,options);
t=[t;t1]; x=[x;x1];

BCs=[6,5,3,1];
params{4}=BCs;
x0=x(end,:);
x0(3)=BCs(2); 
%x0(n-1:n+nenz) = [0 0 0 0 0 0]; %knocking out all enzymes but one
[t1,x1] = ode15s(ODE,[t(end),t(end)+1000],x0,options);
t=[t;t1]; x=[x;x1];

BCs=[6,5,2,1];
params{4}=BCs;
x0=x(end,:);
x0(8)=BCs(3); 
%x0(n-1:n+nenz) = [0 0 0 0 0 0]; %knocking out all enzymes but one
[t1,x1] = ode15s(ODE,[t(end),t(end)+1000],x0,options);
t=[t;t1]; x=[x;x1];

BCs=[6,5,2,0];
params{4}=BCs;
x0=x(end,:);
x0(12)=BCs(4); 
%x0(n-1:n+nenz) = [0 0 0 0 0 0]; %knocking out all enzymes but one
[t1,x1] = ode15s(ODE,[t(end),t(end)+1000],x0,options);
t=[t;t1]; x=[x;x1];


diffxopt = plot_optima_branched(t,x,1);

data = [t x diffxopt];
save ('daes_double_branched_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};

F = double(ftns(x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10),x(11),x(12)));
for i=1:length(F)
  V(i) = e(i)*F(i);
end

dx(1) = 0; % The BCs are incorporated as dx1/dt = 0 and dxn/dt=0
dx(2) = V(1) - V(3);
dx(3) = 0;
dx(4) = V(2) - V(3);
dx(5) = V(3) - V(4);
dx(6) = V(4) - V(5);
dx(7) = V(5) - V(6);
dx(8) = 0;
dx(9) = V(3) - V(7);
dx(10) = V(7) - V(8);
dx(11) = V(8) - V(9);
dx(12) = 0;

dx = dx(:);

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

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


A = params{1};
n = params{3};
BCs = params{4};
sensors = params{5};
S = length(sensors);

x = X(1:12); % metabolites
e = X(13:21); % enzymzes
z = X(22:end); % raw sensor optimum

F = double(ftns(x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10),x(11),x(12)));
for i=1:length(F)
  V(i) = e(i)*F(i);
end

dx(1) = 0; % The BCs are incorporated as dx1/dt = 0 and dxn/dt=0
dx(2) = V(1) - V(3);
dx(3) = 0;
dx(4) = V(2) - V(3);
dx(5) = V(3) - V(4);
dx(6) = V(4) - V(5);
dx(7) = V(5) - V(6);
dx(8) = 0;
dx(9) = V(3) - V(7);
dx(10) = V(7) - V(8);
dx(11) = V(8) - V(9);
dx(12) = 0;

% we make the DAE explicit, by converting it to an ODE
dxs = dx(sensors);
dxsensors = [zeros(n-S,1); dxs(:)];
DZ = double(dObjdxx(z(1),z(2),z(3),z(4),z(5),z(6),z(7),z(8),z(9),z(10),z(11),z(12)));
F = [DZ;zeros(S,n)];
for i=1:S
  F(n-(S-i),sensors(i))=-1; % add S rows that reflect d/dxi (x(sensors)-xi(sensors)==0)
end

dz =  -F \ dxsensors;

Eopt = double(Obj_parts(z(1),z(2),z(3),z(4),z(5),z(6),z(7),z(8),z(9),z(10),z(11),z(12)));
som=sum(Eopt);
for i=1:n-3
  de(i) = 1/100*(Eopt(i)/som - e(i));
end

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

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

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

function dy = dobjdx(x,Obj_parts,BCs);

x = x(:);
Obj = @(t,x) obj(x,Obj_parts,BCs);
dy = numjac(Obj,0,x,Obj(0,x),1e-9);

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

function [obj,y] = obj(x0,Obj_parts,BCs);

% see x0 in metsys IC for correct vector
% [BCs(1) x0(1) BCs(2) x0(2:5)' BCs(3) x0(6:8)' BCs(4)]
y=double(Obj_parts(BCs(1),x0(1),BCs(2),x0(2),x0(3),x0(4),x0(5),BCs(3),...
	x0(6),x0(7),x0(8),BCs(4)));
y = y(:)';
obj = sum(y);

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

function diffxopt = plot_optima_branched(t,y,doprint)

if nargin == 2
  doprint = 0;
end

n = 12; 
x = y(:,1:12);
e = y(:,13:21);
opt = y(:,22:end);
BCs = [1 3 8 12];
rest = setdiff([1:n],BCs);
sensors = [2,4,7,11]; 
rest = setdiff(rest,sensors);
figure(1)
  clf
  hold on
  for i=1:9, 
    p(i) = plot(t,e(:,i),'-','LineWidth',2);
  end
  grid off
  ylim = get(gca,'YLim');
  for i=1:4
    plot([1000*i,1000*i],[0,ylim(2)],'k--')
  end
  set(gca,'Box','on')
  %xlabel('time')
  %ylabel('enzyme concentration')
  hold off
  legend([p],'Enz 1','Enz 2','Enz 3','Enz 4','Enz 5','Enz 6','Enz 7','Enz 8','Enz 9');

figure(2)
  clf
  hold on
  %p1=plot(t,x(:,rest),'r-','LineWidth',2)
  %plot(t,x(:,sensors),'r--','LineWidth',2)
%  p2=plot(t,opt,'b-','LineWidth',2)
  p3=plot(t,x(:,BCs),'Color',[0 0.6 0],'LineWidth',2)
  ylim = get(gca,'YLim');
  for i=1:4
    plot([1000*i,1000*i],[0,ylim(2)],'k--')
  end
  set(gca,'Box','on')
  hold off
  grid off
  set(gca,'YLim',[0,max(y(:))*1.1]);
  %xlabel('time');
  %ylabel('concentration');
  legend([ p3(1) ],'External Metabolites')

figure(3)
  clf
  diffxopt = sum((x-opt).^2,2);
  optend = repmat(opt(end,:),length(t),1);
  diffxoptend=sum((x-optend).^2,2);
  hold on
  plot(t,diffxopt,'b-','LineWidth',2);%,t,diffxoptend)
  ylim = get(gca,'YLim');
  for i=1:4
    plot([1000*i,1000*i],[0,ylim(2)],'k--')
  end
  set(gca,'Box','on')
  grid off  
  hold off
  %xlabel('time');
  %ylabel('| metabolites - predicted optimum |^2');

if doprint
  %print_fig(1, 'pics/double_branched_enz.jpg', [0 0 5 3]);
  %print_fig(2, 'pics/double_branched_metab.jpg', [0 0 5 3]);
  %print_fig(3, 'pics/double_branched_conv.jpg', [0 0 5 3]);
end
