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.05007505416870117
     1     1.762032e+03     1.253034e+03
 * time: 4.190327167510986
     2     1.449413e+03     6.274412e+02
 * time: 5.35685396194458
     3     1.331405e+03     2.099880e+02
 * time: 6.471589088439941
     4     1.310843e+03     1.627399e+02
 * time: 7.6112000942230225
     5     1.295073e+03     1.376531e+02
 * time: 8.847328186035156
     6     1.282335e+03     1.324314e+02
 * time: 10.050434112548828
     7     1.275524e+03     9.521788e+01
 * time: 13.474403142929077
     8     1.271552e+03     5.354883e+01
 * time: 14.838361024856567
     9     1.269315e+03     2.407046e+01
 * time: 16.07956099510193
    10     1.268615e+03     2.368964e+01
 * time: 17.204936981201172
    11     1.268257e+03     2.366048e+01
 * time: 18.364386081695557
    12     1.267352e+03     3.190950e+01
 * time: 19.5216121673584
    13     1.265336e+03     5.399188e+01
 * time: 20.641616106033325
    14     1.260582e+03     8.775888e+01
 * time: 21.692847967147827
    15     1.252929e+03     1.110735e+02
 * time: 22.637086153030396
    16     1.246526e+03     8.328859e+01
 * time: 23.538053035736084
    17     1.240724e+03     8.711807e+01
 * time: 24.38246202468872
    18     1.236881e+03     7.439216e+01
 * time: 25.227177143096924
    19     1.232291e+03     6.878129e+01
 * time: 26.077298164367676
    20     1.227181e+03     5.604034e+01
 * time: 26.90653109550476
    21     1.223093e+03     3.856728e+01
 * time: 27.7409451007843
    22     1.221064e+03     5.186424e+01
 * time: 28.561361074447632
    23     1.220113e+03     4.306874e+01
 * time: 29.352665185928345
    24     1.218936e+03     2.178258e+01
 * time: 30.14394211769104
    25     1.217831e+03     2.144752e+01
 * time: 30.941524982452393
    26     1.216775e+03     2.201560e+01
 * time: 31.74339509010315
    27     1.215938e+03     2.969274e+01
 * time: 32.49789905548096
    28     1.214922e+03     3.438177e+01
 * time: 33.26070499420166
    29     1.212881e+03     4.013028e+01
 * time: 34.03840613365173
    30     1.207231e+03     5.214921e+01
 * time: 34.86295008659363
    31     1.197000e+03     6.567320e+01
 * time: 35.65690612792969
    32     1.189039e+03     6.178092e+01
 * time: 36.50567412376404
    33     1.173765e+03     7.561374e+01
 * time: 37.48718500137329
    34     1.171088e+03     2.385029e+01
 * time: 38.43992209434509
    35     1.170366e+03     2.359334e+01
 * time: 39.27424907684326
    36     1.170030e+03     2.394220e+01
 * time: 40.084580183029175
    37     1.169536e+03     2.287563e+01
 * time: 40.97460103034973
    38     1.168453e+03     1.921641e+01
 * time: 41.806204080581665
    39     1.166697e+03     1.661199e+01
 * time: 42.66723704338074
    40     1.164753e+03     2.146306e+01
 * time: 43.49496603012085
    41     1.163396e+03     1.589737e+01
 * time: 44.32554316520691
    42     1.162810e+03     1.879805e+01
 * time: 45.19176506996155
    43     1.162593e+03     1.922942e+01
 * time: 45.99800109863281
    44     1.162146e+03     1.865903e+01
 * time: 46.78423595428467
    45     1.161045e+03     2.060103e+01
 * time: 47.592491149902344
    46     1.158095e+03     4.011929e+01
 * time: 48.401278018951416
    47     1.150009e+03     6.883383e+01
 * time: 49.251261949539185
    48     1.147406e+03     7.331888e+01
 * time: 50.16746997833252
    49     1.144610e+03     7.504904e+01
 * time: 51.1911940574646
    50     1.138023e+03     6.294975e+01
 * time: 52.08932304382324
    51     1.132705e+03     4.283014e+01
 * time: 52.98611116409302
    52     1.129075e+03     2.327948e+01
 * time: 53.875245094299316
    53     1.128255e+03     2.383372e+01
 * time: 54.798495054244995
    54     1.127697e+03     2.546344e+01
 * time: 55.68593215942383
    55     1.126978e+03     2.642229e+01
 * time: 56.59470200538635
    56     1.126254e+03     2.679647e+01
 * time: 57.58974504470825
    57     1.124269e+03     2.481478e+01
 * time: 58.49312615394592
    58     1.121409e+03     3.925977e+01
 * time: 59.37779211997986
    59     1.118180e+03     3.554728e+01
 * time: 60.273033142089844
    60     1.116311e+03     1.691328e+01
 * time: 61.22349500656128
    61     1.115933e+03     1.725048e+01
 * time: 62.199252128601074
    62     1.115750e+03     1.702130e+01
 * time: 63.18198108673096
    63     1.115157e+03     1.748010e+01
 * time: 64.20085906982422
    64     1.113379e+03     2.744226e+01
 * time: 65.27068901062012
    65     1.107471e+03     4.999759e+01
 * time: 66.46352815628052
    66     1.099217e+03     7.740447e+01
 * time: 68.37810897827148
    67     1.098357e+03     6.707598e+01
 * time: 70.0421371459961
    68     1.093578e+03     3.715572e+01
 * time: 71.55385303497314
    69     1.090146e+03     1.668123e+01
 * time: 73.09645295143127
    70     1.088869e+03     1.082857e+01
 * time: 74.59902715682983
    71     1.087895e+03     9.101157e+00
 * time: 76.26485204696655
    72     1.086966e+03     1.225551e+01
 * time: 77.96845698356628
    73     1.086349e+03     1.200731e+01
 * time: 79.5306830406189
    74     1.086021e+03     9.320681e+00
 * time: 80.95652604103088
    75     1.085858e+03     8.113654e+00
 * time: 82.3368010520935
    76     1.085776e+03     8.127230e+00
 * time: 83.67052912712097
    77     1.085607e+03     7.969659e+00
 * time: 84.9742181301117
    78     1.085219e+03     7.491319e+00
 * time: 86.36067509651184
    79     1.084323e+03     1.248823e+01
 * time: 87.99818301200867
    80     1.082506e+03     2.097609e+01
 * time: 89.38930797576904
    81     1.079392e+03     2.578722e+01
 * time: 90.84561610221863
    82     1.077539e+03     8.916252e+00
 * time: 92.29376006126404
    83     1.077366e+03     1.950846e+00
 * time: 93.71528816223145
    84     1.077341e+03     1.114584e+00
 * time: 95.09158205986023
    85     1.077339e+03     1.131455e+00
 * time: 96.44635605812073
    86     1.077338e+03     1.121808e+00
 * time: 97.77829718589783
    87     1.077335e+03     1.077999e+00
 * time: 99.15225410461426
    88     1.077329e+03     1.045549e+00
 * time: 100.74464797973633
    89     1.077316e+03     1.384916e+00
 * time: 102.23941111564636
    90     1.077292e+03     1.537597e+00
 * time: 103.6275360584259
    91     1.077257e+03     1.232676e+00
 * time: 105.00568795204163
    92     1.077228e+03     5.729974e-01
 * time: 106.41215395927429
    93     1.077219e+03     5.059871e-01
 * time: 107.79241895675659
    94     1.077218e+03     5.033585e-01
 * time: 109.18480610847473
    95     1.077218e+03     5.028759e-01
 * time: 110.55590605735779
    96     1.077218e+03     5.022156e-01
 * time: 111.96864295005798
    97     1.077217e+03     5.008368e-01
 * time: 113.43017506599426
    98     1.077216e+03     4.977187e-01
 * time: 114.91483116149902
    99     1.077212e+03     4.905511e-01
 * time: 116.34241700172424
   100     1.077202e+03     7.604980e-01
 * time: 117.96620416641235
   101     1.077182e+03     1.029997e+00
 * time: 119.52255702018738
   102     1.077153e+03     1.038750e+00
 * time: 121.0212459564209
   103     1.077130e+03     6.003648e-01
 * time: 122.39925503730774
   104     1.077123e+03     5.073911e-01
 * time: 123.91837000846863
   105     1.077122e+03     4.817844e-01
 * time: 125.37690496444702
   106     1.077122e+03     4.675481e-01
 * time: 126.8299970626831
   107     1.077122e+03     4.482693e-01
 * time: 128.2756359577179
   108     1.077121e+03     4.167888e-01
 * time: 129.71754002571106
   109     1.077119e+03     4.159318e-01
 * time: 131.16073298454285
   110     1.077115e+03     7.458266e-01
 * time: 132.5529501438141
   111     1.077106e+03     1.228723e+00
 * time: 133.91985511779785
   112     1.077083e+03     1.875794e+00
 * time: 135.26520609855652
   113     1.077043e+03     2.711417e+00
 * time: 136.62210297584534
   114     1.076990e+03     3.472352e+00
 * time: 137.94868302345276
   115     1.076942e+03     2.744401e+00
 * time: 139.23169016838074
   116     1.076925e+03     9.131364e-01
 * time: 140.52958416938782
   117     1.076924e+03     9.225451e-01
 * time: 141.79222202301025
   118     1.076924e+03     9.211519e-01
 * time: 143.05371403694153
   119     1.076923e+03     9.099674e-01
 * time: 144.29839897155762
   120     1.076920e+03     1.072467e+00
 * time: 145.54616904258728
   121     1.076913e+03     1.901983e+00
 * time: 146.79267001152039
   122     1.076897e+03     2.866979e+00
 * time: 148.03654313087463
   123     1.076863e+03     3.341736e+00
 * time: 149.22282004356384
   124     1.076807e+03     2.150299e+00
 * time: 150.4347381591797
   125     1.076761e+03     2.997934e-01
 * time: 151.6892330646515
   126     1.076749e+03     7.551063e-01
 * time: 152.93441104888916
   127     1.076747e+03     8.968015e-01
 * time: 154.17636799812317
   128     1.076746e+03     7.027920e-01
 * time: 155.38546800613403
   129     1.076744e+03     5.081188e-01
 * time: 156.58852696418762
   130     1.076743e+03     1.635718e-01
 * time: 157.83381700515747
   131     1.076741e+03     1.794518e-01
 * time: 159.03836798667908
   132     1.076740e+03     1.890040e-01
 * time: 160.14994096755981
   133     1.076740e+03     1.923763e-01
 * time: 161.2336790561676
   134     1.076740e+03     1.972342e-01
 * time: 162.29627418518066
   135     1.076740e+03     2.001877e-01
 * time: 163.33357405662537
   136     1.076740e+03     2.031744e-01
 * time: 164.39765095710754
   137     1.076740e+03     2.057411e-01
 * time: 165.57392406463623
   138     1.076740e+03     2.081861e-01
 * time: 166.6868851184845
   139     1.076740e+03     2.083825e-01
 * time: 167.7763409614563
   140     1.076740e+03     2.000834e-01
 * time: 168.88997101783752
   141     1.076739e+03     1.697508e-01
 * time: 170.0373339653015
   142     1.076738e+03     1.576340e-01
 * time: 171.20035696029663
   143     1.076737e+03     7.986184e-02
 * time: 172.37593698501587
   144     1.076737e+03     1.937909e-02
 * time: 173.55480813980103
   145     1.076737e+03     2.909316e-03
 * time: 174.6983940601349
   146     1.076737e+03     8.635752e-04
 * time: 175.84531712532043

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-01-26T15:28:37.266
     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:

  1. Create a population of subjects that has a single covariate (wtbl) that is sampled from a uniform distribution of observed values (40-102 kg).
  2. Simulate an appropriate number of observations (hourly observations will be sufficient).
  3. Repeat the simulation in #2 for a total of 5 simulations.
  4. Store the output in a format that will make post-processing and evaluation as easy as possible.
  5. 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.1936 1 5 0.0 1 0.0 0.0 94.6358 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
