function [t,x,x0,Y0,Finv] = daes_extra_param;

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

format compact
n = 7; % number of metabolites
% x1 -> x2 -> x3 -> x4 -> x5 -> x6 -> x7
%				 ^
%				 |
%				 K3
nenz = 6; % number of reactions/enzymes
% parameters
A = [
3 1 3 2 1 2 3 1 3 2; 
5 2 4 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
K = 40; % the fixed K3 parameter, to which the system needs to be robust as well
sensors = [2,3,6]; % 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} = nenz;
params{8} = K;

loadfunctions=0;
if loadfunctions
  disp('... loading Symbolic functions ...')
  load('Symbolic_functions/daes_extra_param_functions.mat');
else
  %define the object function symbolically, and compute its symbolic Jacobian, 
  % which yields the optimum equations.
  syms x1 x2 x3 x4 x5 x6 x7 K3
  vars=[x1 x2 x3 x4 x5 x6 x7 K3];
  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)/(  K3   + 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);

  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);
  DOBJDXX   = symfun(jacobian(DOBJDX,vars),vars);
  % now make these normal matlab functions: much faster evaluation
  ftns      = matlabFunction(FTNS);
  Obj_parts = matlabFunction(OBJ_PARTS);
  Objsym    = matlabFunction(OBJSYM);
  dObjdx    = matlabFunction(DOBJDX);
  dObjdxx   = matlabFunction(DOBJDXX);

  save('Symbolic_functions/daes_extra_param_functions.mat', 'vars','ftns', 'Obj_parts','Objsym','dObjdx','dObjdxx');
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
% in this case, we just take some arbitrary optimum, using fake BCs, which 
% are twice the original ones. they could be anything (realistic)
BCs0 = [9 3];
K0 = 20;
Obj = @(x0) sum(Obj_parts(BCs0(1),x0(1),x0(2),x0(3),x0(4),x0(5),BCs0(2),K0));
%obj(x,Obj_parts,BCs0,K0);
z0 = [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) K0];
% 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),zopt(8));

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

X0 = fullxvector(x0,BCs);
x0 = [X0(:); e(:); zopt(:)]; 
x0 = x0(:)'

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

plop=1;
if plop
BCs=[5,1];
params{4}=BCs;
x0=x(end,:);
x0(1)=BCs(1);
x0(n)=BCs(2);
[t1,x1] = ode15s(ODE,[t(end),t(end)+5000],x0,options);
t=[t;t1];
x=[x;x1];

BCs=[5,3];
params{4}=BCs;
x0=x(end,:);
x0(1)=BCs(1);
x0(n)=BCs(2);
[t1,x1] = ode15s(ODE,[t(end),t(end)+5000],x0,options);
t=[t;t1];
x=[x;x1];

BCs=[10,3];
params{4}=BCs;
x0=x(end,:);
x0(1)=BCs(1);
x0(n)=BCs(2);
[t1,x1] = ode15s(ODE,[t(end),t(end)+5000],x0,options);
t=[t;t1];
x=[x;x1];

K=100;
params{8}=K;
x0=x(end,:);
x0(1)=BCs(1);
x0(n)=BCs(2);
ODE = @(t,x) ode(t,x,ftns,Obj_parts,dObjdxx,params);
[t1,x1] = ode15s(ODE,[t(end),t(end)+5000],x0,options);
t=[t;t1];
x=[x;x1];
end

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


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

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

e = params{2};
N = params{6};
K = params{8};

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

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

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

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

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

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),K); 
V = e(:).*F(:);

dx = N*V(:);

% we make the DAE explicit, by converting it to an ODE
dxsensors = [zeros(n-nsens+1,1); dx(sensors)];
DZ = dObjdxx(z(1),z(2),z(3),z(4),z(5),z(6),z(7),z(8));
[DZr,DZc] = size(DZ);
F = [DZ;zeros(nsens,n+1)]; % n+1 because of the extra K 
for i=1:nsens,
  F(DZr+i,sensors(i))=-1; 
end
size(dxsensors);
%test if these sensors give an invertible F
%[U,S,V]=svd(F);
%mmin = min(abs(diag(S)));
%if mmin < 1e-4
%  disp(' detecting a singular value close to zero. This choice of sensors %seems to give a noninvertible D_xi G')
%end
F;
dz =  - F \ dxsensors;
%dz(end)=0;

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

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

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

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

