R/timeToEquilibrium.R

In tbuder/CellTrans: Quantifying stochastic cell state transitions

#' Calculation of time to equilibrium.
#' 
#' This function estimates the time from any intitial cell state proportions until the equilibrium proportions. The variable input contains the read data from the function readExperimentalData(). The variable initialDistribution is a vector describing the initial cell state proportions, for example c(0.25,0.25,0.25,0.25) for equal proportions of 4 cell types. The third parameter tol is a tolerance deviation between the cell state proportions of the equilibrium distribution and those of the predicted cell state proportions for each cell state. For the parameter tol, we recommend values between 0.01 and 0.02.
#' @param initialDistribution The starting proportions of the cell states as vector, e.g. c(1,0,..,0).
#' @param input A datalist generated by the function readExperimentaldata().
#' @param tol The tolerance deviation from the equilibrium distribution.
#' @keywords plot
#' @export



timeToEquilibrium <- function(input,initialDistribution,tol) {
  timepoints<-dlgList(title="Data point(s) for estimation",multiple=TRUE, choices=input$timepoints)$res  
  
  trMatrix=calculate_transitionMatrix(input$experimentalData,input$timepoints,timepoints)
  require(markovchain)
  MC <- new("markovchain", states = input$cell_types,
            transitionMatrix = trMatrix,
            name = "Transition matrix")
  equilibrium=steadyStates(MC)
  distribution=initialDistribution
  steps=0
  repeat {
  if (norm(equilibrium-distribution,"M")<tol) { print(paste("Expected time to equilibrium:",steps,input$timeunits))
                                                print("Predicted cell state proportions at this time:")
                                                print(distribution)
                                                break}
    

  
  if (steps==100000) { print(paste("Expected time to equilibrium could not be determined. Increasing the tolerance might help."))
                     break}

  distribution=distribution%*%trMatrix
  steps=steps+1
   }
}