function [p, K, R2_vector, v_vector]=rate_matrix(beta, x)
%------------------- LineShapeKin RATE_MATRIX -----------------------------
%-- version 3.1 
%
% 3/17/08
% Written by Evgenii Kovrigin, Medical College of Wisconsin
%
% This function builds a rate (K) matrix for the spin exchange between the
% sites corresponding to the equation:  
%      P+L=PL, Kd=k2/k1, where k1=kon*[L], k2=Koff
%
%   NOTES:  this is a version with fixed preset Kd (only koff and scaling factor
%   are variable parameters.
%
%
% 3/29/08 - Here I introduced an alternative way of dealing with complex
%                  equilibria: numeric minimization to obtain equilibrium
%                  concentrations of components [L] needed for building of
%                  K-matrix. 
%
%                   THIS IS SIMPLE BINDING EXAMPLE. You may switch it from
%                   numeric to analytical solution.  
%
%                   NOTE: This mode is great for easy introducing new
%                   models of exchange but is approx. 20x slower than
%                   analytical solution.
%
%                   It is advisable to make an effort to obtain analytical solution 
%                   whenever possible. This file serves as example of numerical 
%                   solution setup for a new model but I will turn 'analytical' mode 
%                   on for practical purposes. 

global B_vector Model  X_columns 

mode='analytical';
%mode='numeric';

  
%--- extract fixed (preset) values
%- use fixed preset Kd
Kd = Model.Kd;   

%-- extract variable parameters
koff=beta(B_vector.Koff);




%-- extract experimental parameters:

Ptot = x(X_columns.Ptotal);    
Ltot = Ptot*x(X_columns.LPratio);  


 
%% Analytical solution for [L] out of 
%
%   Kd = [P][L]/[PL]


switch mode
    
    %---- traditional explicit calculation of [L]
    case 'analytical'
            a=1;
            b=-(Ptot+Ltot+Kd);
            c=Ptot*Ltot;
            sola=-b/(2*a);
            solb=sqrt(b^2-(4*a*c))/(2*a);
            sol1=[(sola+solb) (sola-solb)];
            Peq=max(Ptot-sol1); %-- get positive solution
            PLeq=Ptot-Peq;
            Leq=Ltot-PLeq;
            
            
            
            
    %--- minimization to obtain [L]
    case 'numeric'

            

            Leq=0;

            if Ltot>0   %-- to skip 0 concentration point
            %   obtain equilibrium concentration of the ligand by numeric
            %   minimization.

            % --STARTING CONDITIONS. It is essential for numeric stability
            % of the solution to start from just below full concentration.
            % Sometimes it may be unstable - try using different starting
            % conditions.
           
            Leq0=Ltot*0.9;
            params=[Leq0];
            
            % -- I formally define X. This is just one of the constants in
            % the equation. The other constants are supplied directly to
            % the function below.
            x_local=Kd;
            
            % this is Y-array 
            y_local=Ptot;

             %-- physically meaningful bounds
             lower_bounds=[0];  
              upper_bounds=[Ltot];  

               %-- default fitting options. Sometimes TolFun and TolX need
               %    to be adjusted for better convergence. 
              options=optimset('Diagnostics','off', ...
                  'Display','off',...
                  'LargeScale','on',...
                  'LevenbergMarquardt', 'on',...
                  'TolFun',1e-18,...
                  'TolX',1e-18,...
                 'MaxFunEvals', 1e9);
             
             
             %  --- we define the target function using 'anonymous' synthax
             %  and pass additional parameters to it here
             %myfun = @(B,X) ( lnKd1fun(B, [Ptot,X]) );
             myfun = @(B,X) ( Model_AB(B, [X,Ltot]) );
             
             % --- run least-squares minimization protocol to extract
             % equilibrium concentrations we need for K-matrix ('params'
             % vector)
              [Leq, resnorm, residual, exitflag, output,lambda, Jacobian]=...
                  lsqcurvefit(myfun, params, x_local, y_local,...
                  lower_bounds,upper_bounds,options);

  
            end
            
            %-- compute concentrations of other other species
            Peq= Ptot-Ltot+Leq;
            PLeq=Ptot-Peq;

end



% Equilibrium mole fractions of species containting P 
p=[Peq/Ptot,PLeq/Ptot];

%-- relationships of true constants
kon=koff/Kd;

% K matrix elements
kab=kon*Leq;
kba=koff;

K=[-kab,kba;kab,-kba];



%-- prepare relaxation rates and frequencies for the individual species
 R2_vector=[ x(X_columns.R2A),  x(X_columns.R2B) ];  
 v_vector=[ x(X_columns.wA),    x(X_columns.wB)   ];  
 
 
 
