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