% This Code is part of the publication 
% Development of an experimental setup and a validated model for a highly integrated biocatalytic process

%Before use: add experimental data from folder
%Experimental_data_SuppInfo/Lipase kinetic

close all
clear
clc

%Retrieve experimental Data from Excel File 

data_num_in = inputdlg('Enter Number of data sets', 'Data sets input', [1 40]);
n_str = cell2mat(data_num_in);
n = str2num(n_str);
data_exp=cell(n,1);
c0_Alc_all=[];
c0_Ac_all=[];
c0_Est_all=[];
E_all=[];
c_exp_all=[];
c_p_sim_err=[];
v_sim=[];

for i = 1:n

    s1 = "V"; 
    s3 = ".xlsx";
    s2 = num2str(i);
    s = strcat(s1,s2,s3);
    data = readtable(s);

    data_exp{i} = [data];


end



for i = 1:n

    data_expi=table2array(data_exp{i});
    t_exp=data_expi(:,2);
    t1 = max(t_exp);
    cp_exp=data_expi(:,3);
    cs_exp=data_expi(:,6);
    E=data_expi(1,7);
    c0_Alc=data_expi(1,8);
    c0_Ac=data_expi(1,9);
    c0_Est=data_expi(1,10);
    k = [3.558732358114269e-04,25.893017684977547,1.141050521837446,4.519728908114671,10.366590558142567,0.156262707413453,4.063469504805376];
    v= @(t,c)E*k(1).*((c0_Ac-c).*(c0_Alc-c)-(c)./k(2))./(k(3).*k(4)+(c0_Ac-c).*k(4)+k(5).*(c0_Alc-c)+(c0_Alc-c).*(c0_Ac-c)+k(3).*k(4)./k(6).*(c)+k(3).*k(4)./k(7)./k(6).*(c));
    vp= @(c)E*k(1).*((c0_Ac-c).*(c0_Alc-c)-(c)./k(2))./(k(3).*k(4)+(c0_Ac-c).*k(4)+k(5).*(c0_Alc-c)+(c0_Alc-c).*(c0_Ac-c)+k(3).*k(4)./k(6).*(c)+k(3).*k(4)./k(7)./k(6).*(c));
    vp_sim = vp(cp_exp);
    v_sim=[v_sim; vp_sim];
  
tspan=linspace(0, max(t1), 1000);

[t,c] = ode45(v,tspan, c0_Est);

figure;
plot(tspan, c, '*y', t_exp, cp_exp, '.b')
hold on;

t_sim_err = [];
L=length(cp_exp);

for m=1:L
    [value, index] = min(abs(tspan-t_exp(m)));
    t_sim_err = [t_sim_err, tspan(index)];
    c_p_sim_err = [c_p_sim_err, c(index)];
end

c_exp_all=[c_exp_all; cp_exp];
c0_Alc_all = [c0_Alc_all, repmat(c0_Alc, 1, length(cp_exp))];
c0_Ac_all = [c0_Ac_all, repmat(c0_Ac, 1, length(cp_exp))];
E_all = [E_all, repmat(E, 1, length(cp_exp))];

end

% Compute the error
SSE = sum(((c_p_sim_err-transpose(c_exp_all))).^2);
Res1 = abs((c_p_sim_err(:)-c_exp_all(:))./c_exp_all(:));
Res1(isnan(Res1))=0;
Res=v_sim(:).*Res1;

sigma_2 = SSE/(length(c_exp_all)-length(k)); % number of data point minus paramter

% Create Jacobimatrix
% Definition of the partial derivatives and columns of Jacobian matrix

syms c k1 k2 k3 k4 k5 k6 k7 c0_Ac1 c0_Alc1 Em

v = Em*k1.*((c0_Ac1-c).*(c0_Alc1-c)-(c)./k2)./(k3.*k4+(c0_Ac1-c).*k4+k5.*(c0_Alc1-c)+(c0_Alc1-c).*(c0_Ac1-c)+k3.*k4./k6.*(c)+k3.*k4./k7./k6.*(c));

J_sym = jacobian(v, [k1, k2, k3, k4, k5, k6, k7]);  % Symbolic Jacobian matrix

J = zeros(length(c_exp_all), 4);  % Initialize Jacobian matrix

for i = 1:length(c_exp_all)
    c_val = c_exp_all(i);
    c0Ac_val = c0_Ac_all(i);
    c0Alc_val = c0_Alc_all(i);
    E_val = E_all(i);
    J(i, 1) = double(subs(J_sym(1), [c, k1, k2, k3, k4, k5, k6, k7, c0_Ac1, c0_Alc1, Em], [c_val, k(1), k(2), k(3), k(4), k(5), k(6), k(7), c0Ac_val, c0Alc_val, E_val]));
    J(i, 2) = double(subs(J_sym(2), [c, k1, k2, k3, k4, k5, k6, k7, c0_Ac1, c0_Alc1, Em], [c_val, k(1), k(2), k(3), k(4), k(5), k(6), k(7), c0Ac_val, c0Alc_val, E_val]));
    J(i, 3) = double(subs(J_sym(3), [c, k1, k2, k3, k4, k5, k6, k7, c0_Ac1, c0_Alc1, Em], [c_val, k(1), k(2), k(3), k(4), k(5), k(6), k(7), c0Ac_val, c0Alc_val, E_val]));
    J(i, 4) = double(subs(J_sym(4), [c, k1, k2, k3, k4, k5, k6, k7, c0_Ac1, c0_Alc1, Em], [c_val, k(1), k(2), k(3), k(4), k(5), k(6), k(7), c0Ac_val, c0Alc_val, E_val]));
    J(i, 5) = double(subs(J_sym(5), [c, k1, k2, k3, k4, k5, k6, k7, c0_Ac1, c0_Alc1, Em], [c_val, k(1), k(2), k(3), k(4), k(5), k(6), k(7), c0Ac_val, c0Alc_val, E_val]));
    J(i, 6) = double(subs(J_sym(6), [c, k1, k2, k3, k4, k5, k6, k7, c0_Ac1, c0_Alc1, Em], [c_val, k(1), k(2), k(3), k(4), k(5), k(6), k(7), c0Ac_val, c0Alc_val, E_val]));
    J(i, 7) = double(subs(J_sym(7), [c, k1, k2, k3, k4, k5, k6, k7, c0_Ac1, c0_Alc1, Em], [c_val, k(1), k(2), k(3), k(4), k(5), k(6), k(7), c0Ac_val, c0Alc_val, E_val]));

end

for i=1:length(k)
    % Extract confidence intervals
    CI_(i,:) = nlparci(k(i),Res,'jacobian',J(:,i),'alpha',0.32);
    % Calcualte individual confidence levels (in %)
    CIp_(i,:) = (CI_(i,:)-k(i))./k(i).*100;
    CI__(i,:)=(CI_(i,:)-k(i));
    SE(i,:) = nlparci(k (i),Res,'jacobian',J(:,i)); %,'alpha',0.32);
    % Calcualte individual confidence levels (in %)
    SEp(i,:) = (SE(i,:)-k(i))./k(i).*100;
end
