%plotAndCalculateKd_MultipleClones
close all
clc

%This code will fit your titration curves for your scFv binding to lysozyme as measured via flow cytometry. 
%Inputs required: concentration of lysozyme and corresponding fluoresence values for controls and test clones at those concentrations
%Output: plot with titration curves, Kd values

%% Input:

%%%% Input ligand concentration values used for controls and samples %%%%
lysozymeConc = [0.316, 1, 3.16, 10, 31.6, 100, 316, 1000];

%%%% Input raw fluorescence values for controls and clones %%%%
control1 = [72, 255, 831, 3659, 6090, 7756, 9377, 9728];
clone1 = [167, 362, 1098, 2810, 3685, 3828, 4019, 4478];

%%%% Input raw fluorescence values for additonal clones as necessary
clone2 = [33, 101, 310, 521, 753, 689, 817, 705];
clone3 = [20, 34, 39, 46, 59, 198, 274, 327];

%% Solve
lb = [0 0 -5000];
opt = optimoptions('fmincon','Display','off','FunctionTolerance',10^-10);
constantsGuess = [0,0,0];
constantControl = fmincon(@(constants)solveKd(lysozymeConc, control1, constants),constantsGuess, [],[],[],[],lb,[],[],opt);
constantClone1 = fmincon(@(constants)solveKd(lysozymeConc, clone1, constants),constantsGuess, [],[],[],[],lb,[],[],opt);

%% Normalize and plot
X = logspace(log10(min(lysozymeConc)),log10(max(lysozymeConc)),100);
str3 = ' nM';

%plot control
YSolve = constantControl(1).*X./(constantControl(2)+X) + constantControl(3);
normSolveControl1 = ((YSolve-constantControl(3))./constantControl(1));
normControl1 = ((control1-constantControl(3))./constantControl(1));
Control_Kd = round(constantControl(2),2)
str1 = 'Control K_{D} =';
str2 = num2str(Control_Kd);
controlStr = strcat(str1,32,str2,str3);

%plot clone1
YSolve = constantClone1(1).*X./(constantClone1(2)+X) + constantClone1(3);
normSolveClone1 = ((YSolve-constantClone1(3))./constantClone1(1));
normClone1 = ((clone1-constantClone1(3))./constantClone1(1));
Clone1_Kd = round(constantClone1(2),2)
str1 = '\color{red}Clone1 K_{D} =';
str2 = num2str(Clone1_Kd);
Clone1Str = strcat(str1,32,str2,str3);

%Plot Figure 1 (Control vs. Clone 1)
plot(lysozymeConc,normControl1,'ko',X,normSolveControl1,'k-',...
    lysozymeConc,normClone1,'ro',X,normSolveClone1,'r-','LineWidth',2)
str = {controlStr,Clone1Str};
annotation('textbox',[.15 .6 .3 .3],'String',str,'FitBoxToText','on');

%plot specs
set(gca, 'XScale', 'log')
set(gca, 'LineWidth', 1.5)
ylim([0 1.1])
title('scFv Titration Curves','FontSize',16,'FontWeight','bold')
ylabel('Normalized Fluoresence','FontSize',14)
xlabel('Lysozyme Concentration [nM]','FontSize',14)

%% Solve and Plot Additonal Clones
constantClone2 = fmincon(@(constants)solveKd(lysozymeConc, clone2, constants),constantsGuess, [],[],[],[],lb,[],[],opt);
constantClone3 = fmincon(@(constants)solveKd(lysozymeConc, clone3, constants),constantsGuess, [],[],[],[],lb,[],[],opt);

%plot clone2
YSolve = constantClone2(1).*X./(constantClone2(2)+X) + constantClone2(3);
normSolveClone2 = ((YSolve-constantClone2(3))./constantClone2(1));
normClone2 = ((clone2-constantClone2(3))./constantClone2(1));
Clone2_Kd = round(constantClone2(2),2)
str1 = '\color{blue}Clone2 K_{D} =';
str2 = num2str(Clone2_Kd);
Clone2Str = strcat(str1,32,str2,str3);

%plot clone3
YSolve = constantClone3(1).*X./(constantClone3(2)+X) + constantClone3(3);
normSolveClone3 = ((YSolve-constantClone3(3))./constantClone3(1));
normClone3 = ((clone3-constantClone3(3))./constantClone3(1));
Clone3_Kd = round(constantClone3(2),2)
str1 = '\color{magenta}Clone3 K_{D} =';
str2 = num2str(Clone3_Kd);
Clone3Str = strcat(str1,32,str2,str3);

%Plot Figure 2 (Control vs. all Clones)
figure;
plot(lysozymeConc,normControl1,'ko',X,normSolveControl1,'k-',...
    lysozymeConc,normClone1,'ro',X,normSolveClone1,'r-',...
    lysozymeConc,normClone2,'bo',X,normSolveClone2,'b-',...
    lysozymeConc,normClone3,'mo',X,normSolveClone3,'m-','LineWidth',2)
str = {controlStr,Clone1Str,Clone2Str,Clone3Str};
annotation('textbox',[.15 .6 .3 .3],'String',str,'FitBoxToText','on');

%plot specs
set(gca, 'XScale', 'log')
set(gca, 'LineWidth', 1.5)
ylim([0 1.1])
title('scFv Titration Curves','FontSize',16,'FontWeight','bold')
ylabel('Normalized Fluoresence','FontSize',14)
xlabel('Lysozyme Concentration [nM]','FontSize',14)

%% Least Sum Squares Function

function minFunc = solveKd(X,Y,constants)
    Bmax = constants(1);
    Kd = constants(2);
    Background = constants(3);

    residuals = Y - (Bmax.*X./(Kd+X) + Background);
    minFunc = sum(residuals.^2);
end
