modelK: Gene regulation with constant activation

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 considered as one step and are not modeled separately. In PROM, this model is referred to as Model K, therefore it is named genePdmpK and genePolyK here.

Usage

genePdmpK

genePolyK

Format

genePdmpK is an object of class pdmpModel,
genePolyK 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, a = 1, b = 0.06

• k01 = 0.01, k10 = 0.01, a = 1, b = 0.005

• k01 = 0.01, k10 = 0.03, a = 1, b = 0.025

• k01 = 0.03, k10 = 0.01, a = 1, b = 0.025

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 -----

genePdmpK <- new("pdmpModel",
   descr = "Model K: constant activation",
   parms = list(b = 0.005, a = 1, k10 = 0.01, k01 = 0.01),
   init = c(f = 0, d = 1),
   discStates = list(d = 0:1),
   dynfunc = function(t, x, parms) {
     df <- with(as.list(c(x, parms)), a*d - b*f)
     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], 1 - x[2])
   },
   times = c(from = 0, to = 100, by = 0.1),
   solver = "lsodar")

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

genePolyK <- new("polyPdmpModel",
     descr = "Model K: constant activation (polynomial version)",
     parms = list(b = 0.005, a = 1, k10 = 0.01, k01 = 0.01),
     init = c(f = 0, d = 1),
     discStates = list(d = 0:1),
     dynpolys = quote(list(
       list(overall = linear(c(-b, a)))
     )),
     ratepolys = quote(list(
       list(k01, k10)
     )),
     jumpfunc = function(t, x, parms, jtype) {
       c(x[1], 1 - x[2])
     },
     times = c(from = 0, to = 100, by = 0.1),
     solver = "lsodar")

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

identical(sim(genePdmpK, outSlot = FALSE, seed = 20),
          sim(genePolyK, outSlot = FALSE, seed = 20))

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