Pumas.jl Workshop Solutions

Chris Rackauckas, Vijay Ivaturi

Problem 1: Simulate a first-order absorption model with linear elimination after a 100 mg oral dose in 24 subjects

Parameters are: Ka = 1 hr-1, CL = 1 L/hr, V = 20 L/hr.

Part 1: Setup the population

using Pumas, Plots, CSV
single_dose_regimen = DosageRegimen(100, time=0)
first(single_dose_regimen.data)

DataFrameRow

1 rows × 8 columns

timecmtamtevidiiaddlratess
Float64Int64Float64Int8Float64Int64Float64Int8
10.01100.010.000.00

to build a sinlge subject

s1 = Subject(id=1, evs=single_dose_regimen,cvs=(Wt=70,))
Subject
  ID: 1
  Events: 1

let's first define a function to choose body weight randomly

choose_covariates() = (Wt = rand(55:80),)
choose_covariates (generic function with 1 method)

Then, we use generate a population of subjects with a random weight generated from the covariate function above

pop = Population(map(i -> Subject(id = i,evs = single_dose_regimen, cvs =  choose_covariates()),1:24))
Population
  Subjects: 24
  Covariates: Wt

You can view the generated population using by calling a random subject by index and look at the subject's

  • covariates

  • events

  • id numbers

  • observations

  • time

Let us us peek at the first subject's covariates

pop[1].covariates
(Wt = 77,)

Part 2: Write the model

mymodel = @model begin
  @param   begin
    tvcl  RealDomain(lower=0, init = 1.0)
    tvv  RealDomain(lower=0, init = 20)
    tvka  RealDomain(lower = 0, init= 1)
    Ω  PDiagDomain(init=[0.09,0.09, 0.09])
    σ_prop  RealDomain(lower=0,init=0.04)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    CL = tvcl * (Wt/70)^0.75 * exp(η[1])
    V  = tvv * (Wt/70) * exp(η[2])
    Ka = tvka * exp(η[3])
  end
  @covariates Wt

  @dynamics OneCompartmentModel
    #@dynamics begin
    #    Depot' =  -Ka*Depot
    #    Central' =  Ka*Depot - (CL/V)*Central
    #end

  @derived begin
      cp = @. 1000*(Central / V)
      dv ~ @. Normal(cp, sqrt(cp^2*σ_prop))
    end
end
PumasModel
  Parameters: tvcl, tvv, tvka, Ω, σ_prop
  Random effects: η
  Covariates: Wt
  Dynamical variables: Depot, Central
  Derived: cp, dv
  Observed: cp, dv

Note that above, we are using the analytical solution in @dynamics. You can switch to using the differential equation system if you prefer.

Part 3: Simulate

Let's first extract the model parameters

param = init_param(mymodel)
(tvcl = 1.0, tvv = 20.0, tvka = 1.0, Ω = PDMats.PDiagMat{Float64,Array{Floa
t64,1}}(3, [0.09, 0.09, 0.09], [11.1111, 11.1111, 11.1111]), σ_prop = 0.04)

Then using the simobs function, carry out the simulation and visualize the simulation output

obs = simobs(mymodel, pop, param, obstimes=0:1:72)
plot(obs)

where

  • mymodel is the model setup in the Part 2,

  • pop is the population of subjects that was setup in Part 1

  • param is the specified set of model parameters

  • obstimes specifies the simulation time period.

Problem 2: Peform Non-compartmental analysis

We will start by generating a dataframe of the resuls from the simulation step

simdf = DataFrame(obs)
first(simdf, 6)

6 rows × 9 columns

idtimecpdvamtevidcmtrateWt
StringInt64Float64Float64Float64Int8Int64⍰Float64Int64
1100.00.0100.0110.077
2100.00.00.00missing0.077
3112175.852026.860.00missing0.077
4122958.332869.170.00missing0.077
5133189.372446.590.00missing0.077
6143203.973142.270.00missing0.077

For the purpose of NCA, let us use the cp (output without residual error) as our observed value

To prepare the dataset for NCA analysis, let us use the read_nca function. The NCA datasets in Pumas requires a route specification which can either be iv or ev. Since this is an oral drug administration, lets add that to the simdf.

simdf.route = "ev"
"ev"

