Introduction to Pumas

This is an introduction to Pumas, a software system for pharmaceutical modeling and simulation. Pumas supports NLME (non-linear mixed effects models) estimation and simulation, NCA (non-compartmental analysis), Bioequivalence analysis, optimal design, and more.

In this tutorial we focus on NLME. The basic workflow of NLME with Pumas is:

  1. Build a model.
  2. Define subjects or populations to simulate or estimate.
  3. Analyze the results with post-processing and plots.

We will show how to build a multiple-response PK/PD model via the @model macro, define a subject with single doses, and analyze the results of the simulation. Fitting of data using any of the methods available in Pumas is showcased in the tutorials. This tutorial is designed as a broad overview of the workflow and a more in-depth treatment of each section can be found in the tutorials.

Working Example

Let us start by showing complete simulation code, and then break down how it works.

<!– TODO: It would be useful to include a figure that illustrates the model, with the parameter names identified on the figure. –>

using Pumas
using PumasUtilities
using CairoMakie
######### Turnover model

inf_2cmt_lin_turnover = @model begin
    @param begin
        tvcl ∈ RealDomain(; lower = 0)
        tvvc ∈ RealDomain(; lower = 0)
        tvq ∈ RealDomain(; lower = 0)
        tvvp ∈ RealDomain(; lower = 0)
        Ω_pk ∈ PDiagDomain(4)
        σ_prop_pk ∈ RealDomain(; lower = 0)
        σ_add_pk ∈ RealDomain(; lower = 0)
        # PD parameters
        tvturn ∈ RealDomain(; lower = 0)
        tvebase ∈ RealDomain(; lower = 0)
        tvec50 ∈ RealDomain(; lower = 0)
        Ω_pd ∈ PDiagDomain(1)
        σ_add_pd ∈ RealDomain(; lower = 0)
    end

    @random begin
        ηpk ~ MvNormal(Ω_pk)
        ηpd ~ MvNormal(Ω_pd)
    end

    @pre begin
        CL = tvcl * exp(ηpk[1])
        Vc = tvvc * exp(ηpk[2])
        Q = tvq * exp(ηpk[3])
        Vp = tvvp * exp(ηpk[4])

        ebase = tvebase * exp(ηpd[1])
        ec50 = tvec50
        emax = 1
        turn = tvturn
        kout = 1 / turn
        kin0 = ebase * kout
    end

    @init begin
        Resp = ebase
    end

    @vars begin
        conc := Central / Vc
        edrug := emax * conc / (ec50 + conc)
        kin := kin0 * (1 - edrug)
    end

    @dynamics begin
        Central' = -(CL / Vc) * Central + (Q / Vp) * Peripheral - (Q / Vc) * Central
        Peripheral' = (Q / Vc) * Central - (Q / Vp) * Peripheral
        Resp' = kin - kout * Resp
    end

    @derived begin
        dv ~ @. Normal(conc, sqrt((conc * σ_prop_pk)^2 + σ_add_pk^2))
        resp ~ @. Normal(Resp, σ_add_pd)
    end
end

turnover_params = (;
    tvcl = 1.5,
    tvvc = 25.0,
    tvq = 5.0,
    tvvp = 150.0,
    tvturn = 10,
    tvebase = 10,
    tvec50 = 0.3,
    Ω_pk = Diagonal([0.05, 0.05, 0.05, 0.05]),
    Ω_pd = Diagonal([0.05]),
    σ_prop_pk = 0.15,
    σ_add_pk = 0.5,
    σ_add_pd = 0.2,
);
regimen = DosageRegimen(150; rate = 10, cmt = :Central)
pop = map(i -> Subject(; id = i, events = regimen), 1:10)
Population
  Subjects: 10
  Observations: 
sd_obstimes = [
    # single dose observation times
    0.25,
    0.5,
    0.75,
    1,
    2,
    4,
    8,
    12,
    16,
    20,
    24,
    36,
    48,
    60,
    71.9,
];
sims = simobs(inf_2cmt_lin_turnover, pop, turnover_params; obstimes = sd_obstimes)
Simulated population (Vector{<:Subject})
  Simulated subjects: 10
  Simulated variables: dv, resp
sim_plot(
    inf_2cmt_lin_turnover,
    sims;
    observations = [:dv],
    figure = (; fontsize = 18),
    axis = (; xlabel = "Time (hr)", ylabel = "Concentration (ng/mL)"),
)
Example block output
sim_plot(
    inf_2cmt_lin_turnover,
    sims;
    observations = [:resp],
    figure = (; fontsize = 18),
    axis = (; xlabel = "Time (hr)", ylabel = "Response"),
)
Example block output

In this code, we defined a nonlinear mixed effects model by describing the parameters, the random effects, the dynamical model, and the derived (result) values. Then we generated a population of 10 subjects who received a single dose of 150 mg, specified parameter values, simulated the model, and generated a plot of the results. Now let's walk through this process!

Using the Model Macro

First we define the model. The simplest way to do this is via @model. Inside the @model block we have a few subsections. The first of these is a @param block.

Pumas.@paramMacro
@param

Defines the model's fixed effects and corresponding domains, e.g. tvcl ∈ RealDomain(lower = 0). Must be used in an @model block. For example:

@model begin
  @param begin
    tvcl ∈ RealDomain(lower = 0)
    tvv ∈ RealDomain(lower = 0)
    Ω ∈ PDiagDomain(2)
    σ_prop ∈ RealDomain(lower = 0, init = 0.04)
  end
