R/simulateGRC.R

In sRACIPE: Systems biology tool to simulate gene regulatory circuits

#' @export
#' @title Simulate a gene regulatory circuit
#' @import SummarizedExperiment
#' @importFrom utils read.table write.table data
#' @importFrom S4Vectors metadata
#' @description Simulate a gene regulatory circuit using its topology as the
#' only input. It will generate an ensemble of random models.
#' @param circuit data.frame or character. The file containing the circuit or
#' @param config  (optional) List. It contains simulation parameters 
#' like integration method
#' (stepper) and other lists or vectors like simParams, 
#' stochParams, hyperParams, options, thresholds etc. 
#' The list simParams contains values for parameters like the 
#' number of models (numModels), 
#' simulation time (simulationTime), step size for simulations 
#' (integrateStepSize), when to start recording the gene expressions 
#' (printStart), time interval between recordings (printInterval), number of 
#' initial conditions (nIC), output precision (outputPrecision), tolerance for
#' adaptive runge kutta method (rkTolerance), parametric variation (paramRange).
#' The list stochParams contains the parameters for stochastic simulations like
#' the number of noise levels to be simulated (nNoise), the ratio of subsequent
#' noise levels (noiseScalingFactor), maximum noise (initialNoise), whether to
#' use same noise for all genes or to scale it as per the median expression of
#' the genes (scaledNoise), ratio of shot noise to additive noise (shotNoise).
#' The list hyperParams contains the parameters like the minimum and maximum 
#' production and degration of the genes, fold change, hill coefficient etc.
#' The list options includes logical values like annealing (anneal), scaling of 
#' noise (scaledNoise), generation of new initial conditions (genIC), parameters
#' (genParams) and whether to integrate or not (integrate). The user
#' modifiable simulation options can be specified as other arguments. This 
#' list should be used if one wants to modify many settings for multiple 
#' simulations.
#' @param anneal (optional) Logical. Default FALSE. Whether to use annealing
#' for stochastic simulations. If TRUE, the gene expressions at higher noise
#' are used as initial conditions for simulations at lower noise.
#' @param numModels (optional) Integer. Default 2000. Number of random models
#' to be simulated.
#' @param thresholdModels (optional)  integer. Default 5000. The number of
#' models to be used for calculating the thresholds for genes.
#' @param paramRange (optional) numeric (0-100). Default 100. The relative
#' range of parameters (production rate, degradation rate, fold change).
#' @param prodRateMin (optional) numeric. Default 1. Minimum production rate.
#' @param prodRateMax (optional) numeric. Default 100. Maximum production rate.
#' @param degRateMin (optional) numeric. Default 0.1. Minimum degradation rate.
#' @param degRateMax (optional) numeric. Default 1. Maximum degradation rate.
#' @param foldChangeMin (optional) numeric. Default 1. Minimum fold change for
#' interactions.
#' @param foldChangeMax (optional) numeric. Default 100. Maximum fold change for
#' interactions.
#' @param hillCoeffMin (optional) integer. Default 1. Minimum hill coefficient.
#' @param hillCoeffMax (optional) integer. Default 6. Maximum hill coefficient.
#' @param integrateStepSize (optional) numeric. Default 0.02. step size for
#' integration using "EM" and "RK4" steppers.
#' @param simulationTime (optional) numeric. Total simulation time.
#' @param nIC (optional) integer. Default 1. Number of initial conditions to be
#' simulated for each model.
#' @param ... Other arguments
#' @param nNoise (optional) integer. Default 0.
#' Number of noise levels at which simulations
#' are to be done. Use nNoise = 1 if simulations are to be carried out at a
#' specific noise. If nNoise > 0, simulations will be carried out at nNoise
#' levels as well as for zero noise. "EM" stepper will be used for simulations
#' and any argument for stepper will be ignoired.
#' @param simDet (optional) logical. Default TRUE.
#' Whether to simulate at zero noise as well also when using nNoise > 0.
#' @param initialNoise (optional) numeric.
#' Default 50/sqrt(number of genes in the circuit). The initial value of noise
#' for simulations. The noise value will decrease by a factor
#' \code{noiseScalingFactor} at subsequent noise levels.
#' @param noiseScalingFactor (optional) numeric (0-1) Default 0.5.
#' The factor by which noise
#' will be decreased when nNoise > 1.
#' @param shotNoise (optional) numeric. Default 0.
#' The ratio of shot noise to additive
#' noise.
#' @param scaledNoise (optional) logical. Default FALSE. Whether to scale the
#' noise in each gene by its expected median expression across all models. If
#' TRUE the noise in each gene will be proportional to its expression levels.
#' @param outputPrecision (optional) integer. Default 12.
#' The decimal point precison of
#' the output.
#' @param knockOut (optional) List of character or vector of characters.
#' Simulation after knocking out one or more genes. To knock out all the genes
#' in the circuit, use \code{knockOut = "all"}. If it is a vector, then all
#' the genes in the vector will be knocked out simultaneously.
#' @param printStart (optional) numeric (0-\code{simulationTime}). 
#' Default \code{simulationTime}. To be used only when \code{timeSeries} is 
#' \code{TRUE}.
#' The time from which the output should be recorded. Useful for time series
#' analysis and studying the dynamics of a model for a particular initial
#' condition. 
#' @param printInterval (optional) numeric (\code{integrateStepSize}-
#' \code{simulationTime - printStart}). Default 10. The separation between
#' two recorded time points for a given trajectory. 
#' To be used only when \code{timeSeries} is 
#' \code{TRUE}.
#' @param stepper (optional) Character. Stepper to be used for integrating the
#' differential equations. The options include \code{"EM"} for Euler-Maruyama
#' O(1), \code{"RK4"}
#' for fourth order Runge-Kutta O(4) and \code{"DP"} for adaptive stepper based
#' Dormand-Prince algorithm. The default method is \code{"RK4"}
#'  for deterministic
#' simulations and the method defaults to \code{"EM"}
#' for stochastic simulations.
#' @param plots (optional) logical Default \code{FALSE}.
#' Whether to plot the simuated data.
#' @param plotToFile (optional) Default \code{FALSE}. Whether to save the plots
#' to a file.
#' @param genIC (optional) logical. Default \code{TRUE}. Whether to generate
#' the initial conditions. If \code{FALSE}, the initial conditions must be
#' supplied as a dataframe to \code{circuit$ic}.
#' @param genParams (optional) logical. Default \code{TRUE}. Whether to generate
#' the parameters. If \code{FALSE}, the parameters must be
#' supplied as a dataframe to \code{circuit$params}.
#' @param integrate (optional) logical. Default \code{TRUE}. Whether to
#' integrate the differential equations or not. If \code{FALSE}, the function
#' will only generate the parameters and initial conditions. This can be used
#' iteratively as one can fist generate the parameters and initial conditions
#' and then modify these before using these modified values for integration.
#' For example, this can be used to knockOut genes by changing the production
#' rate and initial condition to zero.
#' @param rkTolerance (optional) numeric. Default \code{0.01}. Error tolerance
#' for adaptive integration method.
#' @param timeSeries (optional) logical. Default \code{FALSE}. 
#' Whether to generate time series for a single model instead of performing
#' RACIPE simulations. 
#' @param nCores (optional) integer. Default \code{1}
#' Number of cores to be used for computation. Utilizes \code{multiprocess} from 
#' \code{doFuture} pacakge. Will not work in Rstudio.
#' @return \code{RacipeSE} object. RacipeSE class inherits 
#' \code{SummarizedExperiment} and contains the circuit, parameters, 
#' initial conditions,
#' simulated gene expressions, and simulation configuration. These can be 
#' accessed using correponding getters. 
#' @examples
#' data("demoCircuit")
#' rSet <- sRACIPE::sracipeSimulate(circuit = demoCircuit)
#' @section Related Functions:
#'
#' \code{\link{sracipeSimulate}},  \code{\link{sracipeKnockDown}},
#' \code{\link{sracipeOverExp}},  \code{\link{sracipePlotData}}
#'