2 1 0.0 0.0 94.6358 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 0.0 1 0.0 0.0 94.6358 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
3 1 1.0 0.205142 94.5933 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 1.0 1 4.03962 0.959439 94.5933 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
4 1 2.0 0.899044 93.111 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 2.0 1 0.74919 4.2048 93.111 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
5 1 3.0 1.01501 91.2215 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 3.0 1 0.138945 4.74717 91.2215 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
6 1 4.0 1.02398 89.3729 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 4.0 1 0.0257688 4.78913 89.3729 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
7 1 5.0 1.01329 87.6241 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 5.0 1 0.00477909 4.73914 87.6241 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
8 1 6.0 0.999144 85.9809 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 6.0 1 0.000886331 4.67296 85.9809 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
9 1 7.0 0.984531 84.4398 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 7.0 1 0.000164379 4.60462 84.4398 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
10 1 8.0 0.970009 82.9959 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 8.0 1 3.04858e-5 4.5367 82.9959 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
11 1 9.0 0.955679 81.6442 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 9.0 1 5.65391e-6 4.46968 81.6442 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
12 1 10.0 0.941556 80.3798 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 10.0 1 1.04858e-6 4.40362 80.3798 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
13 1 11.0 0.927641 79.1983 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 47.1936 1 5 11.0 1 1.94469e-7 4.33855 79.1983 -0.370114 -0.140664 0.361095 -0.0207038 0.210987 0.0696353 4.67696 1.68491 94.6358 1.0 1.44937 18.8281 0.0531121 5.02631
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
21059989 100 325.0 2.9301 29.9108 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 13.0 14 3.95531e-13 42.9415 29.9108 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059990 100 326.0 2.88814 29.9292 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 14.0 14 2.94899e-14 42.3266 29.9292 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059991 100 327.0 2.84678 29.9624 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 15.0 14 2.13827e-15 41.7205 29.9624 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059992 100 328.0 2.80602 30.0097 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 16.0 14 -1.77557e-16 41.123 30.0097 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059993 100 329.0 2.76584 30.0704 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 17.0 14 8.45257e-17 40.5342 30.0704 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059994 100 330.0 2.72623 30.144 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 18.0 14 9.06976e-18 39.9538 30.144 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059995 100 331.0 2.68719 30.2298 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 19.0 14 -2.18636e-15 39.3816 30.2298 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059996 100 332.0 2.64871 30.3274 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 20.0 14 -9.90974e-16 38.8177 30.3274 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059997 100 333.0 2.61079 30.4361 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 21.0 14 2.31169e-14 38.2619 30.4361 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059998 100 334.0 2.5734 30.5555 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 22.0 14 4.91953e-14 37.714 30.5555 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21059999 100 335.0 2.53655 30.6851 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 23.0 14 4.85499e-15 37.1739 30.6851 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
21060000 100 336.0 2.50023 30.8244 0 missing 0.0 missing 0.0 0.0 0 0.0 missing 99.8774 200 15 24.0 14 -1.04254e-15 36.6416 30.8244 0.177974 0.251804 0.786492 0.0243247 0.0875191 0.211373 14.6553 2.57825 98.9945 1.0 1.28103 18.8281 0.0531121 5.25781
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 30 78 94
95% CI [22, 38] [69, 85] [89, 98]
Time to Target
Median 124 84.1 63.5
95% CI [102, 150] [75.1, 93.3] [57, 70.4]
Time in TR
Median 95 96.6 97.1
95% CI [86.9, 100] [91.3, 99.4] [91.1, 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.