function plot_optima_branched(t,y,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};
K = params{8};

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


figure(1)
  clf
  hold on
  nenz
  for i=1:nenz
    p(i) = plot(t,e(:,i),'LineWidth',2);
  end
  legend([p],'Enzyme 1', 'Enzyme 2','Enzyme 3','Enzyme 4','Enzyme 5','Enzyme 6')
  grid off
  ylim = get(gca,'YLim')*1.1;
  set(gca,'YLim',[0,ylim(2)]);
  plot([10000,10000],[0,ylim(2)],'k--')
  plot([15000,15000],[0,ylim(2)],'k--')
  plot([20000,20000],[0,ylim(2)],'k--')
  plot([25000,25000],[0,ylim(2)],'k--')
  hold off
  %xlabel('time')
  %ylabel('enzyme concentration')
%  title ('enzymes')
figure(2)
  clf
  hold on
  %p2=plot(t,opt,'b-','LineWidth',2);
  p1=plot(t,x(:,rest(2:end)),'r-','LineWidth',2);
  plot(t,x(:,rest(1)),'r--','LineWidth',2);
  p3=plot(t,x(:,BCs),'Color',[0 0.6 0],'LineWidth',2);
  grid off
%  p1=plot(t,x(:,rest),'r-o','LineWidth',2);
%  p2=plot(t,opt(:,1:end-1),'b-','LineWidth',2);
%  p4=plot(t,[x(:,BCs) ],'Color',[0 0.6 0],'LineWidth',2);
  
  grid off
  ylim = get(gca,'YLim')*1.1;
  set(gca,'YLim',[0,ylim(2)]);
  plot([10000,10000],[0,ylim(2)],'k--')
  plot([15000,15000],[0,ylim(2)],'k--')
  plot([20000,20000],[0,ylim(2)],'k--')
  plot([25000,25000],[0,ylim(2)],'k--')
  hold off
  %xlabel('time');
  %ylabel('metabolite/sensor concentration');
  legend([p1(1) p3(1)],'Internal metabolites', 'External Metabolites')

figure(3)
  clf
  hold on
  T = linspace(0,t(end),30);
  T = [T(1:25) T(25:30)];
  p3=plot(t,opt(:,end),'Color',[204 119 34]/256,'LineWidth',2);
  Ks = [40*ones(1,25) 100*ones(1,6)];
  p4 = plot(T,Ks,'Color',[0 0.6 0],'LineWidth',2);
  ylim = get(gca,'YLim')*1.1;
  set(gca,'YLim',ylim);
  plot([10000,10000],ylim,'k--')
  plot([15000,15000],ylim,'k--')
  plot([20000,20000],ylim,'k--')
  plot([25000,25000],ylim,'k--')
  %xlabel('time');
  %ylabel('K_3')
  legend([p4(1) p3(1)],'true K_3 value','sensor K_3 value')
  hold off

 

if doprint
  print_fig(1, 'pics/linear_chain_kinetic_parameter_dynamics_enz.jpg', [0 0 6 4]);
  print_fig(2, 'pics/linear_chain_kinetic_parameter_dynamics_metab.jpg', [0 0 6 4]);
  print_fig(3, 'pics/linear_chain_kinetic_parameter_dynamics_K3.jpg', [0 0 6 4]);
 
end

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

function recompute_functions(A,K3);
  %define the object function symbolically, and compute its symbolic Jacobian, 
  % which yields the optimum equations.
  syms x1 x2 x3 x4 x5 x6 x7 K3
  vars=[x1 x2 x3 x4 x5 x6 x7 K3];
  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)/(  K3   + 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);

  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);
  DOBJDXX   = symfun(jacobian(DOBJDX,vars),vars);
  % now make these normal matlab functions: much faster evaluation
  ftns      = matlabFunction(FTNS);
  Obj_parts = matlabFunction(OBJ_PARTS);
  Objsym    = matlabFunction(OBJSYM);
  dObjdx    = matlabFunction(DOBJDX);
  dObjdxx   = matlabFunction(DOBJDXX);
