Jose M. G. Vilar
SUMMARY:
CplexA is a software package (available for Mathematica, Python, and Matlab) to study the effects of macromolecular assembly on cellular processes. A main challenge for conventional computational methods to study macromolecular assembly is tackling the exponential increase of the number of configurational states with the number of components. CplexA uses functional programming to efficiently compute probabilities and average properties over such exponentially large number of states from the energetics of the interactions. The package is particularly suited to study gene expression at complex promoters controlled by multiple, local and distal, DNA binding sites for transcription factors. CplexA is freely available at http://sourceforge.net/projects/cplexa/.
REFERENCES:
1. The main reference for CplexA is Vilar, J.M.G. and Saiz, L. (2010) CplexA: a Mathematica package to study macromolecular-assembly control of gene expression, Bioinformatics 26, 2060-2061.
2. The biophysical background of the methods used by CplexA are described in Saiz, L. and Vilar, J.M.G. (2006) Stochastic dynamics of macromolecular-assembly networks, Mol Syst Biol, 2, 2006.0024.
To run the python tutorial, the file CplexA.py must be in the Path, the default list of directories searched by Python to find an external file. A possibility is to place CplexA.py in the same directory as the tutorial file. CplexA uses SymPy, a Python library for symbolic mathematics.
from pylab import *
from numpy import *
from sympy import *
from CplexA import *
The main function of CplexA is AvConf(F,
Note that different units can be used for the thermal energy, for instance kT (AvConf(F,
A transcription factor with concentration n1 binds a DNA site described by the state variable s1. The variable indicates whether (s1=1) or not (s1=0) the transcription factor is bound. S is the set of state variables describing the system (in this case, it consists only of s1); DG is the free energy of the system; DG0 is the standard free energy of binding at 1M; and F is the transcription rate, which is Fmax when the transcription factor is bound (s1=1) or zero when it is not (s1=0).
The example below computes the average transcription rate AF.
S='s1'
DG='(DG0-RT*ln(n1))*s1'
F='Fmax*s1'
AF=AveConf(F,DG,S)
print simplify(AF)
Same example as before, but now the free energy is expressed in terms of the dissociation constant Kd (recall the equivalency
DG='-RT*ln(n1/Kd)*s1'
AF=AveConf(F,DG,S)
print simplify(AF)
Two molecules with concentrations n1 and n2 compete for binding to overlapping sites (binding of one excludes the binding of the other). The variables s1 and s2 indicate whether molecules 1 and 2 are bound, respectively. The expression below computes the probability that molecule 1 is bound, which is given by the average value of s1.
Note the term inf s1 s2 in the free energy. It implements the logical condition that two molecules cannot be bound simultaneously by assigning an infinite free energy to the corresponding state. To obtain the correct thermodynamic result, the symbol "inf" is considered to be finite until the end result, and then the limit to infinity is taken.
AF=AveConf('s1','-RT*ln(n1/K1)*s1-RT*ln(n2/K2)*s2+inf*RT*s1*s2','s1 s2')
print simplify(AF)
Similar to the previous example, but now the binding of one does not exclude the other. The term e12 is the interaction energy between the molecules bound. The following expression computes the probability of having the two sites occupied simultaneously:
AF=AveConf('s1*s2','-RT*ln(n1/K1)*s1-RT*ln(n2/K2)*s2+e12*s1*s2','s1 s2')
print simplify(AF)
The lac repressor is a protein with two DNA binding domains. It can bind to one site or, by looping the intervening DNA, to two distal sites simultaneously.
In addition to the domain-binding state variables to each site, s1 and s2, there is a looping state variable sL that indicates wether the DNA is looped (sL=1) or not (sL=0). Note that the each binding domain at sites 1 and 2 can be from two different lac repressors or from a single one when there is looping. The term inf (1 - s1 s2) sL means that looping can only happen with two domains bound, otherwise the free energy is infinite. The following expression computes the repression level:
AF=AveConf('(1-s1)','-RT*ln(n/K1)*s1-RT*ln(n/K2)*s2\
+(gL+RT*ln(n)*s1*s2)*sL+inf*RT*(1-s1*s2)*sL','s1 s2 sL')
print 1+simplify(1/AF-1)
The quantity gL is the contribution to the free energy due to the conformational change of DNA (free energy of looping).
The repression level can also be expressed in terms of the local concentration cL, defined as
print 1+simplify(subs(1/AF,'gL','-RT*log(cL)')-1)
The CI protein of phage
The free energy of the system DG takes into account interactions between CI dimers bound at neighboring sites and, when there is looping (sL=1), at different operators as well.
S='sOL1 sOL2 sOL3 sOR1 sOR2 sOR3 sL'
DG='-2.5*sOL1*sOL2 - 2.5*sOL2*sOL3 + 2.5*sOL1*sOL2*sOL3\
- 3*sOR1*sOR2 - 3*sOR2*sOR3 + 3*sOR1*sOR2*sOR3\
+ sL*(21 - 21.2*sOL1*sOL2*sOR1*sOR2 - 3*sOL3*sOR3)\
+ sOL1*(-13.8 - RT*log(n)) + sOR1*(-12.7 - RT*log(n))\
+ sOL3*(-12.4 - RT*log(n)) + sOL2*(-12.1 - RT*log(n))\
+ sOR2*(-10.7 - RT*log(n)) + sOR3*(-10.2 - RT*log(n))'
Transcription at the PRM promoter is given by F_rm, which leads to the average AF_rm. There is also another promoter, PR, with transcription given by F_r, which leads to the average AF_r.
F_rm='(45 + (415*(1 - sL) + 195*sL)*sOR2)*(1 - sOR3)'
AF_rm=AveConf(F_rm,DG,S)
AF_rm_n=simplify(subs(AF_rm,'RT',0.6))
print AF_rm_n
F_r='1200*(1 - sOR1)'
AF_r=AveConf(F_r,DG,S)
AF_r_n=simplify(subs(AF_r,'RT',0.6))
print AF_r_n
Note that in both cases it is possible to obtain an exact analytical expression for the effective transcription rates, AF_rm and AF_r, as a function of the CI dimer concentration. The expression NF below gives the CI dimer concentration as function of the normalized CI monomer concentration.
def NF(nn):
return 9.38419e-14 + 7.0922e-10*nn - 1.15373e-11*sqrt(0.0000661586 + 1.*nn)
nnvals=arange(0, 2.5, 0.1)
nvals=array([NF(nn) for nn in nnvals])
Plotted below is AF_rm (solid red line) as a function of the normalized CI monomer concentration. The filled circles correspond the experimentally measured activity of the PRM promoter (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 = array([[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]])
AF_rm_f=pythonFunction(AF_rm_n,'n')
plot(nnvals, AF_rm_f(nvals), 'r--', PRMdata[:,0], PRMdata[:,1],'ob');
xlabel('$[CI]$',fontsize=14); ylabel('$AF_{rm}$',fontsize=14); show()
Plotted below is AF_r (solid red line) as a function of the normalized CI monomer concentration. The filled circles correspond the experimentally measured activity of the PR promoter ( 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).