Model Development
P2 genetic switch
We started modeling the P2 genetic switch by describing the system graphically in terms of C protein dynamics, based on a similar approach applied to bacteriophage λ, as described by Mestivier et al. (3)
Based on the graphical representation, we derived differential equations for each of the reactions that the C protein is involved in. Notation: x - concentration of C; y - concentration of C2 (c protein dimer).
Degradation of C protein
$$ \frac{dx}{dt}=-k_0x$$
Dimerization of C protein
$$ \frac{dx}{dt}=-k_{12}x^2+2k_{21}y;\ \frac{dy}{dt}=k_{12}x^2-k_{21}y$$
Synthesis of C protein
$$ \frac{dx}{dt}=nk_tD_0P$$
Binding of C dimer to Pc promoter at low concentrations
$$ \frac{dy}{dt}=-k_1yD_0+k_{-1}D_1;\ \frac{dD_0}{dt}=-k_1yD_0+k_{-1}D_1;\ \frac{dD_1}{dt}=k_1yD_0-k_{-1}D_1$$
Stimulated synthesis of C protein
$$ \frac{dx}{dt}=nk_{t2}D_1P$$
Binding of C dimer to Pc promoter at high concentrations
$$ \frac{dy}{dt}=-k_2yD_1+k_{-2}D_2;\ \frac{dD_1}{dt}=-k_2yD_1+k_{-2}D_2;\ \frac{dD_2}{dt}=k_2yD_1-k_{-2}D_2\ $$
We then gathered all differential equations into one model, describing the dynamics of C protein:
$$\frac{dx}{dt}=-k_0-k_{12}x^2+2k_{21}y+nk_tD_0P+nk_{t2}D_1P$$
Following that, we made an assumption that dimerization and binding to the operator sites are faster chemical reactions than synthesis, and for these reactions:
$$\frac{dx}{dt}=0$$
This assumption allowed us to neglect some parts of the main equation, as well as derive expressions of some constants and species involved in the reactions. Specifically, for reaction #2 (dimerization):
$$\frac{dx}{dt}=-k_{12}x^2+2k_{21}y=0;\ \frac{dy}{dt}=\ k_{12}x^2-k_{21}y=0\rightarrow k_{12}x^2=k_{21}y\rightarrow y=\ \frac{k_{12}}{k_{21}}x^2=K_Dx^2$$
For reaction #4 (binding):
$$\frac{dy}{dt}=-k_1yD_0+k_{-1}D_1\rightarrow D_1=\frac{k_1}{k_{-1}}yD_0\rightarrow\ D_1=K_1D_0y\rightarrow D_1=K_1D_0\left(K_Dx^2\right)=K_1K_DD_0x^2$$
For reaction #6 (binding):
$$\frac{dy}{dt}=-k_2yD_1+k_{-2}D_2=0\rightarrow\ k_2yD_1=k_{-2}D_2\rightarrow D_2=\frac{k_2}{k_{-2}}yD_1=K_2yD_1=K_D^2K_1K_2x^4D_0$$
Another assumption we made is that when C dimer binds to D0, the affinity constant of C dimer to D0 changes by a factor σ , that is
$${\ K}_2=σK_1$$
This results in expression:
$$D_2=K_D^2σK_1^2x^4D_0$$
The general model of C protein concentration now looks like this:
$$\frac{dx}{dt}=-k_0x+nk_tD_0P+nk_{t2}D_1P$$
We make yet another assumption, that total phage DNA concentration is constant, and derive an expression of D0.
$$D=D_0+D_1+D_2=D_0+K_1K_DD_0x^2+σK_1^2K_D^2x^4D_0$$
Inserting this expression into the general model equation results in an updated model:
$$\frac{dx}{dt}=-k_0x+nk_tD_0P+nk_{t2}K_1K_DD_0x^2P=nD_0P\left(k_t+k_{t2}K_1K_Dx^2\right)-k_0x=\frac{nDP\left(k_t+k_{t2}K_1K_2K_Dx^2\right)}{1+K_1K_Dx^2+σK_1^2K_D^2x^4}-k_0x$$
First we relate the two transcription parameters:
$$k_{t2}=τk_t$$
Then we introduce a new parameter:
$$w=nk_tDP$$
The updated general model now looks like this:
$$\frac{1}{w}\frac{dx}{dt}=\frac{1+τK_1K_Dx^2}{1+K_1K_Dx^2+σK_1^2K_D^2x^4}-\frac{k_0}{w}x$$
Here we introduce a new time unit:
$$t\prime=wt$$
As well as two new parameters:
$$\gamma=\frac{k_0}{w}; u=K_DK_1$$
The final model is:
$$\frac{dx}{dt\prime}=\frac{1+\tau u x^2}{1+ux^2+\sigma u^2x^4}-\gamma x=f(x)$$
We use this equation to simulate C protein dynamics over time using python. To compute the concentration of species x (C protein) versus time t, we use ‘Euler scheme of numerical integration’. We set the initial conditions of species x and the parameter dt to compute small variation in concentration of C protein (dx) over a small time interval (dt) using f (x) ; the variation being f (x).dt.
$$dx/dt = f(x)$$
Hence, the concentration of protein C at a time t + dt is the concentration x at a time t + f(x). dt x(t + dt) = x(t) + dt. f (x) , where x(t) = concentration of x at time t. Having x(t + dt), we can compute the small variation of species x during the next time interval dt, and compute x ((t + dt) + dt) i.e. x(t + 2.dt) and thus follow an iterative procedure further to plot time evolution of x in this case. We also find steady states of our system by plotting f (x) over x. The system achieves a steady state, when C protein would not have a net new concentration i.e. change in x in small time interval dt would be zero.
Hence, the concentration of protein C at a time t + dt is the concentration x at a time t + f(x). dt
$$x(t + dt) = x(t) + dt. f (x)$$
where x(t) = concentration of x at time t. Having x(t + dt), we can compute the small variation of species x during the next time interval dt, and compute x ((t + dt) + dt) i.e. x(t + 2.dt) and thus follow an iterative procedure further to plot time evolution of x in this case.
We also find steady states of our system by plotting f (x) over x. The system achieves a steady state, when C protein would not have a net new concentration i.e. change in x in small time interval dt would be zero.
$$dx/dt =0, f(x)=0$$
Any value of x that nullifies f(x) would be a steady state concentration. The concentration of x at which f(x) crosses the horizontal axis in the plot thus gives us the steady state concentration.
In order to adapt the above general model to a stochastic setting, we use Poisson τ-leap approach - a simple technique to introduce noise into a deterministic model. As described by Tian et al, it is assumed that a number of reactions are occuring in a relatively large time interval interval [t, t +τ). The number of these reactions corresponds to a sample value generated from a Poisson distribution with a mean of a:
$$P(a_j\left({x}\right)\tau); a_j({x})\tau$$
Once this value is generated, the system can be updated in the following way:
$${x}\left(t+\tau\right)={x}\left(t\right)+\sum_{j=1}^{M}{v_jP(a_j\left({x}\right)τ)}$$
We adopt this approach to derive our general stochastic model with a τ value of 0.01. For our general model, the updated equation looks like this:
$$x\left(t+\tau\right)={x}\left(t\right)+\mathcal{P}\left[\frac{1+\tau u x^2}{1+ux^2+\sigma u^2x^4}\tau\right]-\mathcal{P}\left[\gamma x\tau\right]$$
We implement this approach in Python, using numpy.random.poisson function.
Model Plasmid
Based on the graphical representation, we have derived the following differential equations. Notation: x - concentration of C; y- concentration of C2; m - concentration of Cox; n - concentration of Cox4.
Synthesis of C protein
$$\frac{dx}{dt}=nk_{t1}D_1P$$
Degradation of C protein
$$\frac{dx}{dt}=-k_0x$$
Binding of arabinose
$$\frac{d[arabinose]}{dt}=-k_1D_0[arabinose]+k_{-1}D_1$$
$$\frac{dD_0}{dt}=-k_1D_0[arabinose]+k{-1}D_1$$
$$\frac{dD_1}{dt}=k_1D_0[arabinose]-k{-1}D_1$$
Dimerization of C protein
$$\frac{dx}{dt}=-k_{12}x^2+2k_{21}y$$
$$\frac{dy}{dt}=k_{12}x^2-k_{21}y$$
Synthesis of Cox from Pe promoter
$$\frac{dm}{dt}=Nk_{t2}D_1P$$
Tetramerization of Cox
$$\frac{dm}{dt}=-k_{14}m^4+4k_{41}n$$
$$\frac{dn}{dt}=k_{14}m^4-k_{41}n$$
Degradation of Cox
$$\frac{dm}{dt}=-k_0m$$
Binding of C2 at low concentrations to Pe promoter
$$\frac{dy}{dt}=-k_2D_1y+k_{-2}D_3$$
$$\frac{dD_1}{dt}=-k_2D_1y+k_{-2}D_3$$
$$\frac{dD_3}{dt}=k_2D_1y-k_{-2}D_3$$
Slowed synthesis of Cox because of C2 binding
$$\frac{dm}{dt}=Nk_{t3}D_3P$$
Binding of Cox4 at low concentrations to Pe promoter
$$\frac{dn}{dt}=-k_3D_1n+k_{-3}D_2$$
$$\frac{dD_1}{dt}=-k_3D_1n+k_{-3}D_2$$
$$\frac{dD_2}{dt}=k_3D_1n-k_{-3}D_2$$
Slowed synthesis of Cox because of Cox4 binding to Pe promoter
$$\frac{dm}{dt}=Nk_{t4}D_2P$$
Binding of C2 at high concentrations to Pe
$$\frac{dy}{dt}=-k_4D_3y+k_{-4}D_5$$
$$\frac{dD_3}{dt}=-k_4D_3y+k_{-4}D_5$$
$$\frac{dD_5}{dt}=k_4D_3y-k_{-4}D_5$$
Binding of Cox4 leading to inhibited synthesis
$$\frac{dn}{dt}=-k_5D_3n+k_{-5}D_6$$
$$\frac{dD_3}{dt}=-k_5D_3n+k_{-5}D_6$$
$$\frac{dD_6}{dt}=k_5D_3n-k_{-5}D_6$$
Binding of C2 leading to inhibited synthesis
$$\frac{dy}{dt}=-k_6D_2y+k_{-6}D_6$$
$$\frac{dD_2}{dt}=-k_6D_2y+k_{-6}D_6$$
$$\frac{dD_6}{dt}=k_6D_2y-k_{-6}D_6$$
Binding of Cox4 at high concentrations to Pe
$$\frac{dn}{dt}=-k_7D_2n+k_{-7}D_4$$
$$\frac{dD_2}{dt}=-k_7D_2n+k_{-7}D_4$$
$$\frac{dD_4}{dt}=k_7D_2n-k_{-7}D_4$$
We then make an assumption that dimerization, tetramerization, and binding are very fast chemical reactions, faster than synthesis. For these reactions:
From reaction #3 we derive:
$$k_1D_0[arabinose]=k_{-1}D_1$$
$$D_1=\frac{k_1}{k_{-1}}D_0[arabinose]=K_1D_0[arabinose]$$
From reaction #4 we derive:
$$k_{12}x^2=k_{21}y$$
$$y=\frac{k_{12}}{k_{21}}x^2=K_Dx^2$$
From reaction #6 we derive:
$$k_{14}m^4=k_{41}n$$
$$n=\frac{k_{14}}{k_{41}}m^4={K_Tm}^4$$
From reaction #8 we derive:
$$k_2D_1y=k_{-2}D_3$$
$$D_3=\frac{k_2}{k_{-2}}D_1y=K_2K_1D_0[arabinose]K_Dx^2$$
From reaction #10 we derive:
$$k_3D_1n=k_{-3}D_2$$
$$D_2=\frac{k_3}{k_{-3}}D_1n=K_3K_1D_0[arabinose]K_Tm^4$$
From reaction #12 we derive:
$$k_4D_3y=k_{-4}D_5$$
$$D_5=\frac{k_4}{k_{-4}}D_3y=K_4K_2K_1D_0[arabinose]K_D^2m^4$$
From reaction #13 we derive:
$$k_5D_3n=k_{-5}D_6$$
$$D_6=\frac{k_5}{k_{-5}}D_3n=K_5K_2K_1D_0[arabinose]K_Dx^2K_Tm^4$$
From reaction #14 we derive:
$$k_6D_2y=k_{-6}D_6$$
$$D_6=\frac{k_6}{k_{-6}}D_2y=K_6K_3K_1D_0[arabinose]K_Dx^2K_Tm^4$$
From reaction #15 we derive:
$$k_7D_2n=k_{-7}D_4$$
$$D_4=\frac{k_7}{k_{-7}}D_2n=K_7K_3K_1D_0[arabinose]K_T^2m^4$$
The general models with these equations plugged in look like this:
$$\frac{dx}{dt}=Nk_{t1}D_1P-k_0x=Nk_{t1}K_1D_0[arabinose]P-k_0x$$
$$\frac{dm}{dt}=Nk_{t2}D_1P-k_{02}m+Nk_{t3}D_3P+Nk_{t4}D_2P=Nk_{t2}K_1D_0[arabinose]P-k_{02}m+Nk_{t3}K_2K_1D_0[arabinose]K_Dx^2P+Nk_{t4}K_3K_1D_0[arabinose]K_Tm^4P$$
The total plasmid concentration stays the same:
$$D=D_0+D_1+D_2+D_3+D_4+D_5+D_6=D_0+K_1D_0[arabinose]+K_3K_1D_0[arabinose]K_Tm^4+K_2K_1D_0[arabinose]K_Dx^2+K_7K_3K_1D_0[arabinose]K_T^2m^4+K_4K_2K_1D_0[arabinose]K_D^2m^4+K_6K_3K_1D_0[arabinose]K_Dx^2K_Tm^4$$
We introduce new parameters:
$$w=NDP$$
$$\alpha=\frac{k_0}{w};\ \beta=\frac{k_{02}}{w}$$
The final models are:
$$\frac{dx}{dt}=\frac{K_{t1}K_1[arabinose]}{1+[arabinose](K_1+K_3K_1m^4+K_2K_1K_Dx^2+K_7K_3K_5K_T^2m^8+K_4K_2K_1K_D^2x^4+K_5K_2K_1K_Dx^2K_Tm^4)}-αx$$
$$\frac{dm}{dt}=\frac{[arabinose](K_{t2}K_1+k_{t3}K_1K_2K_Dx^2+k_{t4}K_1K_3K_Tm^4)}{1+[arabinose](K_1+K_3K_1m^4+K_2K_1K_Dx^2+K_7K_3K_5K_T^2m^8+K_4K_2K_1K_D^2x^4+K_5K_2K_1K_Dx^2K_Tm^4)}-βm$$
Switch Plasmid
Based on the graphical representation, we have derived the following differential equations for each of the reactions in the network. Notation: x - concentration of C; m - concentration of Cox; r - concentrattion of tetR, q - concentration of tetR2.
Synthesis of C protein
$$\frac{dx}{dt}=Nk_{t1}D_0P$$
Degradation of C protein
$$\frac{dx}{dt}=-k_0x$$
Synthesis of Cox and tetR
$$\frac{dm}{dt}=Nk_{t2}D_1P$$
$$\frac{dr}{dt}=Nk_{t2}D_1P$$
Degradation of Cox protein
$$\frac{dm}{dt}=-k_{02}m$$
Dimerization of tetR
$$\frac{dr}{dt}=-k_{12}r^2+2k_{21}q$$
$$\frac{dq}{dt}=k_{12}r^2-k_{21}q$$
Binding of tetR2 to D0
$$\frac{dq}{dt}=-k_1D_0q+k_{-1}D_2$$
$$\frac{dD_0}{dt}=-k_1D_0q+k_{-1}D_2$$
$$\frac{dD_2}{dt}=k_1D_0q-k_{-1}D_2$$
Slowed synthesis of C
$$\frac{dx}{dt}=Nk_{t3}D_2P$$
Binding of tetR2 to D2
$$\frac{dq}{dt}=-k_2D_2q+k_{-2}D_3$$
$$\frac{dD_2}{dt}=-k_2D_2q+k_{-2}D_3$$
$$\frac{dD_3}{dt}=k_2D_2q-k_{-2}D_3-k_{-2}D_3$$
Binding of arabinose
$$\frac{dD_0}{dt}=-k_3D_0[arabinose]+k-3D1$$
$$\frac{dD_1}{dt}=k_3D_0[arabinose]-k-3D1$$
Degradation of tetR
$$\frac{dr}{dt}=-k_{03}r$$
Then we gathered all the equations into one per each protein involved (general models).
$$\frac{dx}{dt}=Nk_{t1}D_0P-k_0x+Nk_{t3}D_2P$$
$$\frac{dm}{dt}=Nk_{t2}D_1P+k_{02}m$$
$$\frac{dr}{dt}=Nk_{t2}D_1P-k_{12}r^2$$
$$\frac{dq}{dt}=k_{12}r^2-k_{21}q-k_1D_0q-k_{-1}D_2-k_2D_2q+k_{-2}D_3$$
We made an assumption that dimerization and binding are very fast chemical reactions, faster than synthesis. Hence, reaction #5, #6, #8 and #9 are equal to 0.
From reaction #5 we derive:
$$\frac{dr}{dt}=-k_{12}r^2+2k_{21}q=0$$
$$k_{12}r^2=k_{21}q$$
$$q=\frac{k_{12}}{k_{21}}r^2={K_Dr}^2$$
From reaction #6 we derive:
$$\frac{dq}{dt}=-k_1D_0q+k_{-1}D_2=0$$
$$k_1D_0q=k_{-1}D_2$$
$$D_2=\frac{k_1}{k_{-1}}D_0q=K_1D_0q=K_1K_DD_0r^2$$
From reaction #8 we derive:
$$\frac{dq}{dt}=-k_2D_2q+k_{-2}D_3=0$$
$$k_2D_2q=k_{-2}D_3$$
$$D_3=K_2K_1K_DD_0r^2K_Dr^2=K_D^2K_2K_1D_0r^4$$
From reaction #9 we derive:
$$\frac{dD_0}{dt}=-k_3D_0[arabinose]+k_{-3}D_1=0$$
$$k_3D_0[arabinose]=k_{-3}D_1$$
$$D_1=\frac{k_3}{k_{-3}}D_0[arabinose]=K_3D_0[arabinose]$$
Filling in these equations into the general model equations yields:
$$\frac{dx}{dt}=Nk_{t1}D_0P-k_0x+Nk_{t3}D_2P=Nk_{t1}D_0P-k_0x+Nk_{t3}PK_1K_Dr^2D_0$$
$$\frac{dm}{dt}=Nk_{t2}D_1P+k_{02}m=Nk_{t2}K_3D_0[arabinose]P-k_{02}m$$
$$\frac{dr}{dt}=Nk_{t2}D_1P-k_{12}r^2=Nk_{t2}K_3D_0[arabinose]P-k_{03}r$$
The total plasmid concentration stays the same:
$$D=D_0+D_1+D_2+D_3=D_0+K_3[arabinose]D_0+K_1K_Dr^2D_0+K_D^2K_1K_2D_0r^4$$
$$D_0=\frac{D}{D_1+K_3[arabinose]+K_1K_Dr^2+K_D^2K_1K_2r^4}$$
Updated general model now looks like this:
$$\frac{dx}{dt}=\frac{Nk_{t1}DP+Nk_{t3}PK_1K_Dr^2D}{1+K_3[arabinose]+K_1K_Dr^2+K_D^2K_1K_2r^4}-\frac{k_0}{w}x$$
We introduce new parameters:
$$w=NDP$$
$$\alpha=\frac{k_0}{w};\ \beta=\frac{k_{02}}{w};\gamma=\frac{k_{03}}{w}$$
The final models are:
$$\frac{dx}{dt}=\frac{k_{t1}+k_{t3}K_1K_Dr^2}{1+K_3[arabinose]+K_1K_Dr^2+K_D^2K_1K_2r^4}-αx$$
$$\frac{dm}{dt}=\frac{k_{t2}K_3[arabinose]}{1+K_3[arabinose]+K_1K_Dr^2+K_D^2K_1K_2r^4}-βx$$
$$\frac{dr}{dt}=\frac{k_{t2}K_3[arabinose]}{1+K_3[arabinose]+K_1K_Dr^2+K_D^2K_1K_2r^4}-γx$$