This document was made in the RStudio environment, and expects that you have installed the deSolve, rootSolve, coda and FME R-packages beforehand (which can be done by Install Packages in the Tools menu of RStudio). To start a session in Grind first source the grind.R script:

source("grind.R")
## grind.R (03-05-2024) was sourced

Grind needs to sourced to define its functions. This document uses run(), plane(), newton() and continue().

A Lac-operon model

Bacteria use several external substrates for cellular growth and they can switch the corresponding intracellular pathways on and off depending on the availability of a resource in their environment. One of the possible substrates is the sugar lactose, and the regulation of the Lac-operon was one of the first circuits of gene regulation that was revealed (Jacob & Monod, 1961). Regulation of the Lac-operon involves a positive feedback because the sugar allolactose, \(A\), that is formed from lactose intra-cellularly activates gene expression by deactivating a repressive transcription factor. The import of lactose depends on an enzyme permease that is encoded within the Lac-operon. We use \(M\) for the concentration of mRNA, and assume that the concentration of permease is proportional to \(M\). Coining \(L\) for the external lactose concentration, and assuming a steep dependence of the transcription rate \(c\) on \(A\) we propose \(A'=ML-\delta A\) and \(M'=c_0+\frac{cA^n}{1+A^n} -d M\), which is written as

model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dA = M*L - delta*A
    dM = c0 + c*A^n/(1+A^n) - d*M
  return(list(c(dA, dM)))  
  }) 
}
p <- c(L=1,c=1,c0=0.05,d=1,delta=0.5,n=4)
s <- c(A=0,M=0)

with parameters chosen such that the model has alternative steady states, The \(c_0\) parameter defines a basal transcription rate. Next make a phase plane and determine the stability of the three steady states:

plane(xmax=3,ymax=1.2)
low <- newton(s,plot=TRUE)
##         A         M 
## 0.1002016 0.0501008 
## Stable node, eigenvalues:
## [1] -1.0079214 -0.4920786
mid <- newton(c(A=1,M=0.4),plot=TRUE)
##         A         M 
## 0.9068065 0.4534033 
## Saddle point (unstable), eigenvalues:
## [1] -1.8102413  0.3102413
hig <- newton(c(A=2,M=1),plot=TRUE)

##        A        M 
## 1.977162 0.988581 
## Stable node, eigenvalues:
## [1] -1.1732318 -0.3267682

Make a bifurcation diagram by following the state in the middle while varying the external lactose concentration

continue(mid,x="L",y="A",xmax=4,ymax=4)
## Starting at L = 1 with:
##         A         M 
## 0.9068065 0.4534033

## Turning point point at L = 2.70625 
## Turning point point at L = 0.8221875
## NULL

which reveals that the operon is bistable when \(0.82<L<2.71\), and that the Lac-operon is only “on” when there is sufficient lactose in the environment. Heavy red lines depict stable states, and light blue line unstable steady states.

Saddle

The unstable steady state in the middle is a saddle point. One can depict the stable and unstable direction in the saddle by plotting the two eigenvectors. Above we saw that the first and the second eigenvalues defines the stable and unstable direction (manifold) respectively. Calling newton() with the option silent=TRUE returns a list containing the state, the eigenvalues and the eigenvectors. Hence

mid <- newton(c(A=1,M=0.5),silent=TRUE)
vec1 <- mid$vectors[,1]
vec2 <- mid$vectors[,2]
x <- mid$state[1]; y <- mid$state[2]
plane(xmax=1.5,ymax=0.75,vector=TRUE,legend=FALSE)
arrows(x+0.1*vec1[1],y+0.1*vec1[2],x,y,col="darkgreen",length=0.1)
arrows(x-0.1*vec1[1],y-0.1*vec1[2],x,y,col="darkgreen",length=0.1)
arrows(x,y,x+0.2*vec2[1],y+0.2*vec2[2],col="darkorange",length=0.1)
arrows(x,y,x-0.2*vec2[1],y-0.2*vec2[2],col="darkorange",length=0.1)

plots the stable direction in green and the unstable direction in orange.

Separatrix

The saddle point spans up the separatrix defining the basins of attraction of the two stable points. One can draw a separatrix by integrating backwards along the eigenvector pointing in the stable direction. Thus, define a model returning the reverse direction \(-A'\) and \(-M'\):

back <- function(t, state, parms) {
  if (min(state)<0) return(list(c(0,0)))  # Safety valve
  with(as.list(c(state,parms)), {
    dA = M*L - delta*A
    dM = c0 + c*A^n/(1+A^n) - d*M
  return(list(c(-dA, -dM)))       # Negative derivatives
  }) 
}

Plot a phase portrait, and step away from saddle point along the eigen vector in both directions and run

plane(tmax=10,tstep=0.1,xmax=3,ymax=1.1,portrait=TRUE,legend=FALSE)
s <- mid$state + 0.01*vec1
run(10,tstep=0.01,odes=back,traject=TRUE,lwd=3,col=3,pch=".")
##             A             M 
## -2.777740e-06  1.506721e+00
s <- mid$state - 0.01*vec1
run(10,tstep=0.01,odes=back,traject=TRUE,lwd=3,col=3,pch=".")

##             A             M 
##  1.263326e+00 -2.570823e-06

where the pch="." replaces the big bullet at the start of the trajectory.

References

Jacob, F. & Monod, J. (1961). Genetic regulatory mechanisms in the synthesis of proteins. J. Mol. Biol. 3, 318-356.

More information

Check Grind’s homepage for more tutorials