sracipeSimulate <- function( circuit="inputs/test.tpo", config = config,
                      anneal=FALSE, knockOut = NA_character_,numModels=2000,
                      paramRange=100,
                      prodRateMin=1.,prodRateMax=100, degRateMin=0.1,
                      degRateMax=1.,foldChangeMin=1,foldChangeMax=100,
                      hillCoeffMin=1L,hillCoeffMax=6L,integrateStepSize=0.02,
                      simulationTime=50.0,nIC=1L,nNoise=0L,simDet = TRUE,
                      initialNoise=50.0,noiseScalingFactor=0.5,shotNoise=0,
                      scaledNoise=FALSE,outputPrecision=12L,
                      printStart = 50.0,
                      printInterval=10, stepper = "RK4",
                      thresholdModels = 5000, plots = FALSE, 
                      plotToFile = FALSE,
                      genIC = TRUE, genParams = TRUE,
                      integrate = TRUE, rkTolerance = 0.01, timeSeries = FALSE,
                      nCores = 1L,
                      ...) {
 rSet <- RacipeSE()
 metadataTmp <- metadata(rSet)
 configData <- NULL
 data("configData",envir = environment(), package = "sRACIPE")
 configuration <- configData
 

  if(methods::is(circuit,"RacipeSE"))
  {

      rSet <- circuit
      metadataTmp <- metadata(rSet)
      configuration <- metadata(rSet)$config
  }
   if((methods::is(circuit,"character")) | (methods::is(circuit, "data.frame")))
  {
    sracipeCircuit(rSet) <- circuit
    if(methods::is(circuit, "data.frame")){
      metadataTmp$annotation <- deparse(substitute(circuit))
      metadataTmp$nInteractions <- length(circuit[,1])
      }

   }

 if(missing(circuit)){
    message("Please specify a circuit either as a file or as a data.frame")
    return()
  }
if(!missing(config)){
 if(methods::is(config,"list")){
   configuration <- config
 }
  else{
    message("Incorrect config provided!")
    return()
  }
}

 if(!missing(knockOut)){
   if(!(methods::is(knockOut, "list"))){
     knockOut <- list(knockOut)
   }
 }
if(nCores<1) {
  warning("Number of cores, nCores, is less than 1 or not an interger. 
          Using nCores=1")
  nCores=1
  }
 if(!missing(anneal)){
   #print("test1")
   if(anneal)
     configuration$options["anneal"] <- anneal
 }
 if(!missing(numModels)){
   configuration$simParams["numModels"] <- numModels
 }
 
 if(!missing(paramRange)){
   configuration$simParams["paramRange"] <- paramRange
 }
 if(!missing(prodRateMin)){
   configuration$hyperParams["prodRateMin"] <- prodRateMin
 }
 if(!missing(prodRateMax)){
   configuration$hyperParams["prodRateMax"] <- prodRateMax
 }
 
 if(!missing(degRateMin)){
   configuration$hyperParams["degRateMin"] <- degRateMin
 }
 if(!missing(degRateMax)){
   configuration$hyperParams["degRateMax"] <- degRateMax
 }
 if(!missing(foldChangeMin)){
   configuration$hyperParams["foldChangeMin"] <- foldChangeMin
 }
 if(!missing(foldChangeMax)){
   configuration$hyperParams["foldChangeMax"] <- foldChangeMax
 }
 if(!missing(hillCoeffMin)){
   configuration$hyperParams["hillCoeffMin"] <- hillCoeffMin
 }
 if(!missing(hillCoeffMax)){
   configuration$hyperParams["hillCoeffMax"] <- hillCoeffMax
 }
 if(!missing(integrateStepSize)){
   configuration$simParams["integrateStepSize"] <- integrateStepSize
 }
 if(!missing(simulationTime)){
   configuration$simParams["simulationTime"] <- simulationTime
 }
 if(!missing(simDet)){
   configuration$options["simDet"] <- simDet
 } else {
   configuration$options["simDet"] <- TRUE
 }
 if(!missing(nIC)){
   configuration$simParams["nIC"] <- nIC
 }
 if(!missing(outputPrecision)){
   configuration$simParams["outputPrecision"] <- outputPrecision
 }
 
 if(!missing(nNoise)){
   configuration$stochParams["nNoise"] <- nNoise
 }
 if(!missing(initialNoise)){
   configuration$stochParams["initialNoise"] <- initialNoise
 }
 if(!missing(noiseScalingFactor)){
   configuration$stochParams["noiseScalingFactor"] <- noiseScalingFactor
 }
 if(!missing(shotNoise)){
   configuration$stochParams["shotNoise"] <- shotNoise
 }
 if(!missing(scaledNoise)){
   configuration$options["scaledNoise"] <-scaledNoise
 }
 if(!missing(printStart)){
   configuration$simParams["printStart"] <- printStart
 }
 if(!missing(printInterval)){
   configuration$simParams["printInterval"] <- printInterval
   if(configuration$simParams["printInterval"] <
      configuration$simParams["integrateStepSize"]){
     configuration$simParams["printInterval"] <-
       configuration$simParams["integrateStepSize"] 
     warnings("Print Interval cannot be smaller than integration step size. 
              Setting it to integrate step size.")}
   }
 
 if(!missing(genIC)){
   configuration$options["genIC"] <- genIC
 }
 
 if(!missing(genParams)){
   configuration$options["genParams"] <- genParams
 }
 
 if(!missing(integrate)){
   configuration$options["integrate"] <- integrate
 }
 if(missing(printStart)){
   configuration$simParams["printStart"] <-
     configuration$simParams["simulationTime"]
 }
 if(!missing(rkTolerance)){
   configuration$simParams["rkTolerance"] <- rkTolerance
 }
 
 # Apply parameter range
 configuration$hyperParams["prodRateMin"] <- 0.5*(
   configuration$hyperParams["prodRateMin"] +
     configuration$hyperParams["prodRateMax"]) - 0.5*(
       configuration$hyperParams["prodRateMax"] -
         configuration$hyperParams["prodRateMin"])*
   configuration$simParams["paramRange"]/100
 configuration$hyperParams["prodRateMax"] <- 0.5*(
   configuration$hyperParams["prodRateMin"] +
     configuration$hyperParams["prodRateMax"]) + 0.5*(
       configuration$hyperParams["prodRateMax"] -
         configuration$hyperParams["prodRateMin"])*
   configuration$simParams["paramRange"]/100
 
 configuration$hyperParams["degRateMin"] <- 0.5*(
   configuration$hyperParams["degRateMin"] +
     configuration$hyperParams["degRateMax"]) - 0.5*(
       configuration$hyperParams["degRateMax"] -
         configuration$hyperParams["degRateMin"])*
   configuration$simParams["paramRange"]/100
 configuration$hyperParams["degRateMax"] <- 0.5*(
   configuration$hyperParams["degRateMin"] +
     configuration$hyperParams["degRateMax"]) + 0.5*(
       configuration$hyperParams["degRateMax"] -
         configuration$hyperParams["degRateMin"])*
   configuration$simParams["paramRange"]/100
 
 configuration$hyperParams["foldChangeMin"] <- 0.5*(
   configuration$hyperParams["foldChangeMin"] +
     configuration$hyperParams["foldChangeMax"]) - 0.5*(
       configuration$hyperParams["foldChangeMax"] -
         configuration$hyperParams["foldChangeMin"])*
   configuration$simParams["paramRange"]/100
 configuration$hyperParams["foldChangeMax"] <- 0.5*(
   configuration$hyperParams["foldChangeMin"] +
     configuration$hyperParams["foldChangeMax"]) + 0.5*(
       configuration$hyperParams["foldChangeMax"] -
         configuration$hyperParams["foldChangeMin"])*
   configuration$simParams["paramRange"]/100
 

# stepper is not included in configdata. This can be changed
 configuration$stepper <- stepper


  nGenes <- length(rSet)
  geneInteraction <- as.matrix(rowData(rSet))

  thresholdGene <- rep(0, nGenes)

  geneNames <- names(rSet)
  # outFile <- paste0(Sys.Date(),"_",metadata(rSet)$annotation,"_",
  #                  basename(tempfile()))
  outFileGE <- tempfile(fileext = ".txt")
  outFileParams <- tempfile(fileext = ".txt")
  outFileIC <- tempfile(fileext = ".txt")
  if(genParams){
  message("Generating gene thresholds")
  #Rcpp::sourceCpp("src/thresholdGenerator.cpp")
  threshold_test <- generateThresholds(
    geneInteraction,  thresholdGene,  configuration)
  configuration$thresholds <- thresholdGene
  if(threshold_test!=0){
    stop("Error in threshold generation")
    }
  } else {

    if(is.null(colData(rSet))){
    message("Please specify the parameters as genParams is FALSE")
    return(rSet)
  } else {
    params <- as.data.frame(sracipeParams(rSet))

    utils::write.table(params, file = outFileParams,
                sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE)
  }
  }

    if(!genIC){
      if(is.null(colData(rSet))){
      message("Please specify the initial conditions
              as genIC is FALSE")
      return(rSet)
    } else {
      ic <- as.data.frame(t(sracipeIC(rSet)))
      utils::write.table(ic, file = outFileIC,
                  sep = "\t", quote = FALSE, row.names = FALSE,
                  col.names = FALSE)
    }
  }
if(missing(nNoise)){
  if(!missing(initialNoise) & !anneal )
  {
    if(timeSeries)
    # if(timeSeries) 
      {configuration$stochParams["nNoise"] <- 1}
    else {if(simDet) configuration$stochParams["nNoise"] <- 2}
  }
  
}
  if(anneal){
    if(missing(nNoise)) {configuration$stochParams["nNoise"] <- 30L}

  }
  # For time series, default to 1 model and single noise level.
  if(timeSeries){
    configuration$simParams["numModels"] <- 1
    # configuration$stochParams["nNoise"] <- 0
    
    if(missing(printInterval)) 
    configuration$simParams["printInterval"] <- max(
      0.05, configuration$simParams["integrateStepSize"])
    if(missing(printStart)) configuration$simParams["printStart"] <- 0
    
    if(configuration$stepper == "DP") {configuration$stepper <- "RK4"
    warnings("Using RK4 integration method for time series simulation.")
    
    }
  }
  stepperInt <- 1L
  if(configuration$stepper == "RK4"){ stepperInt <- 4L}
  if(configuration$stepper == "DP") {stepperInt <- 5L}
  
  if(configuration$stochParams["nNoise"] > 0) {
    if(stepper != "EM"){
      warnings("Defaulting to EM stepper for stochastic simulations")
      configuration$stepper <- "EM"
      stepperInt <- 1L
    }
    if(missing(initialNoise)) { configuration$stochParams["initialNoise"] <-
      as.numeric(50/sqrt(nGenes))}
  }





  annotationTmp <- outFileGE
#  message("Running the simulations")
  # print(configuration$stochParams["nNoise"])
  if(!configuration$options["integrate"] | (nCores==1) | timeSeries){
  Time_evolution_test<- simulateGRCCpp(geneInteraction, configuration,outFileGE,
                                       outFileParams,outFileIC, stepperInt)
  } else {
    if(nCores>1 & !timeSeries){
      configuration$options["integrate"] <- FALSE
    Time_evolution_test<- simulateGRCCpp(geneInteraction, configuration,outFileGE,
                                         outFileParams,outFileIC, stepperInt)
    configuration$options["integrate"] <- TRUE
    requireNamespace("doFuture")
    multiprocess <- NULL
    registerDoFuture()
    plan(multiprocess)
    
    configList <- list()
    parModel <- floor(configuration$simParams["numModels"]/nCores)
    gEFileList <- character(length = nCores)
    paramFileList <- character(length = nCores)
    iCFileList <- character(length = nCores)

    
    parameters <- utils::read.table(outFileParams, header = FALSE)
    ic <- utils::read.table(outFileIC, header = FALSE)
    
    for(i in seq(1,nCores)){
      parConfig <- configuration

      gEFileList[i] <- tempfile(fileext = ".txt")
      paramFileList[i] <- tempfile(fileext = ".txt")
      iCFileList[i] <- tempfile(fileext = ".txt")
      if(i==nCores){
        parConfig$simParams["numModels"] <- 
          configuration$simParams["numModels"] - (nCores-1)*parModel
        paramPar <- parameters[(((i-1)*parModel+1):(nrow(parameters))),]
        icPar <- ic[(((i-1)*parModel+1):(nrow(parameters))),]
        
      } else {
        parConfig$simParams["numModels"] <- parModel
        paramPar <- parameters[(((i-1)*parModel+1):(i*parModel)),]
        icPar <- ic[(((i-1)*parModel+1):(i*parModel)),]
        
      }
      
      configList[[i]] <- parConfig
      
      utils::write.table(paramPar, file = paramFileList[i],
                         sep = "\t", quote = FALSE, row.names = FALSE,
                         col.names = FALSE)
      utils::write.table(icPar, file = iCFileList[i],
                         sep = "\t", quote = FALSE, row.names = FALSE,
                         col.names = FALSE)
      
    }

    x <- foreach(configurationTmp = configList,outFileGETmp = gEFileList,
                 outFileParamsTmp=paramFileList, outFileICTmp=iCFileList,
                 .export = c("geneInteraction","stepperInt")) %dopar% {

            simulateGRCCpp(
              geneInteraction, configurationTmp,outFileGETmp, outFileParamsTmp,
              outFileICTmp, stepperInt)
            
    
                 }
    geList <- list()
    for(i in seq(1,nCores)){
      
      geList[[i]] <- utils::read.table(iCFileList[i], 
                                       header = FALSE)
    }
    geneExpression <- do.call(rbind, geList)
    
    }
  }
    if(configuration$options["integrate"]){

      
      if(timeSeries){

        timeStamps <- c(seq(
          configuration$simParams["printStart"] +
            configuration$simParams["printInterval"], 
          configuration$simParams["simulationTime"],
          configuration$simParams["printInterval"]),
          configuration$simParams["simulationTime"])
        geneExpression <- utils::read.table(outFileGE, header = FALSE)
        # print(dim(geneExpression))
        geneExpression <- matrix(geneExpression, ncol = nGenes, byrow = TRUE)
        # print(dim(geneExpression))
        parameters <- utils::read.table(outFileParams, header = FALSE)
        paramName <- sracipeGenParamNames(rSet)
        colnames(parameters) <- paramName
        ic <- utils::read.table(outFileIC, header = FALSE)
        if(configuration$stochParams["nNoise"] ==1){
          # print(dim(geneExpression))
          detGeneExp <- geneExpression[(1+dim(geneExpression)[1]/2):dim(geneExpression)[1],]
          geneExpression <- geneExpression[(seq_len(dim(geneExpression)[1]/2)),]
          #print(dim(geneExpression))
          detGeneExp <- t(detGeneExp)
          rownames(detGeneExp) <- geneNames
          colnames(detGeneExp) <- timeStamps[seq_len(dim(detGeneExp)[2])]
          metadataTmp$timeSeriesDet <- NULL
          metadataTmp$timeSeriesDet <- detGeneExp
          #print(dim(detGeneExp))
          
        }
        
        geneExpression <- t(geneExpression)
        #print(dim(geneExpression))
        colnames(ic) <- geneNames
        metadataTmp$normalized <- FALSE

        colData <- cbind(parameters,ic)
        metadataTmp$config <- configuration
      # colnames(geneExpression) <- seq(
      #  configuration$simParams["integrateStepSize"],
      #  configuration$simParams["simulationTime"],
      # configuration$simParams["integrateStepSize"])
        rownames(geneExpression) <- geneNames

        # print(dim(geneExpression))
        # print(length(timeStamps))
        colnames(geneExpression) <- timeStamps[seq_len(dim(geneExpression)[2])]
        metadataTmp$timeSeries <- NULL
      metadataTmp$timeSeries <- geneExpression

        rSet <- RacipeSE(rowData = geneInteraction, 
                         colData = colData,
                         metadata = metadataTmp)
        return(rSet)
      }
      if(!(nCores>1 & !timeSeries)) 
        geneExpression <- utils::read.table(outFileGE, header = FALSE)

    if(
      (configuration$simParams["printStart"] + 
       configuration$simParams["printInterval"]) < 
       configuration$simParams["simulationTime"]){

      timeStamps <- seq(
        configuration$simParams["printStart"] +
          configuration$simParams["printInterval"], 
        configuration$simParams["simulationTime"],
        configuration$simParams["printInterval"])

      tsSimulations <- as.list(timeStamps)
      names(timeStamps) <- timeStamps
      
      for(i in seq_along(timeStamps)){
        tsSimulations[[i]] <- geneExpression[
          , (1+(i-1)*(nGenes)):(i*nGenes)]
        colnames(tsSimulations[[i]]) <- geneNames
      }
      
      geneExpression <- geneExpression[
        ,(1+(length(timeStamps)*nGenes)):
          ((length(timeStamps)+1)*nGenes)]
      
      tsSimulations <- lapply(tsSimulations,t)
      
      geneExpression <- t(geneExpression)
      rownames(geneExpression) <- geneNames
      assayDataTmp <- c(list(deterministic = geneExpression),
                        tsSimulations)
      metadataTmp$tsSimulations <- timeStamps
      
      
    }
    if((configuration$stochParams["nNoise"] == 0) &  
       (configuration$simParams["printStart"] + 
        configuration$simParams["printInterval"] > 
       configuration$simParams["simulationTime"])){
      geneExpression <- t(geneExpression)
      rownames(geneExpression) <- geneNames
    assayDataTmp <- list(deterministic = geneExpression)}
    # colnames(geneExpression) <- geneNames
    if(configuration$stochParams["nNoise"] > 0){
      noiseLevels <- (configuration$stochParams["initialNoise"]*
                        configuration$stochParams["noiseScalingFactor"]^
                        seq(0,configuration$stochParams["nNoise"]-1))
      stochasticSimulations <- as.list(noiseLevels)
      names(stochasticSimulations) <- noiseLevels
      for(i in seq_len(configuration$stochParams["nNoise"])){
        stochasticSimulations[[i]] <- geneExpression[
          , (1+(i-1)*(nGenes)):(i*nGenes)]
        colnames(stochasticSimulations[[i]]) <- geneNames
      }
      stochasticSimulations <- lapply(stochasticSimulations,t)
      if(simDet){
      geneExpression <- geneExpression[
        ,(1+(configuration$stochParams["nNoise"]*nGenes)):
          ((configuration$stochParams["nNoise"]+1)*nGenes)]



      geneExpression <- t(geneExpression)
      rownames(geneExpression) <- geneNames
      assayDataTmp <- c(list(deterministic = geneExpression),
                        stochasticSimulations)
      } else {
        assayDataTmp <- stochasticSimulations
      }
      metadataTmp$stochasticSimulations <- noiseLevels


    }

    paramName <- sracipeGenParamNames(rSet)
    # paramFile <- paste0("tmp/",outFile,"_parameters.txt")
    if(!(nCores>1 & !timeSeries)) 
      parameters <- utils::read.table(outFileParams, header = FALSE)
    colnames(parameters) <- paramName

    # icFile <- paste0("tmp/",outFile,"_IC.txt")
    if(!(nCores>1 & !timeSeries)) ic <- utils::read.table(outFileIC, header = FALSE)
    colnames(ic) <- geneNames
    metadataTmp$normalized <- FALSE

    ## Knockouts

      if(!missing(knockOut)){
        if(knockOut[[1]] == "all"){
          knockOut <- as.list(geneNames)
        }
        knockOutData <- knockOut
        names(knockOutData) <- knockOut

        for(ko in seq_along(knockOut)){
          koGene <- knockOut[[ko]]
          knockOut_number <- which(koGene==geneNames)
          if(length(knockOut_number)==0){
            message("knockOut gene not found in the circuit")
            return(rSet)
          }
          params <- parameters
          params[,knockOut_number] <- 0.0
          
          icTmp <- ic
          icTmp[,knockOut_number] <- 0.0
          
          configTmp <- configuration
          configTmp$options["genIC"] <- FALSE
          configTmp$options["genParams"] <- FALSE

          utils::write.table(params, file = outFileParams,
                      sep = "\t", quote = FALSE, row.names = FALSE,
                      col.names = FALSE)
          utils::write.table(icTmp, file = outFileIC,
                      sep = "\t", quote = FALSE, row.names = FALSE,
                      col.names = FALSE)

          Time_evolution_test<- simulateGRCCpp(geneInteraction, configTmp,
            outFileGE,outFileParams,outFileIC, stepperInt)


          # geFile <- paste0("tmp/",outFileKO,"_geneExpression.txt")
          geneExpression <- utils::read.table(outFileGE, header = FALSE)
          colnames(geneExpression) <- geneNames
          knockOutData[[ko]] <- geneExpression

        }
        knockOutData <- lapply(knockOutData, t)
        assayDataTmp <- c(assayDataTmp,knockOutData)
        metadataTmp$knockOutSimulations <- names(knockOutData)

  }
    } else {
      paramName <- sracipeGenParamNames(rSet)
      if(!(nCores>1 & !timeSeries)) parameters <- utils::read.table(outFileParams, header = FALSE)
      colnames(parameters) <- paramName
      
      ic <- utils::read.table(outFileIC, header = FALSE)
      colnames(ic) <- geneNames
      colData <- (cbind(parameters,ic))
      metadataTmp$config <- configuration
      rSet <- RacipeSE(rowData = geneInteraction, colData = colData,
                       
                       metadata = metadataTmp)
      return(rSet)
  }
  
    colData <- (cbind(parameters,ic))
    metadataTmp$config <- configuration
 # return(list(rowData = geneInteraction, colData = colData,
   #         assays =  assayDataTmp,
  #          metadata = metadataTmp))
    rSet <- RacipeSE(rowData = geneInteraction, colData = colData,
                     assays =  assayDataTmp,
                     metadata = metadataTmp)
##############################
 if(plots){
    rSet <- sracipePlotData(rSet,plotToFile = plotToFile,...)
  }

return(rSet)
}