1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#' Gillespie algorithm for mRNA generating processes
#'
#' Gillespie algorithms allow synthetic data simulation via three different
#' underlying mRNA generating processes: the basic process consists of a
#' simple death-birth model of mRNA transcription and degradation; the
#' switching process considers additionally gene activation and deactivation,
#' with mRNA transcription only happening in active gene states; the
#' bursting process, transcribes mRNA in bursts with geometrically distributed burst sizes.
#' The basic_burst model combines both the basic and the burst model.
#' The IGbasic burst model describes the basic model with non-constant transcription rates, but
#' transcription rates follow an inverse Gaussian distribution governed by one parameter, the mean parameter
#' of the inverse Gaussian distribution. Additionally a burst transcription occures (with NB distributed
#' burst sizes), the whole burst (rate and burst sizes) are determined by the rate parameter.
#'

#' @param n Number of observations
#' @param r.degr mRNA degradation rate (all models)
#' @param r.act DNA activation rate (Switching Model)
#' @param r.deact DNA deactivation rate (Switching Model)
#' @param r.on Transcription rate during gene activation (Switching model)
#' @param r.burst Bursty transcription rate (Bursting model, Basic Burst model and IG Basic Burst model)
#' @param s.burst Mean burst size (Bursting Model and  Basic Burst model)
#' @param r.mu Mean parameter for the inverse Gaussian distribution (IG Basic Burst model)
#' @name gmRNA
#' @rdname gmRNA
#' @export
#' @examples
#' x <- gmRNA_basic(100, 0.75, 0.001)
#' plot(density(x))
gmRNA_basic <- function(n, r.on, r.degr) {
  cpp_gmRNA_basic(n, r.on, r.degr)
}


#' @rdname gmRNA
#' @export
#' @examples
#' x <- gmRNA_switch(100, 0.23, 0.15, 0.75, 0.001)
#' plot(density(x))
gmRNA_switch <- function(n, r.act, r.deact, r.on, r.degr) {
  cpp_gmRNA_switch(n, r.act, r.deact, r.on, r.degr)
}


#' @rdname gmRNA
#' @export
#' @examples
#' x <- gmRNA_burst(10, 0.15, 0.75, 0.001)
#' plot(density(x))
gmRNA_burst <- function(n, r.burst, s.burst, r.degr) {
  cpp_gmRNA_burst(n, r.burst, s.burst, r.degr)
}


#' @rdname gmRNA
#' @export
#' @examples
#' x <- gmRNA_basic_burst(10, 0.75, 0.15, 0.5, 0.001)
#' plot(density(x))
gmRNA_basic_burst <- function(n, r.on, r.burst, s.burst, r.degr) {
    cpp_gmRNA_basic_burst(n, r.on, r.burst, s.burst, r.degr)
}

#' #' @rdname gmRNA
#' #' @export
#' #' @examples
#' #' x <- gmRNA_IGbasic_burst(10, 2, 0.5, 0.001)
#' gmRNA_IGbasic_burst <- function(n, r.mu, r.burst, r.degr) {
#'     cpp_gmRNA_IGbasic_burst(n, r.mu, r.burst, r.degr)
#' }


#' @rdname gmRNA
#' @export
#' @importFrom stats rexp
#' @examples
#' x <- gmRNA_IGbasic_burst(10, 2, 0.5, 0.1)
#' plot(density(x))
gmRNA_IGbasic_burst <- function( n, r.mu, r.burst, r.degr) {
  res <- c()
  t0 = 0
  x0 = 0
  tmax = 20/r.degr
  i = 1

  while( i <= n) {
    x = x0
    tx = t0

    lambda_draw = rInvGaus(1, r.mu, r.mu * r.burst)
    lambda1 = lambda_draw[1]
    lambda2 = r.burst
    lambda3 = r.degr * x
    lambdax = lambda1 + lambda2 + lambda3

    tau_vec = rexp(1, lambdax)
    tau = tau_vec[1]
    tau_stern = min(tau, tmax - tx)
    tx = tx+ tau_stern
    while(tx < tmax) {
      u_vec = runif(1)
      u = u_vec[1]
      if(u <= lambda1/lambdax){
        k = 1
      } else if(u <= (lambda1 + lambda2)/lambdax){
        k = 2} else k = 3


      if(tau <= tau_stern) {
        if(k == 1){
          x = x + 1
        } else if(k == 2) {
          r_vec = rnbinom(1, 1/2, r.burst/(r.burst + 2 * r.mu))
          x = x + r_vec[1]
        } else {
          x= x - 1
        }

      }


      lambda_draw = rInvGaus(1, r.mu, r.mu * r.burst)
      lambda1 = lambda_draw[1]
      lambda2 = r.burst
      lambda3 = r.degr * x
      lambdax = lambda1 + lambda2 + lambda3


      tau_vec = rexp(1, lambdax)
      tau = tau_vec[1]
      tau_stern = min(tau, tmax - tx)
      tx =  tx + tau_stern
    }
    res[i] = x
    i <- i + 1
  }
  res
}