diff --git a/DESCRIPTION b/DESCRIPTION
index 751481e3c..21c7e12f6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -45,7 +45,7 @@ Suggests:
Language: en-US
Encoding: UTF-8
LazyData: true
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
Roxygen: list(markdown = TRUE)
VignetteBuilder: knitr
Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index 720b63dcd..439eeac79 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -87,6 +87,7 @@ export(getQuantity)
export(getQuantityValuesByPath)
export(getSimulationTree)
export(getStandardMoleculeParameters)
+export(getSteadyState)
export(getUnitsForDimension)
export(hasDimension)
export(hasUnit)
diff --git a/NEWS.md b/NEWS.md
index 96152ba74..26e504cc2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,17 @@
# ospsuite 12.0 (development)
+## New features
+
+- Added a function `getSteadyState()` to calculate steady state values for simulations.
+This function is of particular use for models of endogenous substrates, where changing
+a parameter value (e.g., the production rate) will change the steady-state values of
+the substrate.
+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.
+
## Breaking Changes
- The single argument of the `getBaseUnit()` function is now named
diff --git a/R/messages.R b/R/messages.R
index 65a368cd1..64fcc8dc7 100644
--- a/R/messages.R
+++ b/R/messages.R
@@ -85,3 +85,7 @@ messages$wrongUnitForQuantity <- function(quantityPath, unit, dimension) {
messages$invalidDataType <- function(name, dataType) {
paste0("Data type '", dataType, "' specified for data set '", name, "' is not valid. Valid data types are: 'simulated' or 'observed'.")
}
+
+messages$steadyStateTimeNotPositive <- function(value) {
+ paste0("The value of `steadyStateTime` must be > 0, but it is '", paste(value, collapse = ", "), "'")
+}
diff --git a/R/utilities-simulation.R b/R/utilities-simulation.R
index ff9edeaca..5f3291b8e 100644
--- a/R/utilities-simulation.R
+++ b/R/utilities-simulation.R
@@ -749,3 +749,247 @@ getSimulationTree <- function(simulationOrFilePath, quantityType = "Quantity") {
return(pathEnumList)
}
+
+#' 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 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
+#' )
+getSteadyState <- function(simulations,
+ quantitiesPaths = NULL,
+ steadyStateTime = NULL,
+ ignoreIfFormula = TRUE,
+ lowerThreshold = 1e-15,
+ simulationRunOptions = NULL) {
+ # Default time that is added to the time of the last administration for steady-state
+ DELTA_STEADY_STATE <- 3 * 24 * 60 # 3 days in minutes
+
+ ospsuite.utils::validateIsOfType(simulations, type = "Simulation")
+ ospsuite.utils::validateIsString(quantitiesPaths, nullAllowed = TRUE)
+ # Unlisting `steadyStateTime` because it can be a list of values, including NULL
+ ospsuite.utils::validateIsNumeric(unlist(steadyStateTime), nullAllowed = TRUE)
+ simulations <- ospsuite.utils::toList(simulations)
+
+ if (any(unlist(steadyStateTime) <= 0)) {
+ stop(messages$steadyStateTimeNotPositive(steadyStateTime))
+ }
+
+ # If `steadyStateTime` is a list of values, it must be of the same size as
+ # the list of simulations.
+ # Otherwise, repeat the value for the number of simulations
+ if (length(steadyStateTime) > 1) {
+ ospsuite.utils::validateIsSameLength(simulations, steadyStateTime)
+ } else {
+ steadyStateTime <- rep(steadyStateTime, length(simulations))
+ }
+
+ # First prepare all simulations by setting their outputs and time intervals
+ # If no quantities have been specified, the quantities paths may be different
+ # for each simulation and must be stored separately
+ simulationState <- .storeSimulationState(simulations)
+ quantitiesPathsMap <- vector(mode = "list", length = length(simulations))
+ for (idx in seq_along(simulations)) {
+ simulation <- simulations[[idx]]
+ simId <- simulation$id
+ # Set simulation time to the steady-state value.
+ # If the specified steady-state time is NULL, the simulation time is set to
+ # the time of the last application plus a specified delta.
+ clearOutputIntervals(simulation = simulation)
+ if (is.null(steadyStateTime[[idx]])) {
+ latestAdministration <- 0
+ # get the list of all administered molecules in the simulation
+ administeredMolecules <- simulation$allXenobioticFloatingMoleculeNames()
+ for (mol in administeredMolecules) {
+ # Iterate through all applications of each molecule
+ applications <- simulation$allApplicationsFor(mol)
+ for (app in applications) {
+ # if the time of the application is later than the latest administration
+ if (app$startTime$value > latestAdministration) {
+ # set the latest administration to the time of the application
+ latestAdministration <- app$startTime$value
+ }
+ }
+ }
+ # set the simulation time to the time of the last application plus the delta
+ simulation$outputSchema$addTimePoints(timePoints = latestAdministration + DELTA_STEADY_STATE)
+ } else {
+ simulation$outputSchema$addTimePoints(timePoints = steadyStateTime[[idx]])
+ }
+
+ # If no quantities are explicitly specified, simulate all outputs.
+ if (is.null(quantitiesPaths)) {
+ quantitiesPathsMap[[idx]] <- getAllStateVariablesPaths(simulation)
+ } else {
+ quantitiesPathsMap[[idx]] <- quantitiesPaths
+ }
+ names(quantitiesPathsMap)[[idx]] <- simId
+ setOutputs(quantitiesOrPaths = quantitiesPathsMap[[idx]], simulation = simulation)
+ }
+
+ # Run simulations concurrently
+ simulationResults <- runSimulations(
+ simulations = simulations,
+ simulationRunOptions = simulationRunOptions
+ )
+
+ # Iterate through simulations and get their outputs
+ outputMap <- vector(mode = "list", length = length(simulations))
+ for (idx in seq_along(simulations)) {
+ simulation <- simulations[[idx]]
+ simId <- simulation$id
+ simResults <- simulationResults[[simId]]
+
+ allOutputs <- getOutputValues(
+ simResults,
+ quantitiesOrPaths = quantitiesPathsMap[[simId]],
+ addMetaData = FALSE
+ )
+
+ # Get the end values of all outputs
+ endValues <- lapply(quantitiesPathsMap[[simId]], function(path) {
+ # Check if the quantity is defined by an explicit formula
+ isFormulaExplicit <- isExplicitFormulaByPath(
+ path = enc2utf8(path),
+ simulation = simulation
+ )
+
+ if (ignoreIfFormula && isFormulaExplicit) {
+ return(NULL)
+ }
+ value <- tail(allOutputs$data[path][[1]], 1)
+ # If the value is below the cut-off threshold, replace it by 0
+ if (!is.null(lowerThreshold) && value < lowerThreshold) {
+ value <- 0
+ }
+ return(value)
+ })
+
+ # Get the indices for which the outputs have been calculated
+ indices <- which(lengths(endValues) != 0)
+
+ # Reset simulation output intervals and output selections
+ .restoreSimulationState(simulations, simulationState)
+ outputMap[[idx]] <- list(paths = quantitiesPathsMap[[simId]][indices], values = endValues[indices])
+ names(outputMap)[[idx]] <- simId
+ }
+ return(outputMap)
+}
+
+#' Stores current simulation output state
+#'
+#' @description Stores simulation output intervals, output time points,
+#' and output selections in the current state.
+#'
+#' @param simulations List of `Simulation` objects
+#'
+#' @return A named list with entries `outputIntervals`, `timePoints`, and
+#' `outputSelections`. Every entry is a named list with names being the IDs
+#' of the simulations.
+#' @keywords internal
+#' @noRd
+.storeSimulationState <- function(simulations) {
+ simulations <- c(simulations)
+ # Create named vectors for the output intervals, time points, and output
+ # selections of the simulations in their initial state. Names are IDs of
+ # simulations.
+ oldOutputIntervals <-
+ oldTimePoints <-
+ oldOutputSelections <-
+ ids <- vector("list", length(simulations))
+
+ for (idx in seq_along(simulations)) {
+ simulation <- simulations[[idx]]
+ simId <- simulation$id
+ # Have to reset both the output intervals and the time points!
+ oldOutputIntervals[[idx]] <- simulation$outputSchema$intervals
+ oldTimePoints[[idx]] <- simulation$outputSchema$timePoints
+ oldOutputSelections[[idx]] <- simulation$outputSelections$allOutputs
+ ids[[idx]] <- simId
+ }
+ names(oldOutputIntervals) <-
+ names(oldTimePoints) <-
+ names(oldOutputSelections) <- ids
+
+ return(list(
+ outputIntervals = oldOutputIntervals,
+ timePoints = oldTimePoints,
+ outputSelections = oldOutputSelections
+ ))
+}
+
+
+#' Restore simulation output state
+#'
+#' @description Restores simulation output intervals, output time points,
+#' and output selections to the values stored in `simStateList`.
+#' @inheritParams .storeSimulationState
+#' @param simStateList Output of the function `.storeSimulationState`.
+#' A named list with entries `outputIntervals`, `timePoints`, and
+#' `outputSelections`. Every entry is a named list with names being the IDs of
+#' the simulations.
+#'
+#' @keywords internal
+#' @noRd
+.restoreSimulationState <- function(simulations, simStateList) {
+ simulations <- c(simulations)
+ for (simulation in simulations) {
+ simId <- simulation$id
+ # reset the output intervals
+ simulation$outputSchema$clear()
+ for (outputInterval in simStateList$outputIntervals[[simId]]) {
+ addOutputInterval(
+ simulation = simulation,
+ startTime = outputInterval$startTime$value,
+ endTime = outputInterval$endTime$value,
+ resolution = outputInterval$resolution$value
+ )
+ }
+ if (length(simStateList$timePoints[[simId]]) > 0) {
+ simulation$outputSchema$addTimePoints(simStateList$timePoints[[simId]])
+ }
+ # Reset output selections
+ clearOutputs(simulation)
+ for (outputSelection in simStateList$outputSelections[[simId]]) {
+ addOutputs(quantitiesOrPaths = outputSelection$path, simulation = simulation)
+ }
+ }
+}
diff --git a/man/getSteadyState.Rd b/man/getSteadyState.Rd
new file mode 100644
index 000000000..284d2e8d6
--- /dev/null
+++ b/man/getSteadyState.Rd
@@ -0,0 +1,67 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utilities-simulation.R
+\name{getSteadyState}
+\alias{getSteadyState}
+\title{Get the steady-state values of species and state variable parameters.}
+\usage{
+getSteadyState(
+ simulations,
+ quantitiesPaths = NULL,
+ steadyStateTime = NULL,
+ ignoreIfFormula = TRUE,
+ lowerThreshold = 1e-15,
+ simulationRunOptions = NULL
+)
+}
+\arguments{
+\item{simulations}{\code{Simulation} object or a list of \code{Simulation} objects}
+
+\item{quantitiesPaths}{List of quantity paths (molecules and/or parameters)
+for which the steady-state will be simulated. If \code{NULL} (default), all
+molecules and state variable parameters are considered. The same list is
+applied for all simulations.}
+
+\item{steadyStateTime}{Simulation time (minutes). In \code{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
+\code{simulations}. When providing a list, \code{NULL} is allowed to calculate the
+time based on the last application.}
+
+\item{ignoreIfFormula}{If \code{TRUE} (default), species and parameters with
+initial values defined by a formula are not included.}
+
+\item{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 \code{lowerThreshold} is \code{NULL}, no cut-off is applied.
+Default value is 1e-15.}
+
+\item{simulationRunOptions}{Optional instance of a \code{SimulationRunOptions}
+used during the simulation run.}
+}
+\value{
+A named list, where the names are the IDs of the simulations and the
+entries are lists containing \code{paths} and their \code{values} at the end of the
+simulation.
+}
+\description{
+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.
+}
+\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
+)
+}
diff --git a/man/messages.Rd b/man/messages.Rd
index 6af6c3d91..6b1d97d90 100644
--- a/man/messages.Rd
+++ b/man/messages.Rd
@@ -6,7 +6,7 @@
\title{List of functions and strings used to signal error messages
Extends the \code{messages} list from ospsuite.utils}
\format{
-An object of class \code{list} of length 53.
+An object of class \code{list} of length 54.
}
\usage{
messages
diff --git a/man/reexports.Rd b/man/reexports.Rd
index 3ce94825e..35ee57cd1 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -23,6 +23,6 @@ below to see their documentation.
\describe{
\item{ospsuite.utils}{\code{\link[ospsuite.utils:op-null-default]{\%||\%}}, \code{\link[ospsuite.utils]{enum}}, \code{\link[ospsuite.utils]{enumGetKey}}, \code{\link[ospsuite.utils]{enumGetValue}}, \code{\link[ospsuite.utils]{enumHasKey}}, \code{\link[ospsuite.utils]{enumKeys}}, \code{\link[ospsuite.utils]{enumPut}}, \code{\link[ospsuite.utils]{enumRemove}}, \code{\link[ospsuite.utils]{enumValues}}}
- \item{tlf}{\code{\link[tlf]{PlotGridConfiguration}}, \code{\link[tlf]{plotGrid}}}
+ \item{tlf}{\code{\link[tlf]{plotGrid}}, \code{\link[tlf]{PlotGridConfiguration}}}
}}
diff --git a/tests/testthat/_snaps/plot-observed-vs-simulated/custom.svg b/tests/testthat/_snaps/plot-observed-vs-simulated/custom.svg
index 33a9155e0..9731a30d1 100644
--- a/tests/testthat/_snaps/plot-observed-vs-simulated/custom.svg
+++ b/tests/testthat/_snaps/plot-observed-vs-simulated/custom.svg
@@ -191,19 +191,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -227,19 +214,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/tests/testthat/_snaps/plot-observed-vs-simulated/defaults-1-fold-dist.svg b/tests/testthat/_snaps/plot-observed-vs-simulated/defaults-1-fold-dist.svg
index b08a36663..ceb0a9dd8 100644
--- a/tests/testthat/_snaps/plot-observed-vs-simulated/defaults-1-fold-dist.svg
+++ b/tests/testthat/_snaps/plot-observed-vs-simulated/defaults-1-fold-dist.svg
@@ -189,19 +189,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -225,19 +212,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/tests/testthat/_snaps/plot-observed-vs-simulated/defaults-3-fold-dist.svg b/tests/testthat/_snaps/plot-observed-vs-simulated/defaults-3-fold-dist.svg
index a8a8eda1d..ae416135b 100644
--- a/tests/testthat/_snaps/plot-observed-vs-simulated/defaults-3-fold-dist.svg
+++ b/tests/testthat/_snaps/plot-observed-vs-simulated/defaults-3-fold-dist.svg
@@ -193,19 +193,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -229,19 +216,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/tests/testthat/_snaps/plot-observed-vs-simulated/defaults.svg b/tests/testthat/_snaps/plot-observed-vs-simulated/defaults.svg
index e28b02ca3..5138f1c7d 100644
--- a/tests/testthat/_snaps/plot-observed-vs-simulated/defaults.svg
+++ b/tests/testthat/_snaps/plot-observed-vs-simulated/defaults.svg
@@ -187,19 +187,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -223,19 +210,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/tests/testthat/_snaps/plot-observed-vs-simulated/linear-scale.svg b/tests/testthat/_snaps/plot-observed-vs-simulated/linear-scale.svg
index 894a07b1c..af9b9e937 100644
--- a/tests/testthat/_snaps/plot-observed-vs-simulated/linear-scale.svg
+++ b/tests/testthat/_snaps/plot-observed-vs-simulated/linear-scale.svg
@@ -192,20 +192,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -230,20 +216,6 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/tests/testthat/test-utilities-simulation.R b/tests/testthat/test-utilities-simulation.R
index 65bfa11a7..f92565ca5 100644
--- a/tests/testthat/test-utilities-simulation.R
+++ b/tests/testthat/test-utilities-simulation.R
@@ -449,3 +449,26 @@ test_that("it can explore a simulation by instance", {
path <- tree$Organism$Liver$Volume$path
expect_equal(path, "Organism|Liver|Volume")
})
+
+test_that("It calculates steady-state for multiple simulations, single steadyStateTime", {
+ simFilePath <- system.file("extdata", "Aciclovir.pkml", package = "ospsuite")
+ sim1 <- loadSimulation(simFilePath)
+ sim2 <- loadSimulation(simFilePath)
+ output <- getSteadyState(
+ simulations = c(sim1, sim2),
+ steadyStateTime = 1
+ )
+ expect_equal(names(output), c(sim1$id, sim2$id))
+})
+
+test_that("It calculates steady-state for multiple simulations, multiple steadyStateTime", {
+ simFilePath <- system.file("extdata", "Aciclovir.pkml", package = "ospsuite")
+ sim1 <- loadSimulation(simFilePath)
+ sim2 <- loadSimulation(simFilePath)
+ output <- getSteadyState(
+ simulations = c(sim1, sim2),
+ steadyStateTime = list(1, NULL)
+ )
+ expect_equal(names(output), c(sim1$id, sim2$id))
+ expect_equal(names(output[[sim1$id]]), c("paths", "values"))
+})