Supplementary Information for "CplexA: a Mathematica package to study macromolecular-assembly control of gene expression"

 Jose M. G. Vilar1, 2 and Leonor Saiz(( )^3)
1Biophysics Unit (CSIC-UPV/EHU) and Department of Biochemistry and Molecular Biology, University of the Basque Country, P.O. Box 644, 48080 Bilbao, Spain
2IKERBASQUE, Basque Foundation for Science, 48011 Bilbao, Spain
(( )^3)Department of Biomedical Engineering, University of California, 451 East Health Sciences Drive, Davis, CA 95616, USA

SetDirectory@NotebookDirectory[];

<<CplexA`

Detailed calculations for Figure 1

(* Experimental data for the transcriptional activity of the PRM and PR promoters of phage lambda from ref. Dodd,I.B.,Shearwin,K.E.,Perkins,A.J.,Burr,T.,Hochschild, A.and Egan,J.B.(2004) Cooperativity in longrange gene regulation by the lambda CI repressor,Genes Dev,18,344-354 *)
PRMdata:={{0.058, 66.367},{0.373, 189.358},{0.658, 214.233},{0.922, 188.897},{1.168, 173.696},{1.385, 155.271},{1.875, 133.160},{2.281, 136.384}};
PRdata:={{0.063, 1057.627},{0.373, 585.763},{0.658, 204.746},{0.921, 70.508},{1.171, 27.119},{1.384, 18.983},{1.875, 12.203},{2.278, 9.492}};

(* Plot of the experimental data to use later on *)
plotPRM:=ListPlot[PRMdata,PlotMarkers->{Automatic,17},PlotRange->{{0,2.5},{0,250}},ImageSize->300,BaseStyle->{FontFamily->"Helvetica",FontSize->12}];
plotPR:=ListPlot[PRdata,PlotMarkers->{Automatic,17},PlotRange->{{0,2.5},{0,1200}},ImageSize->300,BaseStyle->{FontFamily->"Helvetica",FontSize->12}];

(* Computes the CI dimer concentration as a function of CI monomer concentrations *)
NF[nn_]:=9.38419*10^-14+7.0922*10^-10 nn-1.15373*10^-11 Sqrt[0.0000661586+1. nn]

(* State variables *)
S:={sOL1,sOL2,sOL3,sOR1,sOR2,sOR3,sL}

(* Free energy *)
ΔGλ:=sOR1 (-12.7-RT Log[n])+sOR2 (-10.7-RT Log[n])+sOR3 (-10.2-RT Log[n])-3 sOR1 sOR2-3 sOR2 sOR3+3 sOR1 sOR2 sOR3+sOL1 (-13.8-RT Log[n])+sOL2 (-12.1-RT Log[n])+sOL3 (-12.4-RT Log[n])-2.5 sOL1 sOL2-2.5 sOL2 sOL3+2.5 sOL1 sOL2 sOL3+sL (21-21.2 sOL1 sOL2 sOR1 sOR2-3 sOL3 sOR3);

(* Transcriptional activity of the PRM promoter *)
Γrm:=(1-sOR3)(45+((240-45)sL+(460-45)(1-sL))sOR2);

(* Average transcriptional activity of the PRM promoter *)
AΓrm=AvConf[Γrm,ΔGλ,S];

(* Computed average transcriptional activity of the PRM promoter (red line) and its experimental counterpart (blue circles) as a function of the normalized CI dimer concentration *)
Show[plotPRM,Plot[AΓrm/.{RT->0.6,n->NF[nn]},{nn,0,2.5},PlotStyle->{Red,Thick}]]

(* Transcriptional activity of the PR promoter *)
Γr:=(1-sOR1)1200;

(* Average transcriptional activity of the PR promoter *)
AΓr=AvConf[Γr,ΔGλ,S];

(* Computed average transcriptional activity of the PR promoter (red line) and its experimental counterpart (blue circles) as a function of the normalized CI dimer concentration *)
Show[plotPR,Plot[AΓr/.{RT->0.6,n->NF[nn]},{nn,0,2.5},PlotStyle->{Red,Thick},PlotRange->All]]

Sum over an exponentially large number of states in sub-exponential times

(* Definition of a linear array of binding sites with next-neighbor interactions for molecules with concentration n. The function NNChain returns a list with the state variables, the free energy of the chain, the average occupancy, and the time needed for the computation. NumSites is the number of sites; NNint is the energy of the next-neighbor interaction between bound molecules; the free energy of binding is -1; and the thermal energy is 1. *)
NNChain[NumSites_,NNint_]:=Module[{S,ΔGB,res},
S=Table[ToExpression["s"<>ToString[i]],{i,1,NumSites}];
ΔGB=Simplify[Sum[(-1-Log[n])S[[i]],{i,1,NumSites}]+ NNint Sum[S[[i]]S[[i+1]],{i,1,NumSites-1}]];
res=Timing[AvConf[Mean[S],Simplify[ΔGB],S,1]];{S,ΔGB,res[[2]],res[[1]]}]

(* Analytical expression of the average occupancy of a linear array of 10 binding sites *)
NNChain10=NNChain[10,NNint][[3]]

(* Plot of the analytical expression of the average occupancy of a linear array of 10 binding sites as a function of the energy of interaction and the logarithm of the concentration *)
Plot3D[NNChain10/.n->10^x,{NNint,-5,5},{x,-2,2},AxesLabel->{"NNint","Log[n]","Occupancy","",""}]

TimesTab=Table[{NumSites,NNChain[NumSites,-1][[4]]},{NumSites,1,40}]

(* The CPU time (in seconds) needed to compute the average occupancy is plotted as a function of the number of binding sites of the linear array *)
ListLogPlot[TimesTab,Joined->True,PlotRange->All]

TimeRatio40to20=NNChain[40,-1][[4]]/NNChain[20,-1][[4]]

NumberStatesRatio40to20=2^40/2^20