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_REOD(params)
function [le m] = linked_equilibrium_REOD(params)
    %        Kre         Kre
    %   RO  -----  ROE  ------ ROEE
    %    |          |            |
    %    | Kro      |            |
    %    |   Kre    |     Kre    |
    %    R  -----  RE   ------- REE
    %    |          |            |
    %    | Krr*     |            |
    %    |   Kr*e   |    Kr*e    |
    %   R*  -----  R*E  ------ R*EE
    %    |          |            |
    %    | Kr*d     |            |
    %    |   Kr*e   |     Kr*e   |
    %   R*D ----- R*ED ------- R*EED
    %
    %
    %   input: param.Kro, param.Krstd
    %          param.Kre, param.Krste
    %          param.Krrst
    %          param.Otot, param.Dtot
    %          param.Rtot, param.Etot
    %
    %   output: le (concentrations of all species structure) 
    %           minvalue (optional)
    %
    %   [le (m)] = linked_equilibrium_REOD(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;
    Krstd = params.Krstd;
    Rtot = params.Rtot;
    Otot = params.Otot;
    Etot = params.Etot;
    Dtot = params.Dtot;
    if Rtot == 0 
        le.E = Etot; 
        le.R = 0;
        le.O = Otot;
        le.D = Dtot;
        le.Et = Etot;
        le.Rt = Rtot;
        le.Ot = Otot;
        le.Dt = Dtot;
        le.RE = 0;
        le.RstE = 0;
        le.REE = 0;
        le.RstEE = 0;
        le.RO = 0;
        le.ROE = 0;
        le.ROEE = 0;
        le.RstD = 0;
        le.RstDE = 0;
        le.RstDEE = 0;
        le.params = params;
        m=-inf;
        return
    end
    % Eqns. 21-28
    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*Krstd;
    gamma3 = 2*Krrst*Krstd*Krste;
    delta2 = Krrst*Krstd*(Krste^2);
    
    % Make polynomials. Eqns. 30-32, 36-38
    B1 = DocPolynom([gamma1 beta1 alpha1]);
    B2 = DocPolynom([delta1 gamma2 Kro]);
    B3 = DocPolynom([delta2 gamma3 beta2]);

    A1 = DocPolynom([2*gamma1 beta1]);
    A2 = DocPolynom([2*delta1 gamma2]);
    A3 = DocPolynom([2*delta2 gamma3]);
    % zero out coefficients that would disappear if there is no O or D
    if Dtot == 0
        A3 = 0*A3;
        B3 = 0*B3;        
    end
    if Otot == 0;
        A2 = 0*A2;
        B2 = 0*B2;       
    end
    
    % Make higher order polynomials.  Eqns. 42-44, 45-47
    phi1 = B1*B2*B3;
    phi2 = B1*(B2+B3) + B2*B3*(Otot+Dtot-Rtot);
    phi3 = B1 + B2*(Otot-Rtot) + B3*(Dtot - Rtot);

    E = DocPolynom([1 0]);

    psi1 = E*A1*B2*B3;
    psi2 = E*(B2*B3 + A1*(B2+B3) + A2*B3*Otot + A3*B2*Dtot) - B2*B3*Etot;
    psi3 = E*(B2+B3+A1+A2*Otot + A3*Dtot) - Etot*(B2+B3);

    % make symbolic function
    if Etot == 0
        % E = 0
        % patternsearch is not required if there is no effector
        [m E R O D Et Rt Ot Dt] = evaluate_le(A1,A2,A3,B1,B2,B3,...
            phi1,phi2,phi3,psi1,psi2,psi3,Etot,Rtot,Otot,Dtot,0);
    else
        % search for E that fits equations
        p = @(x) evaluate_le(A1,A2,A3,B1,B2,B3,...
            phi1,phi2,phi3,psi1,psi2,psi3,Etot,Rtot,Otot,Dtot,x);
        x = patternsearch(p,Etot/2,[],[],[],[],0,Etot);
        [m E R O D Et Rt Ot Dt] = p(x);
    end
    le.E = E;
    le.R = R;
    le.O = O;
    le.D = D;
    le.Et = Et;
    le.Rt = Rt;
    le.Ot = Ot;
    le.Dt = Dt;
    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.RstD = R*D*Krrst*Krstd;
    le.RstDE = R*D*E*Krrst*Krstd*Krste;
    le.RstDEE = R*D*E*E*Krrst*Krstd*(Krste^2);
    le.params = params;
    
end

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

end

Contact us