end

For PumasEMModels, the @param block uses a formula syntax for covariate relationships. A distribution from the normal family is used to indicate domain (currently only Normal, LogNormal, and LogitNormal are supported):

@emmodel begin
  @param begin
    CL ~ 1 + wt | LogNormal
    Vc ~ 1 | Normal
  end
end

this gives CL = exp(log(θ.CL[1]) + θ.CL[2] * wt) and VC = θ.VC, where θ is the parameter named tuple. If we only have 1 in the formula, then it is a scalar. If we have covariates, we have a tuple of parameters for variable.

For this model we define structural model parameters of PK and PD and their corresponding variances where applicable:

@param begin
    tvcl ∈ RealDomain(; lower = 0)
    tvvc ∈ RealDomain(; lower = 0)
    tvq ∈ RealDomain(; lower = 0)
    tvvp ∈ RealDomain(; lower = 0)
    Ω_pk ∈ PDiagDomain(4)
    σ_prop_pk ∈ RealDomain(; lower = 0)
    σ_add_pk ∈ RealDomain(; lower = 0)
    # PD parameters
    tvturn ∈ RealDomain(; lower = 0)
    tvebase ∈ RealDomain(; lower = 0)
    tvec50 ∈ RealDomain(; lower = 0)
    Ω_pd ∈ PDiagDomain(1)
    σ_add_pd ∈ RealDomain(; lower = 0)
end

Note the mathematical elegence of the expression of parameters, relationships, and distributions.

