% This optional module, "Fit_TwoKD," is used as a computing and plotting space.
% It calls function "Fit_TwoKD_Func" in order to return fractional saturations, 
% KDs, and Hill coefficients for a model that assumes two transition points,. 

%% INPUT DATA BELOW (SEE ALSO SAMPLE DATA)

% Ligand concentrations (L) should all be in uM.
% Input normalized fluorescence signals (S) rather than raw data. To avoid divide
% by zero errors, put 0.9999 and 0.0001 in place of 1 and 0, respectively.

% This function is sensitive to having data that actually reflects a sample with
% two distinct transitions. The example data below satisfies this requirement - try
% running the program once and looking at the resulting plot before testing your data.

L= [0.0001 0.01 0.05 0.1 0.2 0.4 0.6 0.8 1 2.5 10 100];  
S = [0.98 0.95 0.90 0.80 0.65 0.60 0.62 0.45 0.30 0.15 0.10 0.12];

% Below we initialize an array of ligand concentrations.

L_space = logspace(-4, 4, 10000);    

% Binding fractions are calculated from fluorescence below. 

Y = 1 - S;

% The fitting parameter array is initialized. The fractional saturation at
% transition (f), KD, and n *at each transition* are the six parameters.

C= [0.5 0.1 1 0.2 0.1 1];

% The fitting function is called.

[fKDn, r] = nlinfit(L, Y, @Fit_TwoKd_Func, C);

% Output f, KD and n values so they show up in the command window.

f1 = fKDn(1)
KD1 = fKDn(2)
n1 = fKDn(3)
f2 = fKDn(4)
KD2 = fKDn(5)
n2 = fKDn(6)

% Using the same model as in "Fit_TwoKD_Func," we prepare a vector of 
% values for the binding fraction, using the chosen ligand-space. 

fitY = f1.*(L_space.^n1)./(L_space.^n1 + KD1.^n1) + f2.*(L_space.^n2)./(L_space.^n2 + KD2.^n2);

% Finally, we plot the data.

figure
semilogx(L, Y, 'x')
hold on
semilogx(L_space, fitY)
xlabel('log [L] (\mu M)')
ylabel('binding Y')
legend('data', 'fit')
title('Part 4: Two transitions')