%% Testing U-L2 model 
% by Evgenii Kovrigin, 6/2/2011
%
%
%%


close all 
clear all



%%
% Here you will find all figures
figures_folder='U_L2_testing_figures';


%% Equilbrium concentrations
% Set some meaningful parameters
Rtotal=1e-3; % Receptor concentration, M
LRratio_array=[0 : 0.02 : 1.75]; % Array of L/R 
K_a_A=1e6;  % Binding affinity constant
K_a_B=1e4;  % Dimerization constant

% Set appropriate options for the model (see model file for details)
model_numeric_solver='fminbnd'  ; 

model_numeric_options=optimset('Diagnostics','off', ...
                    'Display','off',...
                    'TolX',1e-9,...   
                    'MaxFunEvals', 1e9);                                
%%
% Important option here is TolX that sets termination tolerance on free ligand
% concentation in molar units. With our solution concentrations in 1e-3
% range TolX should be set to some 1e-9.



%%
% *Compute arrays for populations and plot*
concentrations_array=[];
for counter=1:length(LRratio_array)
    % compute
    [concentrations species_names] = equilibrium_thermodynamic_equations.U_L2_model(...
        Rtotal, LRratio_array(counter), K_a_A, K_a_B,...
        model_numeric_solver, model_numeric_options);
    % collect
    concentrations_array = [concentrations_array ; concentrations];
end

%% 
% Plot
Figure_title= 'U-L2 model';
X_range=[0 max(LRratio_array)+0.1 ]; % extend X just a bit past last point
Y_range=[ ]; % keep automatic scaling for Y
% display figure
figure_handle=equilibrium_thermodynamic_equations.plot_populations(...
    LRratio_array, concentrations_array, species_names, Figure_title, X_range, Y_range);
% save it
results_output.output_figure(figure_handle, figures_folder, 'Concentrations_plot');

%% 
% *Observations*
% The result is exactly what we expect from this model. Upon addition of a
% ligand the RL forms. As saturation is achieved - the free ligand
% accumulates and begins dimerizing when total free L concentration becomes
% appreciable

%% Test influence of very strong dimerization 
K_a_B=1e7;  % Dimerization constant

concentrations_array=[];
for counter=1:length(LRratio_array)
    % compute
    [concentrations species_names] = equilibrium_thermodynamic_equations.U_L2_model(...
        Rtotal, LRratio_array(counter), K_a_A, K_a_B,...
        model_numeric_solver, model_numeric_options);
    % collect
    concentrations_array = [concentrations_array ; concentrations];
end 
% Plot
Figure_title= 'U-L2:strong KB';
X_range=[0 max(LRratio_array)+0.1 ]; % extend X just a bit past last point
Y_range=[ ]; % keep automatic scaling for Y
% display figure
figure_handle=equilibrium_thermodynamic_equations.plot_populations(...
    LRratio_array, concentrations_array, species_names, Figure_title, X_range, Y_range);
% save it
results_output.output_figure(figure_handle, figures_folder, 'Concentrations_plot2');
%%
% By increasing dimerization constant we cause reduction in apparent
% affinity of binding. In addition, it will be difficult to achieve
% saturation because upon increasing L concentration, free L accumulates
% progressively slower. See examples below:


%% Strong dimerization, strong binding 
K_a_B=1e7;  % Dimerization constant
LRratio_array=[0 : 0.02 : 10]; % Array of L/R 

concentrations_array=[];
for counter=1:length(LRratio_array)
    % compute
    [concentrations species_names] = equilibrium_thermodynamic_equations.U_L2_model(...
        Rtotal, LRratio_array(counter), K_a_A, K_a_B,...
        model_numeric_solver, model_numeric_options);
    % collect
    concentrations_array = [concentrations_array ; concentrations];
end 
% Plot
Figure_title= 'U-L2:strong KB';
X_range=[0 max(LRratio_array)+0.1 ]; % extend X just a bit past last point
Y_range=[0 1e-3 ]; % keep automatic scaling for Y
% display figure
figure_handle=equilibrium_thermodynamic_equations.plot_populations(...
    LRratio_array, concentrations_array, species_names, Figure_title, X_range, Y_range);
% save it
results_output.output_figure(figure_handle, figures_folder, 'Concentrations_plot3');
%%
% We see that with strong binding and dimerization we readily bind initial
% amount of receptor and cannot force binding of the remaining portion
% because high concentration of free L drives formation of a dimer removing
% binding-competent L from reaction.
%
% This is in contrast with simple weak binding where we gradually achieve
% high degree of saturation upon addition of excess L (below):


%% Comparison: weak Ka for binding, no dimerization
K_a_A=1e4;  % Binding constant
K_a_B=1e-3;  % Dimerization constant

concentrations_array=[];
for counter=1:length(LRratio_array)
    % compute
    [concentrations species_names] = equilibrium_thermodynamic_equations.U_L2_model(...
        Rtotal, LRratio_array(counter), K_a_A, K_a_B,...
        model_numeric_solver, model_numeric_options);
    % collect
    concentrations_array = [concentrations_array ; concentrations];
end 
% Plot
Figure_title= 'U: weak KA';
X_range=[0 max(LRratio_array)+0.1 ]; % extend X just a bit past last point
Y_range=[0 1e-3 ]; % keep automatic scaling for Y
% display figure
figure_handle=equilibrium_thermodynamic_equations.plot_populations(...
    LRratio_array, concentrations_array, species_names, Figure_title, X_range, Y_range);
% save it
results_output.output_figure(figure_handle, figures_folder, 'Concentrations_plot4');
%%
% Here we eventually saturate the receptor


%%
% *Conclusion*
%
% The model *|equilibrium_thermodynamic_equations.U_L2_model()|* works
% well.




















