%%% --- Biophysics of Macromolecules   --- %%%
%%% --- LMU - Sommersemester 2020      --- %%%
%%% --- Prof. Lipfert                  --- %%%
%%%
%%% Solutions to 
%%% Problem set 3, Problem 3
%%% Ligand receptor binding, free vs. bound ligand

clc;clf;clear;

%%% Read the data from an excel file
%%% Make sure the xlsx file is in the same directory
data = xlsread('Maeda_2000_Figure1A_DATA.xlsx');
c_core = data(:,1);      %%% Core enzyme concentration, nM
p_bound = data(:,2)/100; %%% Fraction bound --- convert from %-bound



%%% Fit simple binding model to the data
%%% This pretend that free ligand = total ligand
[Kd res exitflag] = fminsearch(@(Kd) norm(p_bound - (1./(1+Kd./c_core))), [0.4]);
xfit = 0:0.01:1;
yfit = xfit./(xfit+Kd);


%%% Fit binding model that takes into account the total ligand and receptor concentration
%%% Kdt is the Kd constant determined for the model that takes into account
%%% total concentrations
Rtot = 0.4; %%% in nM
Ltot = c_core; %%% in nM

[Kdt res exitflag] = fminsearch(@(Kdt) norm(p_bound -  ( ((Rtot+Ltot+Kdt)-sqrt((Rtot+Ltot+Kdt).^2 - 4*Rtot*Ltot))./(2*Rtot) ) ), [0.2]);
yfitt = ( ((Rtot+xfit+Kdt)-sqrt((Rtot+xfit+Kdt).^2 - 4*Rtot*xfit))./(2*Rtot) );



%%% Plot the results
figure(1); clf; hold on; box on;
plot(xfit, yfit, 'b--', 'linewidth', 2)
plot(xfit, yfitt, 'b-', 'linewidth', 2)
plot(c_core, p_bound, 'ks', 'linewidth', 2, 'markersize', 20)

legend('L_{tot} = L_{free} model','L_{tot}, R_{tot} model', 'location', 'southeast')

set(gca,'LineWidth', 1,'FontSize', 25, 'FontWeight', 'bold', 'TickLength',[0.02 0.02])
ylabel('Fraction bound sigma');
xlabel('Core enzyme (nM)');
axis([-0.1 1.1 -0.1 1.1])

%%% Show the fitted Kd constants in the title of the graph
title(['Simple Kd = '  num2str(Kd,2) ' nM; Correct Kd = ' num2str(Kdt,2) ' nM' ])

%%% Print the results to eps
print(gcf,'-r600','-depsc', 'Ltot_BindingModels.eps');