Next we can define time, concentration and dose units so the report includes the units for the pharmacokinetic parameters. The general syntax for units are u followed by the unit in quotes "".

timeu = u"hr"
concu = u"mg/L"
amtu  = u"mg"
mg
ncadf = read_nca(simdf, id=:id, time=:time, conc=:cp, amt=:amt,
    route=:route,timeu=timeu, concu=concu, amtu=amtu, lloq=0.4concu)
NCAPopulation (24 subjects):
  ID: ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
 "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24"]
    concentration: mg L^-1
    time:          hr
    auc:           mg hr L^-1
    aumc:          mg hr^2 L^-1
    λz:            hr^-1
    dose:          mg

You can view the concentration-time plots by doing

plot(ncadf)

You can then generate cmax and auc for each subject

auc = NCA.auc(ncadf)

24 rows × 2 columns

idauc
StringUnitful…
1193992.0 mg hr L^-1
222.1821e5 mg hr L^-1
331.38876e5 mg hr L^-1
4464345.1 mg hr L^-1
551.23336e5 mg hr L^-1
661.02345e5 mg hr L^-1
771.13016e5 mg hr L^-1
881.32914e5 mg hr L^-1
991.30535e5 mg hr L^-1
10101.23127e5 mg hr L^-1
111173594.0 mg hr L^-1
121270247.9 mg hr L^-1
131398589.3 mg hr L^-1
14141.32914e5 mg hr L^-1
15151.63776e5 mg hr L^-1
161678982.6 mg hr L^-1
17171.15825e5 mg hr L^-1
181867401.4 mg hr L^-1
191967233.4 mg hr L^-1
20201.12205e5 mg hr L^-1
21211.91297e5 mg hr L^-1
22221.40436e5 mg hr L^-1
23231.2803e5 mg hr L^-1
24241.40759e5 mg hr L^-1
cmax = NCA.cmax(ncadf)

24 rows × 2 columns

idcmax
StringUnitful…
113203.97 mg L^-1
224056.79 mg L^-1
336055.62 mg L^-1
442513.61 mg L^-1
557740.42 mg L^-1
665667.73 mg L^-1
774206.79 mg L^-1
889103.61 mg L^-1
995587.37 mg L^-1
10106366.82 mg L^-1
11114473.4 mg L^-1
12124693.58 mg L^-1
13133520.14 mg L^-1
14149103.61 mg L^-1
15152569.9 mg L^-1
16167324.92 mg L^-1
17171906.77 mg L^-1
18183027.34 mg L^-1
19194291.22 mg L^-1
20207303.68 mg L^-1
21215513.51 mg L^-1
22227661.44 mg L^-1
23232269.45 mg L^-1
24242947.55 mg L^-1

Or generate the entire NCA report using

report = NCAReport(ncadf)
report = NCA.to_dataframe(report)
first(report,6)

6 rows × 38 columns

