% 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);


data_rep = 1;
for i = 1:n

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

    data_exp{i} = [data];


end
c_fit=[];
vi=[];

for i=1:1:n

data_expi=table2array(data_exp{i});
t_exp=data_expi(:,2);
cp_exp=data_expi(:,3);
cs_exp=data_expi(:,6);
c0_s2exp=data_expi(1,9);
time=max(t_exp);
t_sim=linspace(0,time,50);

%Spline interpolation for product concentration over time (cp_fun)
step=time/5;
knots = [0:step:time];
cp_fun= slmengine(t_exp,cp_exp,'knots', knots,'concavedown','on', 'increasing', 'on' ,'minvalue', 0,'result','pp');
cs_fun= slmengine(t_exp,cs_exp,'knots', knots,'concaveup','on', 'decreasing', 'on' ,'minvalue', 0,'result','pp');

%differentiate cp_fun 

dcpdt=fnder(cp_fun,1);
dcsdt=fnder(cs_fun,1);

%define data set for rection rate vip

c_simp=ppval(cp_fun,t_sim);
c_finp=linspace(min(c_simp),max(c_simp),1000);
L=length(c_finp);
tip=[];
for j=1:L
    [value, index] = min(abs(c_simp-c_finp(j)));
    tip = [tip, t_sim(index)];
end

cp_fit=transpose([ppval(cp_fun,tip)]);

c_sims=ppval(cs_fun,t_sim);
c_fins=linspace(min(c_sims),max(c_sims),1000);
L=length(c_fins);
tis=[];
for j=1:L
    [value, index] = min(abs(c_sims-c_fins(j)));
    tis = [tis, t_sim(index)];
end
tis=fliplr(tis);
cs_fit=transpose([ppval(cs_fun,tis)]);

cs2_fit=c0_s2exp-cp_fit;

c_fit=[c_fit; cs_fit, cs2_fit, cp_fit];

vip=[ppval(dcpdt,tip(:))];
vis=[ppval(dcsdt,tis(:))];
vimean=(vip-vis)/2/data_expi(1,7);
vimeanplot=(vip-vis)/2;


vi=[vi; vimean];

 figure 
 plot(tip, vip,'*r', tis, -vis, '-b', tis, vimeanplot, '.g' )

end

%fit to kinetic equation with lsqcurvefit
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt'); 

% k1= vmax,f
% k2= KEq
% k3= Kicac
% k4= Kmcoh
% k5= Kmcac
% k6= Kmcci
% k7= Kicoh

vel_hin = @(k,c)k(1).*(c_fit(:,1).*c_fit(:,2)-c_fit(:,3)./k(2))./(k(3).*k(4)+c_fit(:,2).*k(4)+k(5).*c_fit(:,1)+c_fit(:,1).*c_fit(:,2)+k(3).*k(4)./k(6).*c_fit(:,3)+k(3).*k(4)./k(7)./k(6).*c_fit(:,3));
k0_hin = [0.0063 2.6 13.7 2.8 1.4 0.8 27.4]; %Startwerte Parameter
lb_hin = [0 0 0 0 0 0 0]; %lower Bound Parameter
ub_hin = [inf inf inf inf inf inf inf]; %upper Bound Parameter

[k_mdl_hin,SSE_hin,resid_hin,exitflag_hin,output_hin,lambda_hin,J_hin] = lsqcurvefit(vel_hin, k0_hin, c_fit, vi, lb_hin, ub_hin);