<!– NOTE: Since this is labeled an introduction and a tutorial (even though it's in docs), it would be good to explain, or link to the definition and use of RealDomain and PDiagDomain. Also, to point to some instructions on how to type the 'is an element of' symbol. Also capital omega. These things are very new to non-Pumas/Julia users. –> Next we define our random effects using the @random block.

Pumas.@randomMacro
@random

defines the model's random effects and corresponding distribution, e.g. η ~ MvNormal(Ω). Must be used in an @model block. For example:

@model begin
  @random begin
    η ~ MvNormal(Ω)
  end
end

For PumasEMModels, the @random block is equivalent to the @params block, except that each variable defined here has an additional random effect. For example:

@emmodel begin
  @random begin
    CL ~ 1 + wt | LogNormal
    Vc ~ 1 | Normal
  end
end

this gives CL = exp(log(θ.CL[1]) + θ.CL[2] * wt + η_CL) and VC = θ.VC + η_VC, where θ is the parameter named tuple. If we only have 1 in the formula, then it is a scalar. If we have covariates, we have a tuple of parameters for variable.

The random effects are defined by a distribution from the Distributions package. For more information on defining distributions, please see the Distributions documentation. For this tutorial, we wish to have a multivariate normal of uncorrelated random effects, one for PK and another for PD, so we utilize the syntax:

@random begin
    ηpk ~ MvNormal(Ω_pk)
    ηpd ~ MvNormal(Ω_pd)
end

Now we define our pre-processing step in @pre.

Pumas.@preMacro
@pre

Pre-processes variables for the dynamic system and statistical specification. Must be used in an @model block. For example:

@model begin
  @params begin
    tvcl ∈ RealDomain(lower = 0)
    tvv ∈ RealDomain(lower = 0)
    ωCL ∈ RealDomain(lower = 0)
    ωV ∈ RealDomain(lower = 0)
    σ_prop ∈ RealDomain(lower = 0, init = 0.04)
  end
  @random begin
    ηCL ~ Normal(ωCL)
    ηV ~ Normal(ωV)
  end
  @covariates wt
  @pre begin
    CL = tvcl * (wt / 70)^0.75 * exp(ηCL)
    Vc = tvv * (wt / 70) * exp(ηV)
    ka = CL / Vc
    Q = Vc
  end
end

For our PK/PD model, we define the individual parameters and their distributions, and we give them names as follows:

@pre begin
    CL = tvcl * exp(ηpk[1])
    Vc = tvvc * exp(ηpk[2])
    Q = tvq * exp(ηpk[3])
    Vp = tvvp * exp(ηpk[4])

    ebase = tvebase * exp(ηpd[1])
    ec50 = tvec50
    emax = 1
    turn = tvturn
    kout = 1 / turn
    kin0 = ebase * kout
end

Next we define the @init block which specifies the initial conditions of the dynamic system (i.e., initial compartment amounts).

Pumas.@initMacro
@init

Defines initial conditions of the dynamic system. Must be used in an @model block. For example:

@model begin
  @init begin
    Depot1 = 0.0
    Depot2 = 0.0
    Central = 0.0
  end
end

Any variable not mentioned in this block is assumed to have a zero for its starting value. We wish to only set the starting value for Resp, and thus we use:

@init begin
    Resp = ebase
end

Now we define our dynamics. We do this via the @dynamics block.

Pumas.@dynamicsMacro
@dynamics

Defines the dynamic system by specifying differential equations. Must be used in an @model block.

Note

A @model block with a @dynamics block must not contain a @reactions block.

Examples

A @dynamics block with a system of differential equations:

@model begin
  @dynamics begin
    Depot1' = -ka * Depot1
    Depot2' = ka * Depot1 - ka * Depot2
    Central' = ka * Depot2 - CL / Vc * Central
  end
end

A @dynamics block with a differential equation model:

@model begin
  @dynamics Depots1Central1
end

For our PK/PD model, we specify the three differential equations that describe the PK and PD components of the system:

@dynamics begin
    Central' = -(CL / Vc) * Central + (Q / Vp) * Peripheral - (Q / Vc) * Central
    Peripheral' = (Q / Vc) * Central - (Q / Vp) * Peripheral
    Resp' = kin - kout * Resp
end

Note the convenience of expressing the compartment amounts by name and the value of the derivative equations using the single quote character.

Next we set up alias variables that can be used later in the code. Such alias code can be setup in the @vars block.

Pumas.@varsMacro
@vars

Define variables usable in other blocks. It is recommended to define them in @pre instead if not a function of dynamic variables. Must be used in an @model block. For example:

@model begin
  @vars begin
    conc = Central / Vc
  end
end

We create aliases for the concentration and response as follows:

@vars begin
    conc := Central / Vc
    edrug := emax * conc / (ec50 + conc)
    kin := kin0 * (1 - edrug)
end

Note that the local assignment := can be used to define intermediate statements that will not be carried outside the block. This means that all the resulting data workflows from this model will not contain the intermediate variables defined with :=. This optional use of local assignment improves computational efficiency.

Note that without adjustment the units of conc are the units of the dose divided by the units of the volume parameter. The balancing of conc units with those reported from the assay could occur here, but in this example we make this adjustment below in the @derived block.

Lastly we utilize the @derived block.

Pumas.@derivedMacro
@derived

Defines the error model of the dependent variables. Must be used in an @model block. For example:

@model begin
  @derived begin
    dv ~ @. Normal(1000 * Central / Vc, 1000 * Central / Vc * σ)
  end
end

Our error model includes a component for each type of observation (i.e., conc and Resp):

@derived begin
    dv ~ @. Normal(conc, sqrt((conc * σ_prop_pk)^2 + σ_add_pk^2))
    resp ~ @. Normal(Resp, σ_add_pd)
end

Here dv (i.e., dependent variable) is the drug concentration and resp is the response. At any time that these outputs are sampled dv is normally distributed with mean conc and variance (conc * σ_prop_pk)^2 + σ_add_pk^2, and resp has mean Resp and a constant standard deviation σ_add_pd.

Building a Subject

Now let's build a subject to simulate the model with.

Pumas.SubjectType
Subject

The data corresponding to a single subject:

Fields:

  • id: identifier
  • observations: a named tuple of the dependent variables
  • covariates: a named tuple containing the covariates, or nothing.
  • events: a DosageRegimen, or nothing.
  • time: a vector of time stamps for the observations

When there are time varying covariates, each covariate is interpolated with a common covariate time support. The interpolated values are then used to build a multi-valued interpolant for the complete time support.

From the multi-valued interpolant, certain discontinuities are flagged in order to use that information for the differential equation solvers and to correctly apply the analytical solution per region as applicable.

Constructor

Subject(;id = "1",
         observations = NamedTuple(),
         events = nothing,
         time = observations isa AbstractDataFrame ? observations.time : nothing,
         covariates::Union{Nothing, NamedTuple} = nothing,
         covariates_time = observations isa AbstractDataFrame ? observations.time : nothing,
         covariates_direction = :left)

Subject may be constructed from an <:AbstractDataFrame with the appropriate schema or by providing the arguments directly through separate DataFrames / structures.

Examples:

julia> Subject()
Subject
  ID: 1

julia> Subject(
         id = 20,
         events = DosageRegimen(200, ii = 24, addl = 2),
         covariates = (WT = 14.2, HT = 5.2),
       )
Subject
  ID: 20
  Events: 3
  Covariates: WT, HT

julia> Subject(covariates = (WT = [14.2, 14.7], HT = fill(5.2, 2)), covariates_time = [0, 3])
Subject
  ID: 1
  Covariates: WT, HT

Our model did not make use of covariates, so we will ignore covariate components covariates, covariate_time and covariate_direction for now, and observations are only necessary for fitting parameters to data which will not be covered in this tutorial. Thus, our subject will be defined simply by its dosage regimen.

To do this, we use the DosageRegimen constructor.

Pumas.DosageRegimenType
DosageRegimen(
  amt::Numeric;
  time::Numeric = 0,
  cmt::Union{Numeric,Symbol} = 1,
  evid::Numeric = 1,
  ii::Numeric = zero.(time),
  addl::Numeric = 0,
  rate::Numeric = zero.(amt)./oneunit.(time),
  duration::Numeric = zero(amt)./oneunit.(time),
  ss::Numeric = 0,
  route::NCA.Route = NCA.NullRoute,
)

Lazy representation of a series of events.

Examples

julia> DosageRegimen(100; ii = 24, addl = 6)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss    route
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8  NCA.Route
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1     24.0      6      0.0       0.0     0  NullRoute

julia> DosageRegimen(50; ii = 12, addl = 13)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss    route
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8  NCA.Route
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │     0.0      1     50.0     1     12.0     13      0.0       0.0     0  NullRoute

julia> DosageRegimen(200; ii = 24, addl = 2)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss    route
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8  NCA.Route
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    200.0     1     24.0      2      0.0       0.0     0  NullRoute

julia> DosageRegimen(200; ii = 24, addl = 2, route = NCA.IVBolus)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss    route
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8  NCA.Route
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    200.0     1     24.0      2      0.0       0.0     0  IVBolus

You can also create a new DosageRegimen from various existing DosageRegimens:

evs = DosageRegimen(
  regimen1::DosageRegimen,
  regimen2::DosageRegimen;
  offset = nothing
)

offset specifies if regimen2 should start after an offset following the end of the last event in regimen1.

Examples

julia> e1 = DosageRegimen(100; ii = 24, addl = 6)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss    route
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8  NCA.Route
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1     24.0      6      0.0       0.0     0  NullRoute

julia> e2 = DosageRegimen(50; ii = 12, addl = 13)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss    route
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8  NCA.Route
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │     0.0      1     50.0     1     12.0     13      0.0       0.0     0  NullRoute

julia> evs = DosageRegimen(e1, e2)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss    route
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8  NCA.Route
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1     24.0      6      0.0       0.0     0  NullRoute
   2 │     0.0      1     50.0     1     12.0     13      0.0       0.0     0  NullRoute

julia> DosageRegimen(e1, e2; offset = 10)
DosageRegimen
 Row │ time     cmt    amt      evid  ii       addl   rate     duration  ss    route
     │ Float64  Int64  Float64  Int8  Float64  Int64  Float64  Float64   Int8  NCA.Route
─────┼───────────────────────────────────────────────────────────────────────────────────
   1 │     0.0      1    100.0     1     24.0      6      0.0       0.0     0  NullRoute
   2 │   178.0      1     50.0     1     12.0     13      0.0       0.0     0  NullRoute

For our example, let's assume the dose is a 150 mg infusion administered at the rate of 10 mg/hour into the Central compartment:

regimen = DosageRegimen(150; time = 0, rate = 10, cmt = :Central)
1×10 DataFrame
Rowtimecmtamtevidiiaddlratedurationssroute
Float64SymbolFloat64Int8Float64Int64Float64Float64Int8NCA.Route
10.0Central150.010.0010.015.00NullRoute

Let's define our subject to have id=1 with this single-dose regimen:

subject = Subject(; id = 1, events = regimen)
Subject
  ID: 1
  Events: 1

We can also create a collection of subjects as a Population.

pop = Population(map(i -> Subject(; id = i, events = regimen), 1:10))
Population
  Subjects: 10
  Observations: 

Running a Simulation

The main function for running a simulation is simobs.

Pumas.simobsFunction
simobs(
  model::AbstractPumasModel,
  population::Union{Subject, Population},
  param,
  randeffs=sample_randeffs(model, param, population);
  obstimes=nothing,
  ensemblealg=EnsembleSerial(),
  diffeq_options=NamedTuple(),
  rng=Random.default_rng(),
  simulate_error=Val(true),
)

Simulate random observations from model for population with parameters param at obstimes (by default, use the times of the existing observations for each subject). If no randeffs is provided, then random ones are generated according to the distribution in the model.

Arguments

  • model may either be a PumasModel or a PumasEMModel.

  • population may either be a Population of Subjects or a single Subject.

  • param should be either a single parameter set, in the form of a NamedTuple, or a vector of such parameter sets. If a vector then each of the parameter sets in that vector will be applied in turn. Example: (; tvCL=1., tvKa=2., tvV=1.)

  • randeffs is an optional argument that, if used, should specify the random effects for each subject in population. Such random effects are specified by NamedTuples for PumasModels (e.g. (; tvCL=1., tvKa=2.)) and by Tuples for PumasEMModels (e.g. (1., 2.)). If population is a single Subject (without being enclosed in a vector) then a single such random effect specifier should be passed. If, however, population is a Population of multiple Subjects then randeffs should be a vector containing one such specifier for each Subject. The functions init_randeffs, zero_randeffs, and sample_randeffs are sometimes convenient for generating randeffs:

randeffs = zero_randeffs(model, param, population)
solve(model, population, param, randeffs)

If no randeffs is provided, then random ones are generated according to the distribution in the model.

  • obstimes is a keyword argument where you can pass a vector of times at which to simulate observations. The default, nothing, ensures the use of the existing observation times for each Subject.

  • ensemblealg is a keyword argument that allows you to choose between different modes of parallelization. Options include EnsembleSerial(), EnsembleThreads() and EnsembleDistributed().

  • diffeq_options is a keyword argument that takes a NamedTuple of options to pass on to the differential equation solver.

  • rng is a keyword argument where you can specify which random number generator to use.

  • simulate_error is a keyword argument that allows you to choose whether a sample (Val(true) or true) or the mean (Val(false) or false) of the predictive distribution's error model will be returned for each subject.

simobs(model::PumasModel, population::Population; samples::Int = 10, simulate_error = true, rng::AbstractRNG = default_rng())

Simulates model predictions for population using the prior distributions for the population and subject-specific parameters. samples is the number of simulated predictions returned per subject. If simulate_error is false, the mean of the predictive distribution's error model will be returned, otherwise a sample from the predictive distribution's error model will be returned for each subject. rng is the pseudo-random number generator used in the sampling.

simobs(model::PumasModel, subject::Subject; samples::Int = 10, simulate_error = true, rng::AbstractRNG = default_rng())

Simulates model predictions for subject using the prior distributions for the population and subject-specific parameters. samples is the number of simulated predictions returned. If simulate_error is false, the mean of the predictive distribution's error model will be returned, otherwise a sample from the predictive distribution's error model will be returned for each subject. rng is the pseudo-random number generator used in the sampling.

simobs(b::BayesMCMCResults; subject = nothing, samples::Int = 10, simulate_error = false, rng::AbstractRNG = default_rng())

Simulates model predictions using the posterior samples for the population and subject-specific parameters. The predictions are either for the subject with index subject or the whole population if no index is input. samples is the number of simulated predictions returned per subject. If simulate_error is false, the mean of the predictive distribution's error model will be returned, otherwise a sample from the predictive distribution's error model will be returned for each subject. rng is the pseudo-random number generator used in the sampling.

simobs(pmi::FittedPumasModelInference, population::Population = pmi.fpm.data, randeffs::Union{Nothing, AbstractVector{<:NamedTuple}} = nothing; samples::Int, rng = default_rng(), kwargs...,)

Simulates observations from the fitted model using a truncated multivariate normal distribution for the parameter values when a covariance estimate is available and using non-parametric sampling for simulated inference and bootstrapping. When a covariance estimate is available, the optimal parameter values are used as the mean and the variance-covariance estimate is used as the covariance matrix. Rejection sampling is used to avoid parameter values outside the parameter domains. Each sample uses a different parameter value. samples is the number of samples to sample. population is the population of subjects which defaults to the population associated the fitted model. randeffs can be set to a vector of named tuples, one for each sample. If randeffs is not specified (the default behaviour), it will be sampled from its distribution.

simobs(fpm::FittedPumasModel, [population::Population,] vcov::AbstractMatrix, randeffs::Union{Nothing, AbstractVector{<:NamedTuple}} = nothing; samples::Int, rng = default_rng(), kwargs...)

Simulates observations from the fitted model using a truncated multi-variate normal distribution for the parameter values. The optimal parameter values are used for the mean and the user supplied variance-covariance (vcov) is used as the covariance matrix. Rejection sampling is used to avoid parameter values outside the parameter domains. Each sample uses a different parameter value. samples is the number of samples to sample. population is the population of subjects which defaults to the population associated the fitted model so it's optional to pass. randeffs can be set to a vector of named tuples, one for each sample. If randeffs is not specified (the default behaviour), it will be sampled from its distribution.

Let's simulate subject 1 with a set of chosen parameters:

turnover_params = (;
    tvcl = 1.5,
    tvvc = 25.0,
    tvq = 5.0,
    tvvp = 150.0,
    tvturn = 10,
    tvebase = 10,
    tvec50 = 0.3,
    Ω_pk = Diagonal([0.05, 0.05, 0.05, 0.05]),
    Ω_pd = Diagonal([0.05]),
    σ_prop_pk = 0.15,
    σ_add_pk = 0.5,
    σ_add_pd = 0.2,
)

sims = simobs(inf_2cmt_lin_turnover, subject, turnover_params; obstimes = sd_obstimes)
SimulatedObservations
  Simulated variables: dv, resp
  Time: [0.25, 0.5, 0.75, 1.0, 2.0, 4.0, 8.0, 12.0, 16.0, 20.0, 24.0, 36.0, 48.0, 60.0, 71.9]

We can then plot the simulated observations by using the sim_plot function:

sim_plot(inf_2cmt_lin_turnover, sims; observations = :dv)
Example block output

Note that we can further modify the plot, see Customizing Plots. For example,

sim_plot(
    inf_2cmt_lin_turnover,
    sims,
    observations = :dv,
    figure = (; fontsize = 18),
    axis = (; xlabel = "Time (hr)", ylabel = "Concentration (ng/mL)"),
)
Example block output

When we only give the structural model parameters, the random effects are automatically sampled from their distributions. If we wish to prescribe a value for the random effects, we pass initial values similarly:

vrandeffs = [(; ηpk = rand(4), ηpd = rand()) for _ in pop]
sims =
    simobs(inf_2cmt_lin_turnover, pop, turnover_params, vrandeffs; obstimes = sd_obstimes)
Simulated population (Vector{<:Subject})
  Simulated subjects: 10
  Simulated variables: dv, resp

If a population simulation is required with the inter-subject random effects set to zero, then it is possible to use the zero_randeffs function:

etas = zero_randeffs(inf_2cmt_lin_turnover, turnover_params, pop)

sims = simobs(inf_2cmt_lin_turnover, pop, turnover_params, etas; obstimes = sd_obstimes)

fig = Figure(; resolution = (1000, 600), fontsize = 18)
ax1 = Axis(fig[1, 1]; xlabel = "Time (hr)", ylabel = "Concentration (ng/mL)")
ax2 = Axis(fig[1, 2]; xlabel = "Time (hr)", ylabel = "Response")
sim_plot!(ax1, inf_2cmt_lin_turnover, sims; observations = [:dv])

sim_plot!(ax2, inf_2cmt_lin_turnover, sims; observations = [:resp])

fig
Example block output

We still see variability in the plot above due to the residual variability components in the model. If we set the residual error variance parameters to tiny values as below

turnover_params_wo_sigma =
    (turnover_params..., σ_prop_pk = 1e-20, σ_add_pk = 1e-20, σ_add_pd = 1e-20)
(tvcl = 1.5,
 tvvc = 25.0,
 tvq = 5.0,
 tvvp = 150.0,
 tvturn = 10,
 tvebase = 10,
 tvec50 = 0.3,
 Ω_pk = [0.05 0.0 0.0 0.0; 0.0 0.05 0.0 0.0; 0.0 0.0 0.05 0.0; 0.0 0.0 0.0 0.05],
 Ω_pd = [0.05;;],
 σ_prop_pk = 1.0e-20,
 σ_add_pk = 1.0e-20,
 σ_add_pd = 1.0e-20,)

and perform the simulation without residual random effects we get only the population mean:

etas = zero_randeffs(inf_2cmt_lin_turnover, turnover_params, pop)

sims = simobs(
    inf_2cmt_lin_turnover,
    pop,
    turnover_params_wo_sigma,
    etas;
    obstimes = sd_obstimes,
)

fig = Figure(; resolution = (1000, 600), fontsize = 18)
ax1 = Axis(fig[1, 1]; xlabel = "Time (hr)", ylabel = "Concentration (ng/mL)")
ax2 = Axis(fig[1, 2]; xlabel = "Time (hr)", ylabel = "Response")

sim_plot!(ax1, inf_2cmt_lin_turnover, sims; observations = [:dv])

sim_plot!(ax2, inf_2cmt_lin_turnover, sims; observations = [:resp])

fig
Example block output

Handling the resulting simulation data

The resulting sims object contains the results of the simulation for each subject. The results for the ith subject can be retrieved by sims[i].

We can convert the simulation results to Pumas Subjects using the constructor broadcasted over all elements of sims:

sims_subjects = Subject.(sims)
Population
  Subjects: 10
  Observations: dv, resp

The newly created Subjects will contain the variables defined in the @derived-block as their observations. Then sims_subjects can be used to fit the model parameters based on the simulated dataset. We can also turn the simulated population into a DataFrame with the DataFrame constructor:

DataFrame(sims[1])
16×30 DataFrame
RowidtimedvrespevidamtcmtratedurationssiirouteCentralPeripheralRespCLVcQVpebaseec50emaxturnkoutkin0ηpk_1ηpk_2ηpk_3ηpk_4ηpd_1
StringFloat64Float64?Float64?Int64Float64?Symbol?Float64?Float64?Int8?Float64?NCA.Route?Float64?Float64?Float64?Float64?Float64?Float64?Float64?Float64?Float64?Int64?Int64?Float64?Float64?Float64Float64Float64Float64Float64
110.0missingmissing1150.0Central10.015.000.0NullRoutemissingmissingmissing1.525.05.0150.010.00.31100.11.00.00.00.00.00.0
210.250.0968269.966630missingmissingmissingmissingmissingmissingmissing2.420650.06099929.966631.525.05.0150.010.00.31100.11.00.00.00.00.00.0
310.50.1875979.888390missingmissingmissingmissingmissingmissingmissing4.689930.2382039.888391.525.05.0150.010.00.31100.11.00.00.00.00.00.0
410.750.2727319.784090missingmissingmissingmissingmissingmissingmissing6.818280.5233749.784091.525.05.0150.010.00.31100.11.00.00.00.00.00.0
511.00.3526169.663470missingmissingmissingmissingmissingmissingmissing8.81540.9088439.663471.525.05.0150.010.00.31100.11.00.00.00.00.00.0
612.00.6265319.105460missingmissingmissingmissingmissingmissingmissing15.66333.318139.105461.525.05.0150.010.00.31100.11.00.00.00.00.00.0
714.01.011527.936490missingmissingmissingmissingmissingmissingmissing25.288111.19167.936491.525.05.0150.010.00.31100.11.00.00.00.00.00.0
818.01.430675.957480missingmissingmissingmissingmissingmissingmissing35.766833.23055.957481.525.05.0150.010.00.31100.11.00.00.00.00.00.0
9112.01.658694.525620missingmissingmissingmissingmissingmissingmissing41.467358.2114.525621.525.05.0150.010.00.31100.11.00.00.00.00.00.0
10116.01.470423.528070missingmissingmissingmissingmissingmissingmissing36.760482.72873.528071.525.05.0150.010.00.31100.11.00.00.00.00.00.0
11120.00.817013.105220missingmissingmissingmissingmissingmissingmissing20.425292.56123.105221.525.05.0150.010.00.31100.11.00.00.00.00.00.0
12124.00.5993893.092280missingmissingmissingmissingmissingmissingmissing14.984793.86673.092281.525.05.0150.010.00.31100.11.00.00.00.00.00.0
13136.00.4670293.557060missingmissingmissingmissingmissingmissingmissing11.675787.98343.557061.525.05.0150.010.00.31100.11.00.00.00.00.00.0
14148.00.4267723.897150missingmissingmissingmissingmissingmissingmissing10.669380.96333.897151.525.05.0150.010.00.31100.11.00.00.00.00.00.0
15160.00.3923674.143960missingmissingmissingmissingmissingmissingmissing9.8091674.45594.143961.525.05.0150.010.00.31100.11.00.00.00.00.00.0
16171.90.3610734.360540missingmissingmissingmissingmissingmissingmissing9.0268368.51774.360541.525.05.0150.010.00.31100.11.00.00.00.00.00.0

From there, any Julia tools can be used to analyze these arrays and DataFrames. Or you can conveniently export this to a CSV, Excel, or SAS data file.

Conducting an Estimation

With a model and population you can fit your model to your data. Model fitting in Pumas has the purpose of estimating parameters of a model, given a dataset, and is performed by calling the fit function with the following positional arguments:

  1. a PumasModel constructed via @model, or a PumasEMModel constructed via @emmodel.
  2. a Population or a Subject.
  3. a named tuple of the initial parameter values.
  4. an inference algorithm.

If you want to use the model's initial parameter values declared inside the @param block, you can do so with init_params(model). Note that this will fall back to the parameters' default values if you do not specify an initial value inside the @param block.

Pumas supports the following inference algorithms:

  • Maximum Likelihood Estimation:

    • NaivePooled: first order approximation without inter-subject random-effects.
    • FO: first-order approximation.
    • FOCE: first-order conditional estimation.
    • LaplaceI: second-order Laplace approximation.
    • SAEM: Stochastic Approximation using the Expectation-Maximization algorithm. Note that for this you use @emmodel instead of @model.
  • Bayesian inference:

    • BayesMCMC: MCMC using No-U-Turn Sampler (NUTS).
    • MarginalMCMC: MCMC using No-U-Turn Sampler (NUTS), but marginalizing out the subject-specific parameters.
    • MAP: Maximum A Posteriori estimation.

Now we can simultaneously fit our PK & PD models using one of the simulated populations from the previous section:

turnover_params = (;
    tvcl = 1.5,
    tvvc = 25.0,
    tvq = 5.0,
    tvvp = 150.0,
    tvturn = 10,
    tvebase = 10,
    tvec50 = 0.3,
    Ω_pk = Diagonal([0.05, 0.05, 0.05, 0.05]),
    Ω_pd = Diagonal([0.05]),
    σ_prop_pk = 0.15,
    σ_add_pk = 0.5,
    σ_add_pd = 0.2,
)

sims = simobs(inf_2cmt_lin_turnover, pop, turnover_params; obstimes = sd_obstimes)
Simulated population (Vector{<:Subject})
  Simulated subjects: 10
  Simulated variables: dv, resp
fit_2cmt_lin_turnover =
    fit(inf_2cmt_lin_turnover, Subject.(sims), init_params(inf_2cmt_lin_turnover), FOCE())
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:               Nonlinear ODE
Solver(s):                         (Vern7,Rodas5P)

Log-likelihood value:                   -156.74137
Number of subjects:                             10
Number of parameters:         Fixed      Optimized
                                  0             15
Observation records:         Active        Missing
    dv:                         150              0
    resp:                       150              0
    Total:                      300              0

--------------------------
              Estimate
--------------------------
tvcl           1.1682
tvvc          26.801
tvq            6.3267
tvvp         165.69
Ω_pk₁,₁        0.17137
Ω_pk₂,₂        0.091537
Ω_pk₃,₃        0.055024
Ω_pk₄,₄        2.9089e-74
σ_prop_pk      0.00084244
σ_add_pk       0.53685
tvturn         9.8483
tvebase       11.016
tvec50         0.28956
Ω_pd₁,₁        0.025593
σ_add_pd       0.21128
--------------------------

Calculating confidence intervals and Bootstrap

You can calculate confidence intervals using the infer function (by default 95% CIs):

infer_2cmt_lin_turnover = infer(fit_2cmt_lin_turnover)
Asymptotic inference results

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:               Nonlinear ODE
Solver(s):                         (Vern7,Rodas5P)

Log-likelihood value:                   -156.74137
Number of subjects:                             10
Number of parameters:         Fixed      Optimized
                                  0             15
Observation records:         Active        Missing
    dv:                         150              0
    resp:                       150              0
    Total:                      300              0

------------------------------------------------
              Estimate       SE      95.0% C.I.
------------------------------------------------
tvcl           1.1682      NaN     [NaN; NaN]
tvvc          26.801       NaN     [NaN; NaN]
tvq            6.3267      NaN     [NaN; NaN]
tvvp         165.69        NaN     [NaN; NaN]
Ω_pk₁,₁        0.17137     NaN     [NaN; NaN]
Ω_pk₂,₂        0.091537    NaN     [NaN; NaN]
Ω_pk₃,₃        0.055024    NaN     [NaN; NaN]
Ω_pk₄,₄        2.9089e-74  NaN     [NaN; NaN]
σ_prop_pk      0.00084244  NaN     [NaN; NaN]
σ_add_pk       0.53685     NaN     [NaN; NaN]
tvturn         9.8483      NaN     [NaN; NaN]
tvebase       11.016       NaN     [NaN; NaN]
tvec50         0.28956     NaN     [NaN; NaN]
Ω_pd₁,₁        0.025593    NaN     [NaN; NaN]
σ_add_pd       0.21128     NaN     [NaN; NaN]
------------------------------------------------


Variance-covariance matrix could not be
evaluated. The random effects may be over-
parameterized. Check the coefficients for
variance estimates near zero.

You can also convert the output to a DataFrame with the coeftable function:

coeftable(infer_2cmt_lin_turnover)
15×6 DataFrame
Rowparameterestimateserelative_seci_lowerci_upper
StringFloat64Float64Float64Float64Float64
1tvcl1.16817NaNNaNNaNNaN
2tvvc26.8008NaNNaNNaNNaN
3tvq6.32672NaNNaNNaNNaN
4tvvp165.695NaNNaNNaNNaN
5Ω_pk₁,₁0.171366NaNNaNNaNNaN
6Ω_pk₂,₂0.0915374NaNNaNNaNNaN
7Ω_pk₃,₃0.0550237NaNNaNNaNNaN
8Ω_pk₄,₄2.90889e-74NaNNaNNaNNaN
9σ_prop_pk0.000842441NaNNaNNaNNaN
10σ_add_pk0.536852NaNNaNNaNNaN
11tvturn9.84831NaNNaNNaNNaN
12tvebase11.0164NaNNaNNaNNaN
13tvec500.289555NaNNaNNaNNaN
14Ω_pd₁,₁0.0255929NaNNaNNaNNaN
15σ_add_pd0.211277NaNNaNNaNNaN

The returned table also includes relative standard errors (RSEs). For inference with bootstrap-based confidence intervals you can use Bootstrap as an optional second positional argument in infer:

Pumas.BootstrapType
Bootstrap(; rng=Random.default_rng, samples=200, stratify_by=nothing, ensemblealg=EnsembleThreads())

Constructs a Bootstrap object that allows the user to use bootstrapping with the supplies setting in infer function calls. See ?infer and https://docs.pumas.ai/stable/analysis/inference/ for more information about obtaining inference results.

The keyword arguments are as follows:

  • rng: used to select a pseudo-random number generator to beused for the random resamples of the datasets.
  • samples: controls the number of resampled Populations.
  • resamples: the number of model parameter samples from the original sample distribution to re-sample (by weighed re-sampling without replacement). The number of resamples must be less than that of the original sample. An original sample-to-resample ratio of at least 5:1 is suggested.
  • stratify_by control the stratification used in the re-sampling scheme. Should be a set to a Symbol with the name of a covariate it is possible to stratify by a covariate with a finite number of possible values (e.g :study for a covariate named study).
  • ensemblealg controls the parallelization across fit calls. Each fit is fitted using the ensemblealg that was used in the original fit.

Performing the bootstrap for our model can be performed as:

infer_boot_2cmt_lin_turnover =
    infer(fit_2cmt_lin_turnover, Bootstrap(; ensemblealg = EnsembleThreads()))
Bootstrap inference results

Successful minimization:                      true

Likelihood approximation:                     FOCE
Likelihood Optimizer:                         BFGS
Dynamical system type:               Nonlinear ODE
Solver(s):                         (Vern7,Rodas5P)

Log-likelihood value:                   -156.74137
Number of subjects:                             10
Number of parameters:         Fixed      Optimized
                                  0             15
Observation records:         Active        Missing
    dv:                         150              0
    resp:                       150              0
    Total:                      300              0

--------------------------------------------------------------------------------
              Estimate             SE                         95.0% C.I.
--------------------------------------------------------------------------------
tvcl           1.1682            0.26444           [  0.53796   ;   1.6056    ]
tvvc          26.801             6.5498            [ 19.663     ;  44.861     ]
tvq            6.3267            0.63479           [  5.5201    ;   7.9791    ]
tvvp         165.69             27.436             [119.75      ; 226.47      ]
Ω_pk₁,₁        0.17137           0.19149           [  0.051904  ;   0.64125   ]
Ω_pk₂,₂        0.091537          0.070162          [  4.0211e-9 ;   0.22477   ]
Ω_pk₃,₃        0.055024          0.033284          [  5.8373e-9 ;   0.10965   ]
Ω_pk₄,₄        2.9089e-74        3.9371e-90        [  2.9089e-74;   2.9089e-74]
σ_prop_pk      0.00084244        4.2317e-7         [  0.00084171;   0.00084318]
σ_add_pk       0.53685           0.025487          [  0.48387   ;   0.58449   ]
tvturn         9.8483            0.51322           [  8.4322    ;  10.511     ]
tvebase       11.016             0.5216            [ 10.117     ;  12.141     ]
tvec50         0.28956           0.030112          [  0.24576   ;   0.3583    ]
Ω_pd₁,₁        0.025593          0.0080927         [  0.0067727 ;   0.039601  ]
σ_add_pd       0.21128           0.015438          [  0.17598   ;   0.23344   ]
--------------------------------------------------------------------------------
Successful fits: 200 out of 200
Unique resampled populations: 200
No stratification.

You can also convert the Bootstrap inference results to a DataFrame with coeftable:

coeftable(infer_boot_2cmt_lin_turnover)
15×6 DataFrame
Rowparameterestimateserelative_seci_lowerci_upper
StringFloat64Float64Float64Float64Float64
1tvcl1.168170.2644380.226370.5379571.60558
2tvvc26.80086.549830.2443919.662744.8611
3tvq6.326720.6347940.1003355.520087.97907
4tvvp165.69527.43590.165581119.754226.474
5Ω_pk₁,₁0.1713660.1914921.117450.05190370.641253
6Ω_pk₂,₂0.09153740.07016180.7664834.02113e-90.224768
7Ω_pk₃,₃0.05502370.03328380.60495.83728e-90.109651
8Ω_pk₄,₄2.90889e-743.93713e-901.35348e-162.90889e-742.90889e-74
9σ_prop_pk0.0008424414.23165e-70.0005023080.0008417150.000843182
10σ_add_pk0.5368520.02548670.04747440.4838650.584492
11tvturn9.848310.5132190.05211248.4322110.5107
12tvebase11.01640.5216020.047347910.117112.1409
13tvec500.2895550.03011230.1039950.2457570.358304
14Ω_pd₁,₁0.02559290.008092690.3162090.006772690.0396012
15σ_add_pd0.2112770.01543760.0730680.1759780.233437

Conclusion

This tutorial covered the basic workflow of constructing a model, simulating observations from it, and conducting parameter estimation. You can also see subsequent tutorials that discuss the components in more detail, such as:

  1. More detailed treatment of specifying populations, dosage regimens, and covariates.
  2. Reading in dosage regimens and observations from standard input data.
  3. Fitting models with different estimation methods.