iddoseamtlambda_zhalf_lifetmaxtlagcmaxclastclast_predauclasttlastaucinf_obsvz_f_obscl_f_obsaucinf_predvz_f_predcl_f_predcmax_dauclast_daucinf_d_obsauc_extrap_obsaucinf_d_predauc_extrap_predaumclastaumcinf_obsaumc_extrap_obsaumcinf_predaumc_extrap_predn_samplesrsqrsq_adjustedcorr_xyno_points_lambda_zlambda_z_interceptlambda_z_lowerlambda_z_upperspanroute
StringUnitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Unitful…Float64Unitful…Float64Unitful…Unitful…Float64Unitful…Float64Int64Float64Float64Float64Int64Float64Unitful…Unitful…Float64String
11100.0 mg0.0392181 hr^-117.6742 hr4 hr0 hr3203.97 mg L^-1229.324 mg L^-1229.324 mg L^-188144.6 mg hr L^-172 hr93992.0 mg hr L^-10.0271283 L0.00106392 L hr^-193992.0 mg hr L^-10.0271283 L0.00106392 L hr^-13203.97 mg L^-1881.446 hr L^-1939.92 hr L^-16.22117939.92 hr L^-16.221171.93599e6 mg hr^2 L^-12.5061e6 mg hr^2 L^-122.7492.5061e6 mg hr^2 L^-122.749731.01.01.038.2588470 hr72 hr0.113159EV
22100.0 mg0.0203249 hr^-134.1034 hr4 hr0 hr4056.79 mg L^-11052.35 mg L^-11052.35 mg L^-11.66434e5 mg hr L^-172 hr2.1821e5 mg hr L^-10.0225474 L0.000458273 L hr^-12.1821e5 mg hr L^-10.0225474 L0.000458273 L hr^-14056.79 mg L^-11664.34 hr L^-12182.1 hr L^-123.72772182.1 hr L^-123.72774.72515e6 mg hr^2 L^-11.10005e7 mg hr^2 L^-157.0461.10005e7 mg hr^2 L^-157.046731.01.01.038.4221770 hr72 hr0.0586452EV
33100.0 mg0.0514102 hr^-113.4827 hr3 hr0 hr6055.62 mg L^-1187.095 mg L^-1187.095 mg L^-11.35237e5 mg hr L^-172 hr1.38876e5 mg hr L^-10.0140063 L0.000720067 L hr^-11.38876e5 mg hr L^-10.0140063 L0.000720067 L hr^-16055.62 mg L^-11352.37 hr L^-11388.76 hr L^-12.620521388.76 hr L^-12.620522.5263e6 mg hr^2 L^-12.85912e6 mg hr^2 L^-111.64052.85912e6 mg hr^2 L^-111.6405731.01.01.038.9331570 hr72 hr0.148338EV
44100.0 mg0.0472238 hr^-114.6779 hr4 hr0 hr2513.61 mg L^-1108.884 mg L^-1108.884 mg L^-162039.4 mg hr L^-172 hr64345.1 mg hr L^-10.0329097 L0.00155412 L hr^-164345.1 mg hr L^-10.0329097 L0.00155412 L hr^-12513.61 mg L^-1620.394 hr L^-1643.451 hr L^-13.58334643.451 hr L^-13.583341.24189e6 mg hr^2 L^-11.45673e6 mg hr^2 L^-114.74781.45673e6 mg hr^2 L^-114.7478731.01.01.038.090470 hr72 hr0.136259EV
55100.0 mg0.0748125 hr^-19.26512 hr2 hr0 hr7740.42 mg L^-145.2144 mg L^-145.2144 mg L^-11.22732e5 mg hr L^-172 hr1.23336e5 mg hr L^-10.0108376 L0.00081079 L hr^-11.23336e5 mg hr L^-10.0108376 L0.00081079 L hr^-17740.42 mg L^-11227.32 hr L^-11233.36 hr L^-10.4900171233.36 hr L^-10.4900171.70694e6 mg hr^2 L^-11.75853e6 mg hr^2 L^-12.933871.75853e6 mg hr^2 L^-12.93387731.01.01.039.1979270 hr72 hr0.215863EV
66100.0 mg0.0704468 hr^-19.8393 hr3 hr0 hr5667.73 mg L^-149.9496 mg L^-149.9496 mg L^-11.01636e5 mg hr L^-172 hr1.02345e5 mg hr L^-10.0138699 L0.000977091 L hr^-11.02345e5 mg hr L^-10.0138699 L0.000977091 L hr^-15667.73 mg L^-11016.36 hr L^-11023.45 hr L^-10.6927971023.45 hr L^-10.6927971.53109e6 mg hr^2 L^-11.5922e6 mg hr^2 L^-13.838441.5922e6 mg hr^2 L^-13.83844731.01.01.038.9831870 hr72 hr0.203266EV

Problem 3: Estimate using Non-linear mixed effects

We can use the simulated dataset in the Problem 1 for our estimation. We need a couple of data manipulation steps

  1. missing cmt should be converted to 2 to reflect central compartment

  2. data rows where time = 0, and cp=0 should be removed

simdf.cmt = ifelse.(ismissing.(simdf.cmt), 2, simdf.cmt)
est_df = simdf[.!((simdf.dv .== 0.0) .& (simdf.cmt .==2)),:]
first(est_df,6)

6 rows × 10 columns

idtimecpdvamtevidcmtrateWtroute
StringInt64Float64Float64Float64Int8Int64Float64Int64String
1100.00.0100.0110.077ev
2112175.852026.860.0020.077ev
3122958.332869.170.0020.077ev
4133189.372446.590.0020.077ev
5143203.973142.270.0020.077ev
6153135.394202.210.0020.077ev

