Introduction to Simulations in Pumas

Author

Haden Bunn

using Pumas, PharmaDatasets
using DataFramesMeta, CategoricalArrays
using CairoMakie, AlgebraOfGraphics
using SummaryTables

1 Introduction

Simulation is a fundamental part of most pharmacometrics analyses. Fortunately, Pumas provides a powerful interface for simulating a variety of endpoints without the need for additional software. This tutorial will serve as a basic introduction to simulation in Pumas and will focus on simulating additional observations or alternative dosing scenarios in an existing population. Advanced simulation (e.g., clinical trial simulation) will be discussed in later tutorials.

2 Getting Started

Simulations are performed using the simobs function which has multiple methods, all of which are detailed in the Pumas documentation. Checking ?simobs in the REPL provides the following function signature:

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),
)

The first three positional arguments (model, population, param) are all required.

model expects an AbstractPumasModel, which (for now) refers to a model defined using the @model macro.
population accepts a Population which was discussed in Module 4), or a single Subject which will be discussed later in the current module.
param should be a single parameter set defined as a NamedTuple or a vector of such parameter sets.

The remaining arguments have default values and need not be defined explicitly; however, it is worth knowing how the defaults affect each simulation.

randeffs is used to specify the random effect (“eta”) values for each subject. If left to the default, these values are generated by sampling the associated prior distributions defined in the model.
ensemblealg is used to select the parallelization mode to be used for the simulation.
diffeq_options can be used to pass additional options to the differential equation solver if the model does not have an analytical solution.
rng can be used to specify the random number generator to be used for the simulation.
simulate_error can be used to disable (false) the inclusion of RUV in the value returned by the predictive distribution’s error model.

Reproducibility

Many users are likely familiar with the concept of a random number generator (RNG) and the role they play in computational exercises where values are randomly sampled from a distribution. Using an RNG will make it (nearly) impossible to reproduce the results of a simulation unless steps are taken at the start to ensure reproducibility. In short, your results will differ slightly from those in the tutorial if you are executing the code locally, and that is to be expected. We will discuss this topic in greater detail later; for now, just focus on understanding the simulation workflow.

3 Setup

The examples below were created using the final integrated PK/PD model for warfarin. If you have downloaded this tutorial and are working through it locally, make sure you execute the code in the setup block before continuing. The example code assumes that the warfarin dataset (adeval), model (mdl), and fitted pumas model (“fpm”, myfit) all exist in your current session.

Warfarin Model

#* read dataset from PharmaDatasets
adeval = dataset("pumas/warfarin_pumas")

# population
mypop = read_pumas(adeval; observations = [:conc, :pca], covariates = [:wtbl])

# warfarin model
mdl = @model begin

    @metadata begin
        desc = "Integrated Warfarin PK/PD model"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        # Clearance, L/hr
        tvcl ∈ RealDomain(lower = 0.0, init = 0.134)
        # Volume of distribution, central, L
        tvvc ∈ RealDomain(lower = 0.0, init = 8.11)
        # absorption rate constant, hr^-1
        tvka ∈ RealDomain(lower = 0.0, init = 1.32)
        # absorption lag, hr
        tvalag ∈ RealDomain(lower = 0.0, init = 0.1)

        # PD parameters
        # Baseline, %
        tve0 ∈ RealDomain(lower = 0.0, init = 95, upper = 100)
        # Imax, %
        tvimax ∈ RealDomain(lower = 0.0, init = 0.8, upper = 1)
        # IC50, mg/L
        tvic50 ∈ RealDomain(lower = 0.0, init = 1.0)
        # Turnover
        tvturn ∈ RealDomain(lower = 0.0, init = 14.0)
        # Inter-individual variability
        """
          - Ωcl
          - Ωvc
          - Ωka
        """
        Ωpk ∈ PDiagDomain([0.01, 0.01, 0.01])
        """
          - Ωe0
          - Ωic50
        """
        Ωpd ∈ PDiagDomain([0.01, 0.01])
        # Residual variability
        # proportional error, pk
        σprop_pk ∈ RealDomain(lower = 0.0, init = 0.2)
        # additive error, pk, mg/L
        σadd_pk ∈ RealDomain(lower = 0.0, init = 0.2)
        # additive error, pca, %
        σadd_pd ∈ RealDomain(lower = 0.0, init = 1)
    end

    @random begin
        # mean = 0, covariance = Ωpk
        ηpk ~ MvNormal(Ωpk)
        # mean = 0, covariance = Ωpd
        ηpd ~ MvNormal(Ωpd)
    end

    @covariates wtbl

    @pre begin
        # PK
        cl = tvcl * (wtbl / 70)^0.75 * exp(ηpk[1])
        vc = tvvc * (wtbl / 70) * exp(ηpk[2])
        ka = tvka * exp(ηpk[3])
        # PD
        e0 = tve0 * exp(ηpd[1])
        imax = tvimax
        ic50 = tvic50 * exp(ηpd[2])
        turn = tvturn
        #kout = log(2) / turn
        kout = 1 / turn
        kin = e0 * kout
    end

    @dosecontrol begin
        lags = (depot = tvalag,)
    end

    @init begin
        e = e0
    end

    @vars begin
        # inhibitory model
        imdl := 1 - (imax * (central / vc)) / (ic50 + (central / vc))
    end

    @dynamics begin
        depot' = -ka * depot
        central' = ka * depot - (cl / vc) * central
        e' = kin * imdl - kout * e
    end

    @derived begin
        cp := @. central / vc
        # warfarin concentration, mg/L
        conc ~ @. CombinedNormal(cp, σadd_pk, σprop_pk)
        # prothrombin complex activity, % of normal
        pca ~ @. Normal(e, σadd_pd)
    end
end

# fitted pumas model, fpm
myfit = fit(mdl, mypop, init_params(mdl), FOCE())

