Lac operon MWC linked equilibrium solver

by

 

01 Mar 2013 (Updated )

Finds assumption free solutions to two Monod, Wyman and Changeux (MWC) models of the lac operon.

linked_equilibrium_REO(params)
function [le m] = linked_equilibrium_REO(params)
    %        Kre         Kre
    %   RO  -----  ROE  ------ ROEE
    %    |          |            |
    %    | Kro      |            |
    %    |   Kre    |     Kre    |
    %    R  -----  RE   ------- REE
    %    |          |            |
    %    | Krr*     |            |
    %    |   Kr*e   |    Kr*e    |
    %   R*  -----  R*E  ------ R*EE
    %    |          |            |
    %    | Kr*o     |            |
    %    |   Kr*e   |     Kr*e   |
    %   R*O ----- R*EO ------- R*EEO
    %
    %
    %   input: param.Kro, param.Krsto
    %          param.Kre, param.Krste
    %          param.Krrst
    %          param.Otot, param.Etot
    %          param.Rtot
    %
    %   output: le (concentrations of all species structure) 
    %           minvalue (optional)
    %
    %   [le (m)] = linked_equilibrium_REO(param);
    %
    %   everything should be in the same units, nM, uM or whatever
    %
    %   Requires: DocPolynom library - available from Matlab file exchange
    %             patternsearch - part of the Global Optimization Toolbox
    %
    %   Matthew A. Sochor - Dec 2011
    
    Krrst = params.Krrst;
    Kre = params.Kre;
    Krste = params.Krste;
    Kro = params.Kro;
    Krsto = params.Krsto;
    Rtot = params.Rtot;
    Otot = params.Otot;
    Etot = params.Etot;
    if Rtot == 0
        le.E = Etot; 
        le.R = 0;
        le.O = Otot;
        le.Et = Etot;
        le.Rt = Rtot;
        le.Ot = Otot;
        le.RE = 0;
        le.RstE = 0;
        le.REE = 0;
        le.RstEE = 0;
        le.RO = 0;
        le.ROE = 0;
        le.ROEE = 0;
        le.RstO = 0;
        le.RstOE = 0;
        le.RstOEE = 0;
        le.params = params;
        m=-inf;
        return
    end
    % Make constants
    alpha1 = 1 + Krrst;
    beta1 = 2*Kre + 2*Krrst*Krste;
    gamma1 = (Kre^2) + Krrst*(Krste^2);
    gamma2 = 2*Kro*Kre;
    delta1 = Kro*(Kre^2);
    beta2 = Krrst*Krsto;
    gamma3 = 2*Krrst*Krsto*Krste;
    delta2 = Krrst*Krsto*(Krste^2);
    
    % Build polynomials.  Note these differ from the methods as they are
    % the correct polynomials for the REO equilibrium, not the REOD
    % equilibrium which is worked on in the methods
    B1 = DocPolynom([gamma1 beta1 alpha1]);
    B2 = DocPolynom([delta1+delta2 gamma2+gamma3 Kro+beta2]);

    A1 = DocPolynom([2*gamma1 beta1]);
    A2 = DocPolynom([2*(delta1+delta2) gamma2+gamma3]);
    
    % catch if there is no operator
    if Otot == 0;
        A2 = 0*A2;
        B2 = 0*B2;       
    end
    
    phi1 = B1*B2;
    phi2 = B1 + B2*(Otot-Rtot);

    E = DocPolynom([1 0]);

    
    psi1 = E*A1*B2;
    psi2 = E*(B2+A1+A2*Otot) - Etot*B2;

    % make symbolic function
    if Etot == 0
        % without effector an analytical solution is obtainable without the
        % pattern search
        [m E R O Et Rt Ot] = evaluate_le(A1,A2,B1,B2,...
            phi1,phi2,psi1,psi2,Etot,Rtot,Otot,0);
    else
        % search for E that fits equations
        p = @(x) evaluate_le(A1,A2,B1,B2,...
            phi1,phi2,psi1,psi2,Etot,Rtot,Otot,x);
        x = patternsearch(p,Etot/2,[],[],[],[],0,Etot);
        [m E R O Et Rt Ot] = p(x);
    end
    
    % calculate all remaining species concentrations
    le.E = E;
    le.R = R;
    le.O = O;
    le.Et = Et;
    le.Rt = Rt;
    le.Ot = Ot;
    le.RE = R*E*Kre;
    le.RstE = R*E*Krrst*Krste;
    le.REE = R*E*E*(Kre^2);
    le.RstEE = R*E*E*Krrst*(Krste^2);
    le.RO = R*O*Kro;
    le.ROE = R*E*O*Kre*Kro;
    le.ROEE = R*E*E*O*Kro*(Kre^2);
    le.RstO = R*O*Krrst*Krsto;
    le.RstOE = R*O*E*Krrst*Krsto*Krste;
    le.RstOEE = R*O*E*E*Krrst*Krsto*(Krste^2);
    le.params = params;
    
end

function [m x R O Et Rt Ot] = evaluate_le...
    (A1,A2,B1,B2,phi1,phi2,psi1,psi2,Etot,Rtot,Otot,x)
    
    A = phi1(x);
    B = phi2(x);
    C = -Rtot;
    
    a = psi1(x);
    b = psi2(x);
    c = -Etot+x;
    p1 = DocPolynom([A B C]);
    p2 = DocPolynom([a b c]);
    
    r1 = roots(p1);
    r2 = roots(p2);
    r1pos = r1(r1>0);
    r2pos = r2(r2>0);
    R = mean([r1pos r2pos]);
    
    if Otot == 0
        O = 0;
    else
        O = Otot/(1+R*B2(x));
    end
    if Etot == 0
        Et = 0;
    else    
        Et = x + R*x*A1(x) + R*x*O*A2(x);
    end
        
    Rt = R*B1(x) + R*O*B2(x);
    Ot = O + O*R*B2(x);
    
    % We want to minimize the difference between the known Etot and
    % calculated Et
    m = -1/abs(Et-Etot);

end

Contact us