Steady-state simulations

[PavelBal]

A function that runs a simulation for a sufficiently long time that can be considered as a "steady-state" and returns the end results of all state variables (molecules and parameters). The currently implemented behavior is described in the documentation:

#' Get the steady-state values of species and state variable parameters.
#'
#' @details The steady-state is considered to be the last values of the molecules
#' amounts and state variable parameters in the simulation with sufficiently long
#' simulation time, i.e., where the rates of
#'   the processes do not (significantly) change. The steady-state is NOT analytically
#'   calculated or estimated in any other way than simulating for the given time.
#'
#'
#' @param steadyStateTime Simulation time (minutes). In `NULL` (default), the
#' default simulation time is the start time of the last application plus three
#' days. The simulated time must be long enough for the system to reach a steady-state.
#' Either a single value (will be applied for all simulations), or a list of
#' values specific for each simulation. In latter case, must have equal size as
#' `simulations`. When providing a list, `NULL` is allowed to calculate the
#' time based on the last application.
#' @param quantitiesPaths List of quantity paths (molecules and/or parameters)
#'   for which the steady-state will be simulated. If `NULL` (default), all
#'   molecules and state variable parameters are considered. The same list is
#'   applied for all simulations.
#' @param simulations `Simulation` object or a list of `Simulation` objects
#' @param ignoreIfFormula If `TRUE` (default), species and parameters with
#'   initial values defined by a formula are not included.
#' @param stopIfNotFound Boolean. If `TRUE` (default), an error is thrown when
#'   results for certain species were not generated. This may happen when due to
#'   numerical problems some values cannot be calculated, though the whole
#'   simulation converges. Setting this argument to `FALSE` allows to ignore
#'   such errors. Check the outputs for empty values when using this option.
#' @param lowerThreshold Numerical value (in default unit of the output).
#' Any steady-state values below this value are considered as numerical noise
#' and replaced by 0. If `lowerThreshold` is `NULL`, no cut-off is applied.
#' Default value is 1e-15.
#' @param simulationRunOptions Optional instance of a `SimulationRunOptions`
#'  used during the simulation run.
#'
#' @return A named list, where the names are the IDs of the simulations and the
#'   entries are lists containing `paths` and their `values` at the end of the
#'   simulation.
#' @import ospsuite.utils
#' @export
#' @examples
#' simPath <- system.file("extdata", "Aciclovir.pkml", package = "ospsuite")
#' sim <- loadSimulation(simPath)
#' steadyState <- getSteadyState(simulations = sim)
#' # Set initial values for steady-state simulations
#' setQuantityValuesByPath(
#'   quantityPaths = steadyState[[sim$id]]$paths,
#'   values = steadyState[[sim$id]]$values, simulation = sim
#' )

We also have a function to export such a "steady state" to excel file that can be loaded in MoBi, would it be also of interest for the main package?

[PavelBal]

There is an ongoing discussion about how to deal with molecules/state variables which start values are defined by a formula.

Yuri05
I don't get this: how is a state variable supposed to be defined by explicit formula???

PavelBal
What we are dealing here with is the start value of the parameter, which can be an explicit formula.

Yuri05
ok... but I still don't understand the use case. So if a parameter is a state variable with initial value defined by formula - why would I need to ignore this parameter in the steady state?

PavelBal
I sometimes define the start value of the parameter as a formula that calculates the steady-state of this parameter based on other intputs. In this case, I would not like to overwrite the start value.

Yuri05
While I still don't understand the use case (maybe you can explain it to me during one of the meetings) - I wonder if we should set the default value of ignoreIfFormula to FALSE?
It seems pretty dangerous for me when the steady state of an ODE species is ignored because its initial value is not const. but defined by the formula.

PavelBal
I have been thinking about it for some time now, and I am still not convinced, maybe we should discuss it in a meeting.
The way I see it, if a start value of an ODE parameter (and the same goes for molecules, e.g. start values of proteins) is defined by a formula, then this formula represents the steady-state. Otherwise, there is something off with the equation, don't you think?

In the current proposal, start values will not be returned by default if they are defined by a formula, unless the user specifies the argument ignoreIfFormula to FALSE.

An example of initial values defined by a formula could be an RHS parameter that depends on other parameters (e.g., Stomach|Liquid, but also some effector parameter in a QSP model), but also the initial values of proteins (as they are calculated from expression profiles).

I vote for NOT RETURNING the values for these molecules/parameters by default, as the formula describing the initial value should represent the steady-state by design. Overwriting such a value removes the formula, so the user cannot rely on correct model behavior after applying such a steady-state and changing any parameter after. Overwriting a formula by a constant value is dangerous and I would not like to have this behavior by default.

@Yuri05 @msevestre

@StephanSchaller @wilbertdew @teutonicod @Sulav1 tagging you as you might have some opinion as a power user.