Trials
In Pumas QSP, the different variations of the model to be ran are called the trials. For example, one trial may specify that the model should be solved with a loading dose of 150mg, while the next trial specifies that the loading dose is 250mg. Each trial is then optionally tied to a dataset which, when defined in an inverse problem, specifies a multi-simulation optimization problem that the further functions (calibrate
, vpop
, etc.) generate solutions for. The type of trial is used to signify what the data corresponds to measuring, i.e. whether the trials is used to match data of time series, or steady states, etc.
Trial Types
The following describes the types of trials which can be generated.
PumasQSP.Trial
— TypeTrial(data, model; kwargs...)
The Trial
describes an experiment in which data
was obtained. The dynamics of the investigated system are represented in model
as an ODESystem. The trial is used within the optimization problem, as part of InverseProblem
to fit the unknown model
parameters and initial conditions to data
. In the case of Global Sensitivity Analysis (GSA), the data
is not needed and can be assigned to nothing
.
Required Keyword Arguments
tspan
, which indicates the timespan for which the model equations are solved.
Optional Keywords Arguments
u0
modifies the default values of the initial conditions of some or all of the states, e.g.u0 = [state_name => custom_initial_value]
params
modifies the default values of some or all the parameters, e.g.params = [p1 => specific_value, p2 => other_value]
err
specifies the trial's contribution to the cost function. Defaults to thel2loss
function. The function requires 2 arguments, the solution of the trial and the data and is expected to return a scalar value corresponding to the cost of the trial, i.e.err = (sol, data) -> compute_error
.likelihood
is used forMCMCOpt
methods. This likelihood describes the distribution that generated the data. By default, a multivariate Normal distribution is used, centered at the solution of the trial at each saved timepointu
with a standard deviations
around it, but different distributions can be used, e.g.,likelihood = (u, s) -> MvNormal(u, s)
.noise_priors
is used forMCMCOpt
methods. Describes the prior fors
. Thes
parameter represents any observation noise around the solution of the trial at each timepoint,u
. The default value for this prior distribution is anInverseGamma(2,3)
, but this can be changed, e.g.noise_prior = [s1 => Exponential(1), s2 => Gamma(3,2)]
). If only one prior is given, then it is assumed that it applies to all saved states (save_names
). Otherwise, a vector containing pairs of states and prior distributions with length equal tosaved_names
needs to be provided, to set different priors for each state. In the example above, the model had two statess1
ands2
so two pairs were provided. The noise parameterss
are assumed to be unique for each trial and to be constant across timepoints for each state of each trial.reduction
is used for GSA methods. The output ofreduction
is expected to be the quantity whose sensitivity is being investigated.doses
specifies the doses that occur in a trial. This argument can be equal to an instance ofBolus
,PeriodicBolus
,Infusion
,PeriodicInfusion
or aVector
of multiple instances of any of these dose types.save_names
is used to specify whichmodel
states are saved. The same states are extracted fromdata
.trial_name
, an identifying name. The default name is "Trial".forward_u0=true
, if the trial is part of a collection ofSteadyStateTrials
, thenforward_u0=true
signals that the trial should use the outcome of theSteadyStateTrial
of the same collection as its initial condition.
If additional keywords are passed, they will be forwarded to the solve
call. For example, one can pass alg=Tsit5()
to specify what solver will be used. More information about supported arguments can be found here.
PumasQSP.SteadyStateTrial
— TypeSteadyStateTrial(data, model; kwargs...)
Describes a trial that is ran until a steady state is reached. This object can be initialized in the same way as a Trial
object, with the only difference being that data
needs to be a Vector
here. The data
in this case represents the values of the saved states when the system has reached its steady state.
See the SciML documentation for background information on steady state problems.
PumasQSP.IndependentTrials
— TypeIndependentTrials(trials...)
This trial collection type indicates that each trial can be solved individualy and that there is no interaction between them. This trial type is automatically created it the trials are passed as a Vector
(i.e. [trial1, trial2])
PumasQSP.SteadyStateTrials
— TypeSteadyStateTrials(ss_trial, trials...; postprocess=last)
SteadyStateTrials
are a trial collection that describes a steady state trial (see SteadyStateTrial
) (specified as the first argument) followed by subsequent trials that can continue using the steady state by setting forward_u0=true
. The steady state solution that is passed on can be modified using the postprocess
keyword argument, which accepts a function with a single argument that represents the solution of the first trial and returns the state to be further passed on.
Simulation and Analysis Functions
To better understand and debug trials, the trials come with associated analysis functions to allow for easy investigation of the results in a trial-by-trial form. The following functions help the introspection of such trials.
JuliaSimModelOptimizer.simulate
— Functionsimulate(experiment::AbstractExperiment, prob::InverseProblem, x)
Simulate the given experiment
using optimization-state point x
, which contains values for each parameter and initial condition that is optimized in InverseProblem
prob
.
Loss Functions
By default, the loss function associated with a trial against its data is the standard Euclidian distance, also known as the L2 loss. However, PumasQSP provides alternative loss definitions to allow for customizing the fitting strategy.
JuliaSimModelOptimizer.l2loss
— Functionl2loss(sol, data)
Squared error loss :
$\sum_{i=1}^{M} \sum_{j=1}^{N} \left( \text{sol}_{i,j} - \text{data}_{i,j} \right)^2$
where N is the number of saved timepoints and M the number of measured states in the solution
JuliaSimModelOptimizer.ARMLoss
— FunctionARMLoss(sol, bounds)
Allen-Rieger-Musante (ARM) loss :
$\sum_{i=1}^{M} \sum_{j=1}^{N} \text{max} \left[ \left( \text{sol}_{i,j} - \frac{\text{u}_{i,j} + \text{l}_{i,j}}{2} \right)^2 - \left( \frac{\text{u}_{i,j} - \text{l}_{i,j}}{2} \right)^2, 0 \right]$
where N is the number of saved timepoints, M the number of measured states in the solution and l, u
are the lower and upper bounds of each measured state respectively.
Reference
Allen RJ, Rieger TR, Musante CJ. Efficient Generation and Selection of Virtual Populations in Quantitative Systems Pharmacology Models. CPT Pharmacometrics Syst Pharmacol. 2016 Mar;5(3):140-6. doi: 10.1002/psp4.12063. Epub 2016 Mar 17. PMID: 27069777; PMCID: PMC4809626.