Nothing
#' Simulation of enzyme evolution
#'
#' This function simulates evolution of enzyme concentrations under constraints
#'
#'
#' @details
#' Time step is defined between appearance of two mutation.
#'
#' Time step is simulated with function \code{\link{simul.next.resident}}, which gives values of resident after a time step.
#'
#' To reduce size of result matrix, certain results only are conserved.
#' State of simulation is observed every \code{pasobs} time steps.
#' There is \code{npt} observations, so result matrix has \code{npt} rows.
#' In total, a simulation includes \code{pasobs*npt} time steps. By default, there is \code{125 000} time steps.
#'
#' Chosen equilibrium for \code{pred_enzsim} and \code{pred_contsim} are the theoretical equilibrium for constraints \code{"SC"}, \code{"Comp"} and \code{"RegPos"},
#' and the effective one for constraints \code{"RegNeg"}, \code{"CRPos"} and \code{"CRNeg"}.
#' If there are regulation groups (\code{1<sum(1/B)<n}), the equilibrium is not computed, and \code{pred_} returns \code{NA}.
#'
#'
#'
#' @usage simul.evol.enz.one(E_ini_fun, kin_fun, Keq_fun, N_fun, correl_fun,
#' beta_fun=NULL, X_fun=1, pasobs=250, npt=500, max_mut_size_E=1, max_mut_size_A=1,
#' pmutA=0, typ_E=1, typ_A=1, use.old.mut=FALSE)
#'
#' @param pasobs Numeric. Number of time steps between two successive observations of the system. Default is \code{250}
#' @param npt Numeric. Number of observations. Default is \code{500}
#' @param E_ini_fun Numeric vector of the initial concentrations
#' @param kin_fun Numeric vector of the initial kinetic parameters
#' @inheritParams simul.next.resident
#'
#'
#'
#' @return Invisible list of 4 elements:
#' \itemize{
#' \item \code{$res_sim}: numeric matrix of \code{npt} rows and \code{3*n+3} columns.
#' Each row corresponds to state at each observation step (i.e. after \code{pasobs} mutations),
#' and columns are respectively concentrations (\code{1:n}), kinetic parameters (\code{n+1:2n}), total concentration (\code{2n+1}), total kinetic (\code{2n+2}), flux/fitness (\code{2n+3}), activities (\code{2n+4:3n+3});
#' \item \code{$pred_enzsim}: numeric matrix of \code{npt} rows and \code{n} columns, corresponding to relative concentrations at equilibrium, for each observation step (in rows);
#' \item \code{$pred_contsim}: same as \code{$pred_enzsim}, but for response coefficients
#' \item \code{$param}: list of input parameters: \itemize{
#' \item \code{n}: number of enzymes,
#' \item \code{E0}: numeric vector of initial concentrations,
#' \item \code{kin}: numeric vector of initial kinetic parameters,
#' \item \code{Keq}: numeric vector of constant equilibrium,
#' \item \code{beta}: matrix of co-regulation coefficients,
#' \item \code{B}: numeric vector of global co-regulation coefficients,
#' \item \code{correl}: character string indicating the constraint abbreviation,
#' \item \code{N}: population size,
#' \item \code{pasobs}: number of steps between two system observations,
#' \item \code{npt}: number of system observations,
#' \item \code{X}: parameter for flux computation,
#' \item \code{pmutA}: probability for activity mutation,
#' \item other input parameters}
#' }
#'
#' Note that \code{n} is the number of enzymes, which is the length of \code{E_ini_fun}.
#'
#'
#' @seealso
#' See function \code{\link{simul.next.resident}} to see how works each time steps.
#'
#' Use function \code{\link{simul.evol.enz.multiple}} to compute several simulations.
#'
#'
#'
#'
#'
#' @export
# Simulations of enzyme evolution
simul.evol.enz.one <- function(E_ini_fun, kin_fun, Keq_fun, N_fun, correl_fun, beta_fun=NULL,
X_fun=1, pasobs=250, npt=500, max_mut_size_E=1, max_mut_size_A=1,
pmutA=0, typ_E=1, typ_A=1, use.old.mut=FALSE) {
########### Settings
n_fun <- length(E_ini_fun)
if (length(beta_fun)==0) {
beta_fun <- diag(rep(1,n_fun))
}
B_fun <- apply(beta_fun,1,sumbis)
#Number of time steps/generations = observation steps (i.e. number of time step between two observations) x number of observations
ngene <- pasobs*npt
########## Test
is.correl.authorized(correl_fun)
is.beta.accurate(beta_fun,n_fun,correl_fun)
if (any(E_ini_fun<0)|any(kin_fun<0)|N_fun<0) {
stop("Enter only positive numeric for E_ini_fun, A_fun and N_fun")
}
if (pasobs<0|npt<0) {
stop("Simulations go only to the future. Enter positive value for pasobs and npt.")
}
########## Initialization
#Matrix for results of simulations
#Each row of resobs = one observation, after pasobs time steps
#E, kin, sum(E), sum(kin), J, A (npt lignes x 3*n+3 colonnes)
resobs<-matrix(0,nrow=npt,ncol=(3*n_fun+3))
#First resident
E_res <- E_ini_fun
kin_res <- kin_fun
#Activities and flux
A_fun <- activities(kin_fun,Keq_fun)
J_fun <- flux(E_ini_fun,A_fun,X_fun)
# First data
resobs[1,]<-c(E_ini_fun,kin_fun,sum(E_ini_fun),sum(kin_fun),J_fun,A_fun)
## Predictions for theoretical equilibrium
pred_enz<-matrix(NA,nrow=npt,ncol=n_fun) # for enzyme concentrations e* or tilde{e}
pred_cont<-matrix(NA,nrow=npt,ncol=n_fun) # for response coefficients R*
#First data
if (sum(1/B_fun)==1|sum(1/B_fun)==n_fun) {
#if ALL enzymes are co-regulated or independent
star <- predict_th(A_fun,correl_fun,B_fun)
pred_cont[1,] <- star$pred_r
if (correl_fun=="RegNeg"|correl_fun=="CRPos"|correl_fun=="CRNeg") {
star <- predict_eff(E_ini_fun,B_fun,A_fun,correl_fun)
pred_cont[1:npt,1:n_fun] <- 0 #null response coefficient
}
pred_enz[1,] <- star$pred_e
}
############## Simulation
#For each time step
for (l in (pasobs+1):ngene) {
#beginning from pasobs+1, because if l=pasobs, (l/pasobs)=1 and initial values are modified
#so there is (ngene-pasobs) mutation time steps
#Computes next resident values
next_gen <- simul.next.resident(E_res,kin_res,Keq_fun,N_fun,correl_fun,beta_fun,
X_fun,max_mut_size_E,max_mut_size_A,pmutA,typ_E,typ_A,use.old.mut)
A_next <- activities(next_gen$kin_next,Keq_fun)
#Set new resident values
E_res <- next_gen$E_next
kin_res <- next_gen$kin_next
###### Registration values every 'pasobs' time steps, named observation steps
# => accumulates mutations until (l/pasobs) is an integer
#test if it is an observation step, by verifying if (l/pasobs) is an integer
if(l/pasobs==floor(l/pasobs)) {
#register resident values
resobs[l/pasobs,] <- c(next_gen$E_next,next_gen$kin_next,next_gen$Etot_next,next_gen$kintot_next,next_gen$w_next,A_next)
#register resident predictions, when all enzymes are co-regulated of all independent
if (sum(1/B_fun)==1|sum(1/B_fun)==n_fun) {
star <- predict_th(A_next,correl_fun,B_fun)
pred_cont[l/pasobs,] <- star$pred_r
if (correl_fun=="RegNeg"|correl_fun=="CRPos"|correl_fun=="CRNeg") {
star <- predict_eff(E_ini_fun,B_fun,A_next,correl_fun)
}
pred_enz[l/pasobs,] <- star$pred_e
}
# Note: probably simpler but longer with rbind()...
}
}
return(invisible(list(res_sim=resobs,pred_enzsim=pred_enz,pred_contsim=pred_cont,
param=list(E0=E_ini_fun,kin=kin_fun,Keq=Keq_fun,n=n_fun,beta=beta_fun,B=B_fun,
correl=correl_fun,N=N_fun,pasobs=pasobs,npt=npt,
X=X_fun,max_mut_size_E=max_mut_size_E,max_mut_size_A=max_mut_size_A,
pmutA=pmutA,typmut_E=typ_E,typmut_A=typ_A,nsim=1))))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.