Code covered by the BSD License

# Lac operon MWC linked equilibrium solver

### Matthew Sochor (view profile)

01 Mar 2013 (Updated )

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

```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)
%
%
%   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```