%%This program fits competition between a cooperative and a isodesmic pathway %to a series of melting curves at different concentration in a single %solvent, using 1 CD and 1 UV channel. clear all close all clc %% load and extract data %select the folder containing this script. Ensure the experimental data is %saved in a folder named 'data', Similarly, create two folders named 'figures' %and 'results' respectively. These folders will hold the files outputted by the %system. Currentpath = pwd; %pwd returns the current directory %extract filenames from 'data'-folder data = []; conc = []; datasets = {'DataPhePyTA.txt','DataPhePyTA_40.txt','SelectedDataPhePyTA.txt'} for zz = 1:3 dataset = datasets{zz}; data = dlmread(dataset); tic for z = 1:3 mkdir([dataset,'_',num2str(z)]); %% input parameters and constants R = 8.3145; %gas constant J/mole/K PL = 1; %cm nSamps = 200; %number of parameter sets Emono = 0; % Molar ellipticity monomers Eseq = 0; % Molar ellipticity sequestrated Epol_a = 63.2586/30e-6; % Molar ellipticity polymer (Taken from data at 26.67 ?C) Epol_b = 40/300e-6; % Molar ellipticity polymer (Taken from data at 26.67 ?C) H_BPTA = -86e3; S_BPTA = -155; NP_BPTA = 15.4e3; H_PhePyTA = -94e3; S_PhePyTA = -220; NP_PhePyTA = NP_BPTA; c_tot = 30e-6; model = z; %% generate fit parameters %input sample range H_max = -70e3; %estimated Gibbs free energy at 0% H_min = -100e3; S_min = -100; % m value for the aggregation at low vol% S_max = -200; %Generate sample space us latin hypercube sampling which generates evenly %space values between 0 and 1 par = lhsdesign(nSamps,2); par(:,1) = H_min + (H_max-H_min).*par(:,1); par(:,2) = S_min + (S_max-S_min).*par(:,2); %generate enthalpies from estimated Gibbs free energies and entropies labels = {'H', 'S'}; %define fitting settings Constants = [R, PL, Emono, Eseq, Epol_a, Epol_b, c_tot, H_BPTA, S_BPTA, NP_BPTA, H_PhePyTA, S_PhePyTA, NP_PhePyTA, model]; max_iter = 50; r_check = 10e99; check = []; AllResnorm = []; AllParFin= []; AllRes = []; AllJacob = cell(1,1); Freport = []; Sreport = []; lowerbound = [H_min*ones(nSamps,1), S_min.*ones(nSamps,1)]; upperbound = [H_max*ones(nSamps,1), S_max.*ones(nSamps,1)]; n = 1; j = 1; tic %% Fitting while n < nSamps Currentpar = par(n,:);%load current parameters %define optimisation options outputFun= @(x,optimValues,state,varargin)interuptFun(x,optimValues,state,varargin,r_check,max_iter); Options = optimset('MaxIter', 200, 'Display', 'iter', 'MaxFunEvals', 2000, ... 'TolFun', 1e-15, 'TolX', 1e-15, 'Algorithm', 'levenberg-marquardt', 'OutputFcn', outputFun, 'UseParallel', true); %Options = optimset('MaxIter', 200, 'Display', 'iter', 'MaxFunEvals', 2000, ... %'TolFun', 1e-15, 'TolX', 1e-15, 'Algorithm', 'levenberg-marquardt', 'UseParallel', false); lb = [];%lowerbound(n,:); ub = [];%upperbound(n,:); try [par_fin,resnorm,residual,exitflag,output,lambda,Jacobian] =... lsqnonlin(@SolvDep_2pathways_cost,Currentpar, lb,ub, Options, ... data, Constants); %check whether Aiso is within bounds % if any([par_fin(6)<0, par_fin(6)>1]) % error('Ai is out of bounds'); % end %save characteristics of current best fit if r_check > resnorm r_check = resnorm; Bestindex = n; check = [check, n]; BestJacob = Jacobian; BestResidual = residual; BestExitFlag = exitflag; end clc display(strcat('Progress: ', num2str(n/nSamps*100), '%')); display(['Resnorm: ',num2str(r_check)]) toc %save all succesfully optimised parameter sets AllResnorm = [AllResnorm; resnorm]; Sreport = [Sreport; Currentpar]; if ~isrow(par_fin) AllParFin = [AllParFin; par_fin']; else AllParFin = [AllParFin; par_fin]; end AllRes = [AllRes; residual']; AllJacob{1,n} = Jacobian; n = n + 1; j = 1; catch t = toc; if t >5400 continue end Freport = [Freport; Currentpar]; newpar = lhsdesign(1,2); newpar(:,1) = H_min + (H_max-H_min).*newpar(:,1); newpar(:,2) = S_min + (S_max-S_min).*newpar(:,2); par(n,:) = newpar; j = j + 1; if j > 1 j = 1; n = n + 1; end end end %% Analysis of fit try BestPar_Fin = AllParFin(Bestindex,:); BestRes = AllRes(Bestindex,:); catch error('No fit was found') end storage.BestPar_Fin = BestPar_Fin; BestFitIndices = find(AllResnorm == r_check); %Check for multiple distinct minima if length(BestFitIndices) > 1 all = AllParFin(BestFitIndices,:); diff = all - ones(size(all),1)*BestPar_Fin; if sum(sum(diff)) ~= 0 warning('Multiple minima detected for' + num2str(conc) + 'C'); end end %Find all good fits GoodIndices = find(AllResnorm <= 1.01*r_check); GoodFits = AllParFin(GoodIndices,:); GoodFits_log = sign(GoodFits).*log10(abs(GoodFits)); storage.GoodFits = GoodFits; %Generate some figures if length(GoodIndices) == 1 warning(['No boxplot is generated because there are no fits within 5% of the best fit']); elseif ~isempty(GoodFits_log) && size(GoodFits_log,1)>2 figure; %Create a Boxplot of goodfits boxplot(GoodFits_log, 'Labels', {labels}); title('Boxplot of logarithms of fit parameters within 1% of the best fit'); savefig([dataset,'_',num2str(z),'\2Pathways_box.fig']) close(gcf); figure; %Create matrix plot which compares the optimized parameters. [S,AX,BigAx,H,HAx] = plotmatrix(GoodFits); title(BigAx,'comparisons of fit parameters') for j = 1:size(GoodFits,2) AX(j,1).YLabel.String = labels(j); AX(size(GoodFits,2),j).XLabel.String = labels(j); end savefig([dataset,'_',num2str(z),'\2Pathways_cor.fig']); close(gcf); else disp(['Only one distinct fit found for ', num2str(conc),'C, no boxplot generated']); end %Extract best fit parameters H_fit = BestPar_Fin(1); S_fit = BestPar_Fin(2); %Calculate bestfit curves uniquefracs = unique(data(:,1)); for i = 1:length(uniquefracs) frac = uniquefracs(i); simT = linspace(0,100,100); simdat = [ones(1,100)*frac;simT]'; [CalcCD, CalcA, CalcB, CalcP_a, CalcP_b, CalcSeq_a, CalcSeq_b] = SolvDep_2pathways_simulator(BestPar_Fin, Constants,simdat); figure(1); hold on plot(simT, CalcCD, '-'); xlabel('Temperature (C)'); ylabel('Absorbance (a.u.)'); figure(2); hold on plot(simT, CalcA,'r',simT, CalcB,'m',simT,CalcP_a,'g',simT,CalcP_b,'y',simT,CalcSeq_a,'b',simT,CalcSeq_a,simT,(CalcA+CalcP_a+CalcSeq_a),'k',simT,CalcB+CalcP_b+CalcSeq_b); legend('MonomerA','MonomerB', 'PolymerA', 'PolymerB','SequestratedA', 'SequestratedB','TotalA','TotalB') xlabel('Temperature (?C)'); ylabel('Concentration'); %save the simulated signal FileID = fopen(strcat([dataset,'_',num2str(z),'\2Pathways_signal_c=_', num2str(frac), '.txt']),'wt'); txt = sprintf('Fraction \t Absorbance signal'); fprintf(FileID, '%s\t\r\n' ,txt); for j = 1:length(simT) DATA = [simT(j),CalcCD(j)]; fprintf(FileID, '%5f\t%5f\n', DATA); end %species FileID = fopen(strcat([dataset,'_',num2str(z),'\PhePyTA_sequestration_ratio=', num2str(frac), '.txt']),'wt'); txt = sprintf('Temperature \t MonomerA \t MonomerB \t PolymerA \t PolymerB \t SequestratedA \t SequestratedB \t TotalA \t TotalB'); fprintf(FileID, '%s\t\r\n', txt); for j = 1:length(simT) DATA = [simT(j), CalcA(j),CalcB(j), CalcP_a(j), CalcP_b(j), CalcSeq_a(j),CalcSeq_b(j), (CalcA(j) + CalcP_a(j) + CalcSeq_a(j)),(CalcB(j) + CalcP_b(j) + CalcSeq_b(j))]; fprintf(FileID, '%5e\t%5e\t%5e\t%5e\t%5e%5e\t%5e\t%5e\t\r\n', DATA); end end figure(1) plot(data(:,2),data(:,3),'o','MarkerSize',2) savefig(strcat([dataset,'_',num2str(z),'\2Pathways_Fit_c=_', num2str(frac), '.fig'])); figure(2) savefig(strcat([dataset,'_',num2str(z),'\2Pathways_Speciation_c=_', num2str(frac), '.fig'])); %perform statistics on the fit parameters df = -length(BestPar_Fin); %number of degrees of freedom for i = 1:size(data,2) df = df + length(data(:,1)); end alpha = 0.32; %one standard deviation try ci = nlparci(BestPar_Fin, BestResidual, BestJacob, alpha); %calculates the 2*alpha% confidence interval, returns upper and lower boundaries per parameter tsd = tinv(1-alpha./2,df); x_sd = (ci(:,2)-ci(:,1))./(2*tsd); % calculates the standard deviation for each parameter catch x_sd = zeros(size(BestPar_Fin)); end storage.x_sd = x_sd; storage.Jacob = BestJacob; storage.Constants = Constants; %an alternative method of determining standard deviations. cov = r_check./df.*inv(transpose(BestJacob)*BestJacob); %generates the covariant matrix xsd = sqrt(diag(cov)); %generates alternative standard deviations corr = zeros(size(cov)); for i = 1:size(cov,1) %generates correlation matrix for j = 1:size(cov,2) corr(i,j) = cov(i,j)./(sqrt(cov(i,i).*cov(j,j))); end end storage.cov = cov; storage.corr = corr; storage.xsd = xsd; save('storage_iso.mat', '-struct','storage'); %export parameters FileID = fopen(strcat([dataset,'_',num2str(z),'\2Pathways_parameters.txt']), 'wt'); txt = sprintf('H= %5e\t+/-\t%5e\r\nS= %5e\t+/-\t%5e\r\n\nresnorm= %5e', ... H_fit, x_sd(1),S_fit,x_sd(2),r_check); fprintf(FileID,'%s\t\r\n', txt); fclose(FileID); close all end end