First, load the libraries needed for the workflow
library(tidyverse)
library(CoRC)Put in all variables needed to customise the workflow
Put in the path to the model (.cps), as local path (as string) and load the model. Use loadSBML() for sbml models.
modelPath <- "./Models/Aston2018Breakthrough_model.cps"
model <- loadModel(modelPath)Load the data into the variable data.
data <- read.csv("Data/Aston2018_data2.csv", sep = " ")Define experiment(s) for Parameter Estimation. For multiple experiments use list of defineExperiments().
pe_experiments <- defineExperiments(
data = data,
type = c("time", "dependent"),
mapping = c(NA, species("V", reference = "Concentration")),
weight_method = "mean_square"
)Use the addParameterEstimationParameter() function to add Parameter Estimation Parameters. Alternatively, use the add_pe_param_global() function defined in utils for adding global Quantities as parameters.
lapply(getGlobalQuantities()$key, add_pe_param_global, lower = 100, upper = 100)Set Parameter Estimation Settings.
setParameterEstimationSettings(
method = "LevenbergMarquardt"
)Set Sigma Point Method Settings. Choose alpha, beta, kappa, and the measurement error (var)
alpha <- 0.5
beta <- 2
kappa <- 3
measurement_error <- 0.1Do you want to save the result as a RDS object? Also give the file, where you want to save it.
save <- TRUE
file <- "docs/SP_Aston"This should generally not need to be changed, only if the workflow is to be changed more generally.
# Run Sigma Point Method
result <-
runSigmaPoint(
alpha = alpha,
beta = beta,
kappa = kappa,
var = measurement_error,
experiments = pe_experiments,
mean_fit_as_basis = TRUE
)
# Save result
if (save){
saveRDS(result,
file = file)
}This workflow produces the following output:
result$sp_mean
#> Values[s].InitialValue Values[r_T_Tmax].InitialValue Values[r_T_d].InitialValue
#> -3.806190e-02 2.577607e+04 2.931675e-03
#> Values[R].InitialValue Values[D].InitialValue Values[beta*].InitialValue
#> 2.399477e+05 5.612107e+02 7.715511e-09
#> Values[pstar].InitialValue Values[c].InitialValue
#> 8.629680e+01 3.484531e+00
result$sp_cov_matrix
#> Values[s].InitialValue Values[r_T_Tmax].InitialValue
#> Values[s].InitialValue 6.349302e-03 4.781556e+03
#> Values[r_T_Tmax].InitialValue 4.781556e+03 4.207596e+09
#> Values[r_T_d].InitialValue 7.937299e-06 6.785358e+00
#> Values[R].InitialValue -3.898922e+04 -3.091264e+10
#> Values[D].InitialValue 1.715471e+00 1.462854e+06
#> Values[beta*].InitialValue -6.248146e-11 -5.318492e-05
#> Values[pstar].InitialValue 7.139463e-01 6.068387e+05
#> Values[c].InitialValue -9.102983e-03 -7.702077e+03
#> Values[r_T_d].InitialValue Values[R].InitialValue
#> Values[s].InitialValue 7.937299e-06 -3.898922e+04
#> Values[r_T_Tmax].InitialValue 6.785358e+00 -3.091264e+10
#> Values[r_T_d].InitialValue 1.124024e-08 -5.248890e+01
#> Values[R].InitialValue -5.248890e+01 2.968827e+11
#> Values[D].InitialValue 2.418738e-03 -1.146623e+07
#> Values[beta*].InitialValue -8.602220e-14 4.019171e-04
#> Values[pstar].InitialValue 9.917123e-04 -4.685236e+06
#> Values[c].InitialValue -1.267008e-05 5.869745e+04
#> Values[D].InitialValue Values[beta*].InitialValue
#> Values[s].InitialValue 1.715471e+00 -6.248146e-11
#> Values[r_T_Tmax].InitialValue 1.462854e+06 -5.318492e-05
#> Values[r_T_d].InitialValue 2.418738e-03 -8.602220e-14
#> Values[R].InitialValue -1.146623e+07 4.019171e-04
#> Values[D].InitialValue 5.234561e+02 -1.860457e-08
#> Values[beta*].InitialValue -1.860457e-08 6.828131e-19
#> Values[pstar].InitialValue 2.142627e+02 -7.756204e-09
#> Values[c].InitialValue -2.742847e+00 9.887801e-11
#> Values[pstar].InitialValue Values[c].InitialValue
#> Values[s].InitialValue 7.139463e-01 -9.102983e-03
#> Values[r_T_Tmax].InitialValue 6.068387e+05 -7.702077e+03
#> Values[r_T_d].InitialValue 9.917123e-04 -1.267008e-05
#> Values[R].InitialValue -4.685236e+06 5.869745e+04
#> Values[D].InitialValue 2.142627e+02 -2.742847e+00
#> Values[beta*].InitialValue -7.756204e-09 9.887801e-11
#> Values[pstar].InitialValue 8.887328e+01 -1.122692e+00
#> Values[c].InitialValue -1.122692e+00 1.487077e-02Creates a function for easy Parameter Estimation Parameter Usage:
make_pe_param <- function(name, upper, lower) {
value_ref <- parameter_strict(name, reference = "Value")
value <- getValue(value_ref)
defineParameterEstimationParameter( ref = value_ref,
start_value = value,
lower_bound = value*lower,
upper_bound = value*upper)
}
make_pe_param_global <- function(name, upper, lower){
value <- getGlobalQuantities(key = name)$initial_value
value_ref <- getGlobalQuantityReferences(key = name)$initial_value
defineParameterEstimationParameter(ref = value_ref,
start_value = value,
lower_bound = value-value*lower,
upper_bound = value+value*upper)
}
add_pe_param_global <- function(name, upper, lower){
value <- getGlobalQuantities(key = name)$initial_value
value_ref <- getGlobalQuantityReferences(key = name)$initial_value
addParameterEstimationParameter( ref = value_ref,
start_value = value,
lower_bound = value-value*lower,
upper_bound = value+value*upper)
}