get.sim.tl {RestoreNet}

\tau-leaping simulation algorithm

Description

Simulate a trajectory of length S for a stochastic reaction network.

Usage

get.sim.tl(Yt, theta, S, s2 = 0, tau = 1, rct.lst, verbose = TRUE)

Arguments

Yt

starting point of the trajectory

theta

vector parameter for the reactions.

S

length of the simulated trajectory.

s2

noise variance (defaults to 0).

tau

time interval length (defaults to 1).

rct.lst

list of biochemical reactions.

verbose

(defaults to TRUE) Logical value. If TRUE, then information messages on the simulation progress are printed to the console.

Details

This function allows to simulate a trajectory of a single clone given an initial conditions Y_0 for the cell counts, and obeying to a particular cell differentiation network defined by a net-effect (stoichiometric) matrix V and an hazard function h(). The function allows to consider only three cellular events, such as cell duplication (Y_{it} \rightarrow 1), cell death (Y_{it} \rightarrow \emptyset) and cell differentiation (Y_{it} \rightarrow Y_{jt}) for a clone-specific time counting process

Y_t = (Y_{1t},\dots,Y_{Nt})

observed in $N$ distinct cell lineages. In particular, the cellular events of duplication, death and differentiation are respectively coded with the character labels \texttt{"A->1"}, \texttt{"A->0"}, and \texttt{"A->B"}, where \texttt{A} and \texttt{B} are two distinct cell types. The output is a 3-dimensional array Y whose ijk-entry Y_{ijk} is the number of cells of clone k for cell type j collected at time i. More mathematical details can be found in the vignette of this package.

Value

A S \times p dimensional matrix of the simulated trajectory.

Examples

rcts <- c("A->1", "B->1", "C->1", "D->1",
          "A->0", "B->0", "C->0", "D->0",
          "A->B", "A->C", "C->D") ## set of reactions
ctps <- head(LETTERS,4)
nC <- 3 ## number of clones
S <- 10 ## trajectory length
tau <- 1 ## for tau-leaping algorithm
u_1 <- c(.2, .15, .17, .09*5,
         .001, .007, .004, .002,
         .13, .15, .08)
u_2 <- c(.2, .15, .17, .09,
         .001, .007, .004, .002,
         .13, .15, .08)
u_3 <- c(.2, .15, .17*3, .09,
         .001, .007, .004, .002,
         .13, .15, .08)
theta_allcls <- cbind(u_1, u_2, u_3) ## clone-specific parameters
rownames(theta_allcls) <- rcts
s20 <- 1 ## additional noise
Y <- array(data = NA,
           dim = c(S + 1, length(ctps), nC),
           dimnames = list(seq(from = 0, to = S*tau, by = tau),
                           ctps,
                           1:nC)) ## empty array to store simulations
Y0 <- c(100,0,0,0) ## initial state
names(Y0) <- ctps
for (cl in 1:nC) { ## loop over clones
  Y[,,cl] <- get.sim.tl(Yt = Y0,
                        theta = theta_allcls[,cl],
                        S = S,
                        s2 = s20,
                        tau = tau,
                        rct.lst = rcts,
                        verbose = TRUE)
}
Y