Part 1: Read datasets for NLME estimation

We can use the read_pumas function to prepare the dataset for NLME estimation

data = read_pumas(est_df ,cvs = [:Wt], dvs=[:dv])
Population
  Subjects: 24
  Covariates: Wt
  Observables: dv

where

  • cvs takes an array of covariates

  • dvs takes an array of the dependent variables

  • since the dataframe has time as the variable, the function does not need a specific input

Part 2: Perform a model fit

We now use the

  • mymodel model that we wrote earlier

  • the set of parameters specified in param as initial estimates

  • data that was read in using the read_pumas function

to fit the model.

res = fit(mymodel,data,param,Pumas.FOCEI())
FittedPumasModel

Successful minimization:                true

Likelihood approximation:        Pumas.FOCEI
Objective function value:            11205.8
Total number of observation records:    1728
Number of active observation records:   1728
Number of subjects:                       24

-----------------
       Estimate
-----------------
tvcl    0.93689
tvv    19.93
tvka    0.97408
Ω₁,₁    0.10496
Ω₂,₂    0.20987
Ω₃,₃    0.077789
σ_prop  0.037023
-----------------

Part 3: Infer the results

infer provides the model inference

infer(res)
Calculating: variance-covariance matrix. Done.
FittedPumasModelInference

Successful minimization:                true

Likelihood approximation:        Pumas.FOCEI
Objective function value:            11205.8
Total number of observation records:    1728
Number of active observation records:   1728
Number of subjects:                       24

-------------------------------------------------
       Estimate      RSE          95.0% C.I.
-------------------------------------------------
tvcl    0.93689    6.7391[ 0.81315  ;  1.0606  ]
tvv    19.93       9.0617[16.39     ; 23.469   ]
tvka    0.97408    7.9653[ 0.82201  ;  1.1261  ]
Ω₁,₁    0.10496   22.65  [ 0.058367 ;  0.15156 ]
Ω₂,₂    0.20987   23.142 [ 0.11468  ;  0.30506 ]
Ω₃,₃    0.077789  47.259 [ 0.0057359;  0.14984 ]
σ_prop  0.037023   2.5185[ 0.035196 ;  0.038851]
-------------------------------------------------

Part 4: Inspect the results

inspect gives you the

  • model predictions

  • residuals

  • Empirical Bayes estimates

preds = DataFrame(predict(res))
first(preds, 6)

6 rows × 6 columns

idtimeWtpredipredpred_approx
StringFloat64Int64Float64Float64Pumas…
111.0772667.192044.85FOCEI()
212.0773587.222827.5FOCEI()
313.0773831.093081.79FOCEI()
414.0773817.83116.55FOCEI()
515.0773710.173061.64FOCEI()
616.0773570.972971.66FOCEI()
resids = DataFrame(wresiduals(res))
first(resids, 6)

6 rows × 6 columns

idtimeWtwresiwreswres_approx
StringFloat64Int64Float64Float64Pumas…
111.077-0.604593-0.0457289FOCEI()
212.077-0.02076410.0765949FOCEI()
313.077-0.905388-1.0712FOCEI()
414.0770.357970.0429013FOCEI()
515.0771.91781.9361FOCEI()
616.077-0.585382-0.367784FOCEI()
ebes = DataFrame(empirical_bayes(res))
first(ebes, 6)

6 rows × 7 columns

idtimeWtebe_1ebe_2ebe_3ebes_approx
StringFloat64Int64Float64Float64Float64Pumas…
111.0770.07000940.233745-0.118514FOCEI()
212.0770.07000940.233745-0.118514FOCEI()
313.0770.07000940.233745-0.118514FOCEI()
414.0770.07000940.233745-0.118514FOCEI()
515.0770.07000940.233745-0.118514FOCEI()
616.0770.07000940.233745-0.118514FOCEI()

There is an inspect function that provides all the results at once

Note that this function below fails to convert into a dataframe due to a bug. Will be fixed soon

resout = DataFrame(inspect(res))
first(resout, 6)

Problem 4: Validate your model

Finally validate your model with a visual predictive check

vpc(res,200) |> plot

or stratify it based on bodyweight

vpc(res,200, stratify_on=[:Wt]) |> plot