modelK2: Gene regulation with constant activation and translation

In CharlotteJana/pdmppoly: Moment Approximation for Polynomial PDMPs

Description

This PDMP models the most simple situation of gene regulation, where we have one gene and a constant activation rate without a further regulation mechanism. Transcription and translation are modeled separately which leads to a model wit two continous variables (the first (f1) representing the mRNA and the second (f2) representing the protein arising from translation). In PROM, this model is referred to as Model K2, therefore it is named genePdmpK2 and genePolyK2 here.

Usage

genePdmpK2

genePolyK2

Format

genePdmpK2 is an object of class pdmpModel,
genePolyK2 is an object of class polyPdmpModel.

Simulation

The simulations in PROM were done with slot times set to

• from = 0, to = 1000, by = 0.1.

The following parameter sets were simulated:

• k01 = 0.01, k10 = 0.01, a1 = 1, b1 = 0.06, a2 = 0.5, b2 = 0.02

• k01 = 0.01, k10 = 0.01, a1 = 1, b1 = 0.025, a2 = 0.5, b2 = 0.02

• k01 = 0.01, k10 = 0.03, a1 = 1, b1 = 0.025, a2 = 0.5, b2 = 0.0025

Source

The model, including most of the parameter sets, are described in [RajCo2006] and [Zeiser2009]. The parameter values do not rely on real data.

Examples

library(spray)
#------ code to generate the pdmpModel version -----

genePdmpK2 <- new("pdmpModel",
   descr = "Model K2: constant activation with translation",
   parms = list(b1 = 0.025, a1 = 1, k10 = 0.01, k01 = 0.01, a2 = 0.5, b2 = 0.02),
   init = c(f1 = 0.5, f2 = 0, d = 1), 
   discStates = list(d = 0:1),
   dynfunc = function(t, x, parms) {
     df <- with(as.list(c(x, parms)), c(a1*d - b1*f1, a2*f1-b2*f2))
     return(c(df, 0))
   }, 
   ratefunc = function(t, x, parms) {
     return(with(as.list(c(x, parms)), switch(d + 1, k01, k10)))
   }, 
   jumpfunc = function(t, x, parms, jtype) {
     c(x[1:2], 1 - x[3])
   }, 
   times = c(from = 0, to = 100, by = 0.1), 
   solver = "lsodar")

#------ code to generate the polyPdmpModel version -----

genePolyK2 <- new("polyPdmpModel",
   descr = "Model K2: constant activation with translation (polynomial version)",
   parms = list(b1 = 0.025, a1 = 1, k10 = 0.01, k01 = 0.01, a2 = 0.5, b2 = 0.02),
   init = c(f1 = 0.5, f2 = 0, d = 1), 
   discStates = list(d = 0:1),
   dynpolys = quote(list(
     list(overall = linear(c(-b1, 0, a1))), #df1/dt = -b1f1 + a1d
     list(overall = linear(c(a2, -b2, 0)))  #df2/dt = a2f1 - b2f2
   )),
   ratepolys = quote(list(
     list(k01,k10)
   )),
   jumpfunc = function(t, x, parms, jtype) {
     c(x[1:2], 1 - x[3])
   }, 
   times = c(from = 0, to = 100, by = 0.1),
   solver = "lsodar")

#------- comparison of the models --------------

identical(sim(genePdmpK2, outSlot = FALSE, seed = 20),
          sim(genePolyK2, outSlot = FALSE, seed = 20))

CharlotteJana/pdmppoly documentation built on Sept. 4, 2019, 4:40 p.m.