[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter     Function value   Gradient norm 
     0     5.536943e+03     7.052675e+03
 * time: 0.041929006576538086
     1     1.762032e+03     1.253034e+03
 * time: 4.214653015136719
     2     1.449413e+03     6.274412e+02
 * time: 5.1580259799957275
     3     1.331405e+03     2.099880e+02
 * time: 6.136675834655762
     4     1.310843e+03     1.627399e+02
 * time: 7.11238694190979
     5     1.295073e+03     1.376531e+02
 * time: 8.13167405128479
     6     1.282335e+03     1.324314e+02
 * time: 9.108187913894653
     7     1.275524e+03     9.521788e+01
 * time: 10.01130485534668
     8     1.271552e+03     5.354883e+01
 * time: 10.979696035385132
     9     1.269315e+03     2.407046e+01
 * time: 11.95639705657959
    10     1.268615e+03     2.368964e+01
 * time: 12.954828977584839
    11     1.268257e+03     2.366048e+01
 * time: 13.931599855422974
    12     1.267352e+03     3.190950e+01
 * time: 14.926499843597412
    13     1.265336e+03     5.399188e+01
 * time: 15.84483790397644
    14     1.260582e+03     8.775888e+01
 * time: 16.696285009384155
    15     1.252929e+03     1.110735e+02
 * time: 17.467474937438965
    16     1.246526e+03     8.328859e+01
 * time: 18.174316883087158
    17     1.240724e+03     8.711807e+01
 * time: 18.864006996154785
    18     1.236881e+03     7.439216e+01
 * time: 19.602386951446533
    19     1.232291e+03     6.878129e+01
 * time: 20.304532051086426
    20     1.227181e+03     5.604034e+01
 * time: 20.992873907089233
    21     1.223093e+03     3.856728e+01
 * time: 21.642499923706055
    22     1.221064e+03     5.186424e+01
 * time: 22.241209030151367
    23     1.220113e+03     4.306874e+01
 * time: 22.866899967193604
    24     1.218936e+03     2.178258e+01
 * time: 23.468889951705933
    25     1.217831e+03     2.144752e+01
 * time: 24.06534504890442
    26     1.216775e+03     2.201560e+01
 * time: 24.668293952941895
    27     1.215938e+03     2.969273e+01
 * time: 25.258220911026
    28     1.214922e+03     3.438177e+01
 * time: 25.86903405189514
    29     1.212881e+03     4.013028e+01
 * time: 26.46351385116577
    30     1.207231e+03     5.214921e+01
 * time: 27.079733848571777
    31     1.197000e+03     6.567320e+01
 * time: 27.689594984054565
    32     1.189039e+03     6.178092e+01
 * time: 28.297005891799927
    33     1.173765e+03     7.561375e+01
 * time: 28.902008056640625
    34     1.171088e+03     2.385028e+01
 * time: 29.517082929611206
    35     1.170366e+03     2.359334e+01
 * time: 30.124361991882324
    36     1.170030e+03     2.394220e+01
 * time: 30.718626022338867
    37     1.169536e+03     2.287562e+01
 * time: 31.308125972747803
    38     1.168453e+03     1.921641e+01
 * time: 31.906085968017578
    39     1.166697e+03     1.661199e+01
 * time: 32.52202796936035
    40     1.164753e+03     2.146306e+01
 * time: 33.16405701637268
    41     1.163396e+03     1.589737e+01
 * time: 33.8055579662323
    42     1.162810e+03     1.879805e+01
 * time: 34.45263195037842
    43     1.162593e+03     1.922942e+01
 * time: 35.090128898620605
    44     1.162146e+03     1.865903e+01
 * time: 35.74095106124878
    45     1.161045e+03     2.060103e+01
 * time: 36.3992018699646
    46     1.158095e+03     4.011929e+01
 * time: 37.02961206436157
    47     1.150009e+03     6.883383e+01
 * time: 37.66959095001221
    48     1.147406e+03     7.331888e+01
 * time: 38.44691705703735
    49     1.144610e+03     7.504905e+01
 * time: 39.15776491165161
    50     1.138023e+03     6.294975e+01
 * time: 39.871328830718994
    51     1.132705e+03     4.283014e+01
 * time: 40.554584980010986
    52     1.129075e+03     2.327948e+01
 * time: 41.29772686958313
    53     1.128255e+03     2.383372e+01
 * time: 42.03863286972046
    54     1.127697e+03     2.546344e+01
 * time: 42.74056005477905
    55     1.126978e+03     2.642229e+01
 * time: 43.47938084602356
    56     1.126254e+03     2.679647e+01
 * time: 44.22636103630066
    57     1.124269e+03     2.481478e+01
 * time: 44.94812798500061
    58     1.121409e+03     3.925976e+01
 * time: 45.65275192260742
    59     1.118180e+03     3.554727e+01
 * time: 46.338696002960205
    60     1.116311e+03     1.691328e+01
 * time: 47.0559139251709
    61     1.115933e+03     1.725048e+01
 * time: 47.75369906425476
    62     1.115750e+03     1.702130e+01
 * time: 48.462493896484375
    63     1.115157e+03     1.748010e+01
 * time: 49.17825388908386
    64     1.113379e+03     2.744225e+01
 * time: 49.93651485443115
    65     1.107471e+03     4.999759e+01
 * time: 50.821863889694214
    66     1.099217e+03     7.740446e+01
 * time: 52.21660399436951
    67     1.098357e+03     6.707598e+01
 * time: 53.477942943573
    68     1.093578e+03     3.715572e+01
 * time: 54.570740938186646
    69     1.090146e+03     1.668128e+01
 * time: 55.70665097236633
    70     1.088869e+03     1.082857e+01
 * time: 56.87196397781372
    71     1.087895e+03     9.101157e+00
 * time: 57.997851848602295
    72     1.086966e+03     1.225553e+01
 * time: 59.092018842697144
    73     1.086349e+03     1.200735e+01
 * time: 60.06861400604248
    74     1.086021e+03     9.320698e+00
 * time: 60.98715901374817
    75     1.085858e+03     8.113654e+00
 * time: 61.90627098083496
    76     1.085776e+03     8.127231e+00
 * time: 62.82839798927307
    77     1.085607e+03     7.969663e+00
 * time: 63.7496018409729
    78     1.085219e+03     7.491324e+00
 * time: 64.68053698539734
    79     1.084323e+03     1.248809e+01
 * time: 65.6185450553894
    80     1.082506e+03     2.097591e+01
 * time: 66.59559392929077
    81     1.079392e+03     2.578708e+01
 * time: 67.59074687957764
    82     1.077539e+03     8.916359e+00
 * time: 68.5765769481659
    83     1.077366e+03     1.950845e+00
 * time: 69.53595900535583
    84     1.077341e+03     1.114594e+00
 * time: 70.57854795455933
    85     1.077339e+03     1.131454e+00
 * time: 71.60485696792603
    86     1.077338e+03     1.121807e+00
 * time: 72.64227986335754
    87     1.077335e+03     1.078000e+00
 * time: 73.69576406478882
    88     1.077329e+03     1.045586e+00
 * time: 74.67245483398438
    89     1.077316e+03     1.384979e+00
 * time: 75.63800287246704
    90     1.077292e+03     1.537679e+00
 * time: 76.60296487808228
    91     1.077257e+03     1.232746e+00
 * time: 77.59374785423279
    92     1.077228e+03     5.729795e-01
 * time: 78.58315396308899
    93     1.077219e+03     5.059871e-01
 * time: 79.57540392875671
    94     1.077218e+03     5.033585e-01
 * time: 80.56820106506348
    95     1.077218e+03     5.028759e-01
 * time: 81.52740001678467
    96     1.077218e+03     5.022155e-01
 * time: 82.4892029762268
    97     1.077217e+03     5.008367e-01
 * time: 83.4532060623169
    98     1.077216e+03     4.977185e-01
 * time: 84.43835496902466
    99     1.077212e+03     4.905506e-01
 * time: 85.4375410079956
   100     1.077202e+03     7.605415e-01
 * time: 86.4151120185852
   101     1.077182e+03     1.030054e+00
 * time: 87.38476085662842
   102     1.077153e+03     1.038802e+00
 * time: 88.35336494445801
   103     1.077130e+03     6.003895e-01
 * time: 89.362233877182
   104     1.077123e+03     5.073914e-01
 * time: 90.80477499961853
   105     1.077122e+03     4.817840e-01
 * time: 92.52817797660828
   106     1.077122e+03     4.675479e-01
 * time: 94.28085494041443
   107     1.077122e+03     4.482687e-01
 * time: 95.63090205192566
   108     1.077121e+03     4.167878e-01
 * time: 96.72366285324097
   109     1.077119e+03     4.159409e-01
 * time: 97.74860095977783
   110     1.077115e+03     7.458466e-01
 * time: 98.830894947052
   111     1.077106e+03     1.228748e+00
 * time: 99.86471700668335
   112     1.077083e+03     1.875833e+00
 * time: 100.90498304367065
   113     1.077043e+03     2.711460e+00
 * time: 101.96883392333984
   114     1.076990e+03     3.472380e+00
 * time: 103.03019595146179
   115     1.076942e+03     2.744371e+00
 * time: 104.08783888816833
   116     1.076925e+03     9.131372e-01
 * time: 105.15365886688232
   117     1.076924e+03     9.225451e-01
 * time: 106.2352979183197
   118     1.076924e+03     9.211518e-01
 * time: 107.31336498260498
   119     1.076923e+03     9.099676e-01
 * time: 108.40698194503784
   120     1.076920e+03     1.072485e+00
 * time: 109.48189783096313
   121     1.076913e+03     1.902001e+00
 * time: 110.57382488250732
   122     1.076897e+03     2.866993e+00
 * time: 111.64716506004333
   123     1.076863e+03     3.341730e+00
 * time: 112.72743487358093
   124     1.076807e+03     2.150254e+00
 * time: 113.79986596107483
   125     1.076761e+03     2.997759e-01
 * time: 114.87197303771973
   126     1.076749e+03     7.551212e-01
 * time: 115.96878790855408
   127     1.076747e+03     8.967956e-01
 * time: 117.07808899879456
   128     1.076746e+03     7.027883e-01
 * time: 118.14830684661865
   129     1.076744e+03     5.081113e-01
 * time: 119.22464895248413
   130     1.076743e+03     1.635722e-01
 * time: 120.33754396438599
   131     1.076741e+03     1.794520e-01
 * time: 121.48296093940735
   132     1.076740e+03     1.890041e-01
 * time: 122.6086778640747
   133     1.076740e+03     1.923763e-01
 * time: 123.71921491622925
   134     1.076740e+03     1.972343e-01
 * time: 124.84014391899109
   135     1.076740e+03     2.001877e-01
 * time: 125.95519185066223
   136     1.076740e+03     2.031744e-01
 * time: 127.11352896690369
   137     1.076740e+03     2.057411e-01
 * time: 128.3254270553589
   138     1.076740e+03     2.081860e-01
 * time: 129.48734283447266
   139     1.076740e+03     2.083823e-01
 * time: 130.57218194007874
   140     1.076740e+03     2.000831e-01
 * time: 131.65658283233643
   141     1.076739e+03     1.697498e-01
 * time: 132.7576768398285
   142     1.076738e+03     1.576128e-01
 * time: 133.86483192443848
   143     1.076737e+03     7.984908e-02
 * time: 134.97701406478882
   144     1.076737e+03     1.937515e-02
 * time: 136.08467388153076
   145     1.076737e+03     2.909238e-03
 * time: 137.18372201919556
   146     1.076737e+03     8.635316e-04
 * time: 138.26713705062866

FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             32

Observation records:         Active        Missing
    conc:                       251             47
    pca:                        232             66
    Total:                      483            113

Number of parameters:      Constant      Optimized
                                  0             16

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -1076.7369

----------------------
           Estimate
----------------------
tvcl        0.13551
tvvc        7.9849
tvka        1.1742
tvalag      0.87341
tve0       96.616
tvimax      1.0
tvic50      1.1737
tvturn     18.828
Ωpk₁,₁      0.069236
Ωpk₂,₂      0.021885
Ωpk₃,₃      0.83994
Ωpd₁,₁      0.0028192
Ωpd₂,₂      0.18162
σprop_pk    0.08854
σadd_pk     0.4189
σadd_pd     4.1536
----------------------

4 Simulation Basics

Simulation complexity increases as the conditions in a given scenario diverge from the set of conditions (e.g., population, dosage regimen) used to develop the underlying model.
We examine four scenarios with increasing complexity to introduce the user to simulations in Pumas.

4.1 Scenario 1

We begin with the simple goal of generating complete profiles for each subject in the original dataset.
- Simple because underlying population and dosage regimen are unchanged.
Two approaches, predict and simobs
- predict not technically simulation, but end result is the same, primary difference is the lack of RUV compared to simobs. The random effects will additionally be set to the empirical bayes estimates.

# sampling times from studies
stimes = [0.5, 1, 1.5, 2, 3, 6, 9, 12, 24, 36, 48, 72, 96, 120, 144]

# predictions based on individual EBEs and obstimes
# if passed, obstimes are merged with existing observation times vector
mypred = predict(myfit, myfit.data; obstimes = stimes)

Prediction
  Subjects: 32
  Predictions: conc, pca
  Covariates: wtbl

Application: NCA

With a few simple modifications, the mypred object can used for NCA. Refer to the documentation for more information on performing NCA in Pumas.

mynca = @chain DataFrame(mypred) begin
    # merge conc and conc_ipred into new column, conc_new
    transform([:conc, :conc_ipred] => ByRow((c, i) -> coalesce(c, i)) => :conc_new)
    # update route column for NCA
    @transform :route = "ev"
    # create nca population, specify obs col as conc_new
    read_nca(_; observations = :conc_new)
    run_nca(_)
end

NCA Report
     Timestamp: 2026-04-02T18:20:47.571
     Version number: 0.1.0

Output Parameters DataFrame
32×38 DataFrame
 Row │ id      dose     tlag     tmax     cmax      tlast    clast      clast_ ⋯
     │ String  Float64  Float64  Float64  Float64   Float64  Float64    Float6 ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ 1         100.0      0.5      9.0  10.8        144.0  0.124115    0.124 ⋯
   2 │ 2         100.0      0.5      9.0  11.2304     144.0  1.10587     1.093
   3 │ 3         120.0      0.5      9.0  14.4        144.0  2.25022     2.241
   4 │ 4          60.0      0.5      6.0  11.9        144.0  1.48714     1.505
   5 │ 5         113.0      0.5      3.0   8.93322    144.0  1.22051     1.200 ⋯
   6 │ 6          90.0      0.5      3.0  13.4        144.0  0.0272853   0.044
   7 │ 7         135.0      0.5      2.0  17.6        144.0  0.71204     0.873
   8 │ 8          75.0      0.5      9.0  12.9        144.0  0.808874    0.719
  ⋮  │   ⋮        ⋮        ⋮        ⋮        ⋮         ⋮         ⋮          ⋮  ⋱
  26 │ 27        120.0      0.5      6.0  15.3014     144.0  1.40167     1.548 ⋯
  27 │ 28        120.0      0.5      6.0  12.3473     144.0  1.58808     1.535
  28 │ 29        153.0      0.5      6.0  11.5899     144.0  1.24201     1.379
  29 │ 30        105.0      0.5      6.0  12.4077     144.0  1.35182     1.407
  30 │ 31        125.0      0.5      6.0  12.0114     144.0  1.59054     1.663 ⋯
  31 │ 32         93.0      0.5      6.0  11.3265     144.0  1.82322     1.816
  32 │ 33        100.0      0.5      6.0  11.638      144.0  1.75404     1.728
                                                  31 columns and 17 rows omitted

An example of a similar analysis using simobs.

simobs(myfit.model, myfit.data, coef(myfit), empirical_bayes(myfit); obstimes = stimes)

Simulated population (Vector{<:Subject})
  Simulated subjects: 32
  Simulated variables: conc, pca

Example post-processing for mypred

@chain DataFrame(mypred) begin
    # only observations
    filter(df -> df.evid == 0, _)
    # if conc missing, 1, else 0
    transform(:conc => ByRow(c -> ismissing(c) ? 1 : 0) => :isnew)
    # merge conc and conc_ipred into new column, conc_new
    transform([:conc, :conc_ipred] => ByRow((c, i) -> coalesce(c, i)) => :conc_new)
    # CT scatter plot
    data(_) *
    mapping(
        :time => "Time After Dose, hours",
        :conc_new => "Warfarin Concentration, mg/L",
        group = :id => nonnumeric,
        color = :isnew => renamer(0 => "No", 1 => "Yes") => "Predicted",
    ) *
    visual(Scatter)
    draw(
        _;
        axis = (;
            title = "Individual Concentration-Time Profiles",
            subtitle = "Warfarin ~ TAD",
        ),
    )
end

4.2 Scenario 2

In this scenario we generate observation (conc, pca) time profiles for the trial population following a one-time LD of 0.75 mg/kg PO.
- Same population, different dosage regimen lets us examine the Subject and DosageRegimen constructors without the additional complexity of creating a virtual population from scratch.

Key Concept: Subject Constructor

The Subject constructor is a fundamental part of most simulation workflows in Pumas. If you have not reviewed the corresponding tutorial, we recommend doing so before proceeding here.

Often, the best approach to building a simulation in Pumas is to focus on a single subject workflow, then, once everything is working, use a repeated-evaluation construct to complete the analysis.
In the initial setup, we showcase the mutating Subject syntax by accessing data from the first Subject stored in myfit.data[1].
- Converting the mg/kg dose to mg requires wtbl which we extract from the covariates field.

# first subject in population used in model fit
sub01 =
    Subject(myfit.data[1]; events = DosageRegimen(0.75 * myfit.data[1].covariates(0).wtbl))

sim01 = simobs(
    mdl,                        # model
    sub01,                      # subject or population of subjects                      
    coef(myfit),                # parameter estimates
    empirical_bayes(myfit)[1];  # random effects (i.e., EBEs)
    obstimes = stimes,          # obstimes for full study profile
    simulate_error = false,      # set RUV=0
)

SimulatedObservations
  Simulated variables: conc, pca
  Time: [0.5, 1.0, 1.5, 2.0, 3.0, 6.0, 9.0, 12.0, 24.0, 36.0, 48.0, 72.0, 96.0, 120.0, 144.0]

With a working single subject simulation, we can move on to simulating observations for a population.
Here, we take a slightly different approach by creating a dataframe of subject-level covariates and iterating over each row to create our population.
- Note, we could have iterated over all Subjects store in myfit.data and modified them as we did above; this syntax below shows an equivalent approach that may be more intuitive to new users.

# df with one row for each unique patient in original dataset
_patients = combine(groupby(adeval, :id), first)

# iterate over _patients creating 1 subject per row
pop02 = map(eachrow(_patients)) do r
    Subject(
        id = r.id[1],
        events = DosageRegimen(0.75 * r.wtbl[1]),
        covariates = (; wtbl = r.wtbl[1]),
    )
end

sim02 = simobs(
    mdl,
    pop02,
    coef(myfit),
    empirical_bayes(myfit);
    obstimes = stimes,
    simulate_error = false,
)

Simulated population (Vector{<:Subject})
  Simulated subjects: 32
  Simulated variables: conc, pca

A bit of additional post-processing
The figure below shows individual warfarin CT profiles for all subjects.

@chain DataFrame(sim02) begin
    data(_) *
    mapping(
        :time => "Time After Dose, hours",
        :conc => "Warfarin Concentration, mg/L",
        group = :id => nonnumeric,
    ) *
    visual(Scatter)
    draw(
        _;
        axis = (;
            title = "Individual Concentration-Time Profiles",
            subtitle = "0.75 mg/kg x1",
        ),
    )
end

Individual CT profiles can be created by stratifying the data using the layout kwarg in mapping, and then separated using the paginate function.
x|yticklabelsize was adjusted to improve readability along those axes.

# plot layers
_plt = @chain DataFrame(sim02) begin
    data(_) *
    mapping(
        :time => "Time After Dose, hours",
        :conc => "Warfarin Concentration, mg/L",
        group = :id => nonnumeric,
        layout = :id => nonnumeric,
    ) *
    visual(Scatter)
end

# draw(paginate(...)) returns a vector of `FigureGrid` objects
_pgrid = draw(
    paginate(_plt, layout = 16);
    figure = (;
        size = (6.5, 6.5) .* 96,
        title = "Individual Concentration-Time Profiles",
        subtitle = "0.75 mg/kg x1",
    ),
    axis = (; xticks = 0:24:144, xticklabelsize = 12, yticklabelsize = 12),
)

2-element Vector{AlgebraOfGraphics.FigureGrid}:
 FigureGrid()
 FigureGrid()

The result is a Vector{FigureGrid} with figures that can be accessed via indexing. The first 16 subjects are shown in the panel below.

_pgrid[1]

4.3 Scenario 3

In this scenario we assess the impact of augmented clearance on target attainment after a one-time LD of 1.5 mg/kg PO.
- We use this scenario to showcase creating a Subject from scratch along with the zero_randeffs helper function.

Introduction to Julia Callback

These last two scenarios should reinforce why Julia fundamentals are so important and why they were chosen for Module 1. We encourage the reader to revisit that tutorial if any of the code that follows is unclear.

# final parameter estimates
_params = coef(myfit)

sim03 = map([0.8, 1, 1.2]) do i
    simobs(
        mdl,
        Subject(
            id = "CL: $i",
            events = DosageRegimen(1.5 * 70),
            covariates = (; wtbl = 70),
        ),
        merge(_params, (; tvcl = _params.tvcl * i)),
        zero_randeffs(mdl, _params);
        obstimes = 0.5:0.5:144,
        simulate_error = false,
    )
end

Simulated population (Vector{<:Subject})
  Simulated subjects: 3
  Simulated variables: conc, pca

Visualize PK profile

@chain DataFrame(sim03) begin
    data(_) *
    mapping(
        :time => "Time After Dose, hours",
        :conc => "Warfarin Concentration, mg/L",
        group = :id => nonnumeric,
        color = :id => "Scenario",
    ) *
    visual(Lines)
    draw(
        _;
        axis = (;
            title = "Population Concentration-Time Profiles with Augmented CL",
            subtitle = "1.5 mg/kg x1",
        ),
    )
end

Visual PD profile
More complex figures in AoG can be easier to manage if their respective layers are stored in separate variables.

# band for therapeutic range
tr_layer = mapping(0:144, 20, 35) * visual(Band; color = (:blue, 0.2))

# profiles
profiles =
    data(DataFrame(sim03)) *
    mapping(
        :time => "Time After Dose, hours",
        :pca => "PCA, % of Normal",
        group = :id => nonnumeric,
        color = :id => "Scenario",
    ) *
    visual(Lines)

draw(
    tr_layer + profiles;
    axis = (;
        title = "Population Concentration-Time Profiles with Augmented CL",
        subtitle = "1.5 mg/kg x1",
    ),
)

4.4 Scenario 4

In this scenario we combine the concepts discussed above to evaluate three alternative dosage regimens: 5, 10, or 15 mg PO daily for 14 days.
- The estimated population half-life for warfarin per our model is ~41 hours which means it should take roughly 9 days (on average) to achieve steady-state; we extend this to 14 days to ensure each of our virtual subjects is at SS prior to evaluation.
We will generate a Population of 100 Subjects and use it simulate a total of 600 trials (200 per dosage regimen).
- We will not include RUV, since most variability comes from BSV and RUV can make results difficult to interpret.
We are interested in three metrics.
- The probability of obtaining a pca within the therapeutic range (20-35%) at any time during treatment.
- The time needed to reach the first therapeutic pca value.
- The total time spent in the TR as a percentage of the dosing interval (i.e., 24 hours) at SS (Day 14).

4.4.1 Setup

We begin, as before, by developing the code for a single simulation that we can then reuse for the remaining dosage regimens. While working through the setup, we will keep our code simple by limiting our “population” and replicates to 5. We will also focus on one dosage regimen, 5 mg PO daily for 14 days. This will allow us to spot check our code and the results to ensure the output it what we expect instead of trying to troubleshoot for the full population and profile. In order, we must:

Create a population of subjects that has a single covariate (wtbl) that is sampled from a uniform distribution of observed values (40-102 kg).
Simulate an appropriate number of observations (hourly observations will be sufficient).
Repeat the simulation in #2 for a total of 5 simulations.
Store the output in a format that will make post-processing and evaluation as easy as possible.
Process and evaluate the result then present our findings in a meaningful way.

4.4.1.1 Population

We can combine map with a Range between 1:n and a do-block to create a vector of virtual subjects (i.e., a Population)
For simplicity, we will also use scalar literal values for the range of wtbl in the Uniform call instead of obtaining them programmatically with extrema or some other function.

pop03 = map(1:5) do i
    _wtbl = rand(Uniform(40, 120))
    Subject(
        id = i,
        events = DosageRegimen(5; ii = 24, addl = 13),
        covariates = (; wtbl = _wtbl),
    )
end

Population
  Subjects: 5
  Covariates: wtbl
  Observations:

4.4.1.2 Simulation

Simulating observations for 2 subjects is comparatively simple, and we can repeat that simulation by mapping over a range 1:n as we did when creating pop03.
The resulting vector of SimulatedObservation objects can be concatenated into a single object using the reduce(vcat, myvectorofsims), then converted into a data frame for post-processing.
Since we need to summarize values from each simulation, we will need to include a variable to track the iteration number for each simulation. We can do this by leveraging the mutating Subject syntax to add a rep_id as a covariate for each subject pop03 inside the map before we call simobs.

goodsim = map(1:5) do i
    # rebuild pop03 using mutating Subject to add rep_id
    _pop = map(pop03) do s
        Subject(s; covariates = (; rep_id = i))
    end
    # simulate values for _pop 
    simobs(mdl, _pop, coef(myfit); obstimes = 0:1:(24*14), simulate_error = false)
end

5-element Vector{Vector{Pumas.SimulatedObservations{PumasModel{(tvcl = 1, tvvc = 1, tvka = 1, tvalag = 1, tve0 = 1, tvimax = 1, tvic50 = 1, tvturn = 1, Ωpk = 3, Ωpd = 2, σprop_pk = 1, σadd_pk = 1, σadd_pd = 1), 5, (:depot, :central, :e), (:cl, :vc, :ka, :e0, :imax, :ic50, :turn, :kout, :kin), ParamSet{@NamedTuple{tvcl::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvvc::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvka::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvalag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tve0::RealDomain{Float64, Int64, Int64}, tvimax::RealDomain{Float64, Int64, Float64}, tvic50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvturn::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, Ωpk::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, Ωpd::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σprop_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pd::RealDomain{Float64, TransformVariables.Infinity{true}, Int64}}}, ParamSet{@NamedTuple{ηpk::Pumas.RandomEffect{false, var"#1#12"}, ηpd::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{}, Pumas.ConstantInterpolationStructArray{Vector{Float64}, @NamedTuple{wtbl::Vector{Float64}, rep_id::Vector{Int64}}}, DosageRegimen{Float64, Symbol, Float64, Float64, Float64, Float64}, Nothing}, StepRange{Int64, Int64}, @NamedTuple{tvcl::Float64, tvvc::Float64, tvka::Float64, tvalag::Float64, tve0::Float64, tvimax::Float64, tvic50::Float64, tvturn::Float64, Ωpk::PDMats.PDiagMat{Float64, Vector{Float64}}, Ωpd::PDMats.PDiagMat{Float64, Vector{Float64}}, σprop_pk::Float64, σadd_pk::Float64, σadd_pd::Float64}, @NamedTuple{ηpk::Vector{Float64}, ηpd::Vector{Float64}}, @NamedTuple{wtbl::Vector{Float64}, rep_id::Vector{Int64}}, @NamedTuple{saveat::StepRange{Int64, Int64}, ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, SciMLBase.ODESolution{Float64, 2, Vector{StaticArraysCore.SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{StaticArraysCore.SVector{3, Float64}}}, Nothing, Any, Any, Any, SciMLBase.DEStats, Vector{Int64}, Nothing, Nothing, Nothing}, @NamedTuple{lags::Vector{@NamedTuple{depot::Float64}}}, @NamedTuple{cl::Vector{Float64}, vc::Vector{Float64}, ka::Vector{Float64}, e0::Vector{Float64}, imax::Vector{Float64}, ic50::Vector{Float64}, turn::Vector{Float64}, kout::Vector{Float64}, kin::Vector{Float64}}, @NamedTuple{conc::Vector{Float64}, pca::Vector{Float64}}, @NamedTuple{}}}}:
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca

We used the label goodsims because this approach, while valid, is redundant because we are recreating the Population from scratch with each simulation.
We could simplify the code by using a nested map call to create the population during each iteration.

bettersim = map(1:5) do i
    _pop = map(1:5) do s
        Subject(
            id = s,
            events = DosageRegimen(5, ii = 24, addl = 13),
            covariates = (; wtbl = rand(Uniform(40, 102)), rep_id = i),
        )
    end
    simobs(mdl, _pop, coef(myfit); obstimes = 0:1:(24*14), simulate_error = false)
end

5-element Vector{Vector{Pumas.SimulatedObservations{PumasModel{(tvcl = 1, tvvc = 1, tvka = 1, tvalag = 1, tve0 = 1, tvimax = 1, tvic50 = 1, tvturn = 1, Ωpk = 3, Ωpd = 2, σprop_pk = 1, σadd_pk = 1, σadd_pd = 1), 5, (:depot, :central, :e), (:cl, :vc, :ka, :e0, :imax, :ic50, :turn, :kout, :kin), ParamSet{@NamedTuple{tvcl::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvvc::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvka::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvalag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tve0::RealDomain{Float64, Int64, Int64}, tvimax::RealDomain{Float64, Int64, Float64}, tvic50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvturn::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, Ωpk::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, Ωpd::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σprop_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pd::RealDomain{Float64, TransformVariables.Infinity{true}, Int64}}}, ParamSet{@NamedTuple{ηpk::Pumas.RandomEffect{false, var"#1#12"}, ηpd::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{}, Pumas.ConstantCovar{@NamedTuple{wtbl::Float64, rep_id::Int64}}, DosageRegimen{Float64, Symbol, Float64, Float64, Float64, Float64}, Nothing}, StepRange{Int64, Int64}, @NamedTuple{tvcl::Float64, tvvc::Float64, tvka::Float64, tvalag::Float64, tve0::Float64, tvimax::Float64, tvic50::Float64, tvturn::Float64, Ωpk::PDMats.PDiagMat{Float64, Vector{Float64}}, Ωpd::PDMats.PDiagMat{Float64, Vector{Float64}}, σprop_pk::Float64, σadd_pk::Float64, σadd_pd::Float64}, @NamedTuple{ηpk::Vector{Float64}, ηpd::Vector{Float64}}, @NamedTuple{wtbl::Vector{Float64}, rep_id::Vector{Int64}}, @NamedTuple{saveat::StepRange{Int64, Int64}, ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, SciMLBase.ODESolution{Float64, 2, Vector{StaticArraysCore.SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{StaticArraysCore.SVector{3, Float64}}}, Nothing, Any, Any, Any, SciMLBase.DEStats, Vector{Int64}, Nothing, Nothing, Nothing}, @NamedTuple{lags::Vector{@NamedTuple{depot::Float64}}}, @NamedTuple{cl::Vector{Float64}, vc::Vector{Float64}, ka::Vector{Float64}, e0::Vector{Float64}, imax::Vector{Float64}, ic50::Vector{Float64}, turn::Vector{Float64}, kout::Vector{Float64}, kin::Vector{Float64}}, @NamedTuple{conc::Vector{Float64}, pca::Vector{Float64}}, @NamedTuple{}}}}:
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 5, variables: conc, pca

The bettersim syntax provides a reasonable solution for simulating a single dose, now we just need to abstract that code out into a function so that we can apply it to our two remaining doses.
We create a function, simulate_warfarin that accepts a single positional argument, dose that we can use along with map.
- dose was also added as a covariate in the Subject constructor so that we can use it for stratification during post-processing.
- The values for number of subjects (100), and number of samples (200) were hard-coded for simplicity. In a real-world application it would be better to pass those parameters as arguments to simulate_warfarin to improve its overall utility.

function simulate_warfarin(dose)
    #! using literal values for n samples
    _sim = map(1:200) do s
        # create a population
        #! using literal for n subjects
        _pop = map(1:100) do p
            Subject(
                id = p,
                #! using literals for dosing frequency and duration
                events = DosageRegimen(dose; ii = 24, addl = 13),
                #! using literals for wtbl range, adding rep_id and dose as covariates
                covariates = (; wtbl = rand(Uniform(40, 102)), rep_id = s, dose),
            )
        end
        # simulation
        simobs(
            mdl,
            _pop,
            coef(myfit);
            #! using literal for n days in obstimes range
            obstimes = 0:1:(24*14),
            #! no RUV
            simulate_error = false,
        )
    end
end

simulate_warfarin (generic function with 1 method)

Lastly, we perform the simulations using a mapreduce call and save the result in a variable, sim03.
- mapreduce allows us to combine the map and reduce(vcat) steps into a single function call.

sim03 = mapreduce(simulate_warfarin, vcat, [5, 10, 15])

600-element Vector{Vector{Pumas.SimulatedObservations{PumasModel{(tvcl = 1, tvvc = 1, tvka = 1, tvalag = 1, tve0 = 1, tvimax = 1, tvic50 = 1, tvturn = 1, Ωpk = 3, Ωpd = 2, σprop_pk = 1, σadd_pk = 1, σadd_pd = 1), 5, (:depot, :central, :e), (:cl, :vc, :ka, :e0, :imax, :ic50, :turn, :kout, :kin), ParamSet{@NamedTuple{tvcl::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvvc::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvka::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvalag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tve0::RealDomain{Float64, Int64, Int64}, tvimax::RealDomain{Float64, Int64, Float64}, tvic50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, tvturn::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, Ωpk::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, Ωpd::PDiagDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σprop_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pk::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σadd_pd::RealDomain{Float64, TransformVariables.Infinity{true}, Int64}}}, ParamSet{@NamedTuple{ηpk::Pumas.RandomEffect{false, var"#1#12"}, ηpd::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{}, Pumas.ConstantCovar{@NamedTuple{wtbl::Float64, rep_id::Int64, dose::Int64}}, DosageRegimen{Float64, Symbol, Float64, Float64, Float64, Float64}, Nothing}, StepRange{Int64, Int64}, @NamedTuple{tvcl::Float64, tvvc::Float64, tvka::Float64, tvalag::Float64, tve0::Float64, tvimax::Float64, tvic50::Float64, tvturn::Float64, Ωpk::PDMats.PDiagMat{Float64, Vector{Float64}}, Ωpd::PDMats.PDiagMat{Float64, Vector{Float64}}, σprop_pk::Float64, σadd_pk::Float64, σadd_pd::Float64}, @NamedTuple{ηpk::Vector{Float64}, ηpd::Vector{Float64}}, @NamedTuple{wtbl::Vector{Float64}, rep_id::Vector{Int64}, dose::Vector{Int64}}, @NamedTuple{saveat::StepRange{Int64, Int64}, ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, SciMLBase.ODESolution{Float64, 2, Vector{StaticArraysCore.SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{StaticArraysCore.SVector{3, Float64}}}, Nothing, Any, Any, Any, SciMLBase.DEStats, Vector{Int64}, Nothing, Nothing, Nothing}, @NamedTuple{lags::Vector{@NamedTuple{depot::Float64}}}, @NamedTuple{cl::Vector{Float64}, vc::Vector{Float64}, ka::Vector{Float64}, e0::Vector{Float64}, imax::Vector{Float64}, ic50::Vector{Float64}, turn::Vector{Float64}, kout::Vector{Float64}, kin::Vector{Float64}}, @NamedTuple{conc::Vector{Float64}, pca::Vector{Float64}}, @NamedTuple{}}}}:
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 ⋮
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca
 Simulated population (Vector{<:Subject}), n = 100, variables: conc, pca

4.4.1.3 Post-processing

The starting point for post-processing will depend on the output needed to answer the question of interest; in this case a simple tabular summary of metrics with 95% CIs and a graphical summary of 90% PIs for each regimen will suffice.
We will start from a data frame (sim03df)

#! takes ~3-5min on 16vCPU
sim03df = DataFrame(reduce(vcat, (sim03)))

21060000×35 DataFrame

21059975 rows omitted

Row	id	time	conc	pca	evid	lags_depot	amt	cmt	rate	duration	ss	ii	route	wtbl	rep_id	dose	tad	dosenum	depot	central	e	ηpk₁	ηpk₂	ηpk₃	ηpd₁	ηpd₂	cl	vc	ka	e0	imax	ic50	turn	kout	kin
	String	Float64	Float64?	Float64?	Int64?	Float64?	Float64?	Symbol?	Float64?	Float64?	Int8?	Float64?	NCA.Route?	Float64?	Int64?	Int64?	Float64?	Int64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?	Float64?
1	1	0.0	missing	missing	1	0.873414	5.0	depot	0.0	0.0	0	0.0	NullRoute	47.737	1	5	0.0	1	0.0	0.0	94.708	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
2	1	0.0	0.0	94.708	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	0.0	1	0.0	0.0	94.708	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
3	1	1.0	0.149753	94.6588	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	1.0	1	4.31344	0.685737	94.6588	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
4	1	2.0	0.788619	92.8783	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	2.0	1	1.34302	3.61118	92.8783	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
5	1	3.0	0.973958	90.5105	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	3.0	1	0.418161	4.45987	90.5105	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
6	1	4.0	1.01834	88.1335	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	4.0	1	0.130198	4.66308	88.1335	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
7	1	5.0	1.01907	85.8575	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	5.0	1	0.040538	4.66644	85.8575	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
8	1	6.0	1.00645	83.7079	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	6.0	1	0.0126218	4.60866	83.7079	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
9	1	7.0	0.989913	81.6877	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	7.0	1	0.00392991	4.53293	81.6877	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
10	1	8.0	0.972384	79.7933	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	8.0	1	0.00122361	4.45266	79.7933	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
11	1	9.0	0.954774	78.0189	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	9.0	1	0.000380979	4.37202	78.0189	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
12	1	10.0	0.93736	76.3588	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	10.0	1	0.000118621	4.29228	76.3588	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
13	1	11.0	0.920226	74.8071	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	47.737	1	5	11.0	1	3.69335e-5	4.21382	74.8071	-0.184456	-0.173256	-0.00634348	-0.0199407	-0.276096	0.0845647	4.57912	1.16681	94.708	1.0	0.89052	18.8281	0.0531121	5.03014
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
21059989	100	325.0	4.81291	23.7686	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	13.0	14	1.76682e-5	28.2518	23.7686	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059990	100	326.0	4.70712	23.786	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	14.0	14	5.73158e-6	27.6308	23.786	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059991	100	327.0	4.60365	23.8234	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	15.0	14	1.85933e-6	27.0234	23.8234	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059992	100	328.0	4.50245	23.88	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	16.0	14	6.0317e-7	26.4294	23.88	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059993	100	329.0	4.40347	23.9549	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	17.0	14	1.95669e-7	25.8484	23.9549	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059994	100	330.0	4.30668	24.0475	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	18.0	14	6.34753e-8	25.2802	24.0475	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059995	100	331.0	4.21201	24.1571	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	19.0	14	2.05915e-8	24.7245	24.1571	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059996	100	332.0	4.11942	24.283	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	20.0	14	6.6799e-9	24.181	24.283	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059997	100	333.0	4.02886	24.4246	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	21.0	14	2.16697e-9	23.6494	24.4246	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059998	100	334.0	3.9403	24.5813	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	22.0	14	7.02967e-10	23.1296	24.5813	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21059999	100	335.0	3.85369	24.7526	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	23.0	14	2.28043e-10	22.6211	24.7526	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515
21060000	100	336.0	3.76897	24.9378	0	missing	0.0	missing	0.0	0.0	0	0.0	missing	58.0694	200	15	24.0	14	7.39763e-11	22.1239	24.9378	0.102255	-0.12084	-0.0421485	-0.00709995	0.308127	0.130474	5.87	1.12577	95.932	1.0	1.59723	18.8281	0.0531121	5.09515

Sanity Check

We can quickly check the observation-time profile(s) for a single subject to limit the risk of down-stream errors as we continue our analysis. The profile in the figure below appears reasonable.

@chain sim03df begin
    filter(df -> df.id == "1" && df.rep_id == 1, _)
    filter(df -> df.evid == 0, _)
    select(:time, :conc, :pca)
    # default col names for stack are variable and value
    stack(Not([:time]))
    data(_) * mapping(:time, :value, row = :variable) * visual(Lines)
    draw(_; facet = (; linkyaxes = false))
end

Our tabular summary will include three metrics of interest (TA, TTA, TTR) which will be evaluated for each subject, in each simulation.
We will take the average for each metric per simulation and then report the relevant percentiles (2.5, 50, 97.5).
Since there is effectively one assessment per subject, we can make use of the split-apply-combine design for data frames.
The war_metrics function is a custom analysis for table_one; see the SummaryTables.jl documentation for details.

# custom analysis function for table_one
function war_metrics(col)
    all(ismissing, col) && return ("-" => "Median", "-" => "95% CI")
    (
        median(col) => "Median",
        Concat("[", quantile(col, 0.025), ", ", quantile(col, 0.975), "]") => "95% CI",
    )
end

@chain sim03df begin
    # drop records where pca is missing
    dropmissing(_, :pca)
    # first combine step evaluates metrics for individual subjects
    combine(groupby(_, [:dose, :rep_id, :id])) do gdf
        #! metrics 1 and 2
        # find index of first pca value in TR; returns index or nothing 
        i = findfirst(x -> 20 ≤ x < 35, gdf.pca)
        # if `i` was found, return 1 (true), else 0 (false)
        ta_i = Int(!isnothing(i))
        # if no index was found, return missing, else return corresponding time
        tta_i = isnothing(i) ? missing : gdf.time[i]

        #! metric 3
        # temporary df of SS obs from start of Day 14 (312 hours) that are in TR
        _ssdf = filter(df -> df.time >= 312 && 20 ≤ df.pca < 35, gdf)
        # if no obs found, return missing, else return TTR as percentage of ii
        ttr_i =
            iszero(nrow(_ssdf)) ? missing :
            ((last(_ssdf.time) - first(_ssdf.time)) / 24) * 100

        # return a named tuple of the 3 metrics for each subject
        return (; ta_i, tta_i, ttr_i)
    end
    # second combine summarizes each metric per simulation (rep_id)
    combine(
        groupby(_, [:dose, :rep_id]),
        # mean(0|1) * 100 = TA percentage
        :ta_i => (x -> mean(x) * 100) => :ta,
        # applies anonymous function to tta_i and ttr_i cols
        # possible all values could be missing, else could have just used `mean`
        [:tta_i, :ttr_i] .=> function (c)
            all(ismissing, c) && return missing
            mean(skipmissing(c))
        end .=> [:tta, :ttr],
    )
    # from SummaryTables.jl
    table_one(
        _,
        [
            :ta => war_metrics => "Probability of TA",
            :tta => war_metrics => "Time to Target",
            :ttr => war_metrics => "Time in TR",
        ],
        sort = false,
        groupby = :dose => "Dose, mg",
        show_total = false,
    )
end


	Dose, mg
	5	10	15

Probability of TA
Median	31	78	93
95% CI	[22, 39]	[69, 84]	[90, 98]
Time to Target
Median	123	84	62.7
95% CI	[103, 148]	[74.7, 96.6]	[57, 70.9]
Time in TR
Median	95.4	96.8	96.9
95% CI	[86.3, 100]	[92.4, 99.3]	[90.6, 100]

The tabular summary focuses on average response, the graphical summary provide a better understanding of the range of predicted values that we might expect.
We will summarize the relevant pca percentiles (5, 50, 90%) at each time point following the Day 14 dose.

_tbl = @chain sim03df begin
    dropmissing(_, :pca)
    # Day 14 (SS) observations only
    filter(df -> df.time >= 312, _)
    # Summarize by dose and tad; l, m, h are 5th,50th,95th percentile
    combine(
        groupby(_, [:dose, :tad]),
        :pca => (x -> quantile(x, [0.05, 0.5, 0.9])') => [:l, :m, :h],
    )
end

# plot layers
median_layer = mapping(:tad, :m) * visual(Lines; linewidth = 2)
pi_layer = mapping(:tad, :l, :h) * visual(Band, alpha = 0.2)
tr_layer =
    mapping([20, 35], color = "TR" => AlgebraOfGraphics.scale(:secondary)) *
    visual(HLines; linestyle = :dash, alpha = 0.5)

# color and facet map
cf_map = mapping(
    color = :dose => renamer(5 => "5mg", 10 => "10mg", 15 => "15mg") => "",
    col = :dose => renamer(5 => "5mg", 10 => "10mg", 15 => "15mg"),
)

# combine layers and draw
(data(_tbl) * (pi_layer + median_layer) * cf_map) + tr_layer |> draw(
    scales(;
        Y = (; label = "PCA, % of normal"),
        X = (; label = "Time after previous dose, hours"),
        secondary = (; palette = [:gray30]),
    );
    figure = (;
        size = (6, 4) .* 96,
        title = "Predicted PCA-Time Profiles at Steady-state (Day 14) by Dose",
        subtitle = "Median (line), 90%PI (band), TR (dash-line)",
        titlealign = :left,
    ),
    axis = (;
        limits = (0, 24, 0, 80),
        xticks = [0, 12, 24],
        xlabelpadding = 10,
        yticks = 0:20:80,
    ),
    legend = (; orientation = :horizontal, framevisible = false, position = :bottom),
)

4.5 Evaluation

10 mg PO daily dosage regimen offers a reasonable balance of TA, TTA, TTR.

5 Conclusion

Presented the basics of simulation in Pumas using several examples that utilize built-in functionality and user-defined functions.

Reuse

CC BY-SA 4.0