Pumas.jl Workshop Solutions

PumasAI
August 2020

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, Random
Random.seed!(0)
Random.MersenneTwister(UInt32[0x00000000], Random.DSFMT.DSFMT_state(Int32[7
48398797, 1073523691, -1738140313, 1073664641, -1492392947, 1073490074, -16
25281839, 1073254801, 1875112882, 1073717145  …  943540191, 1073626624, 109
1647724, 1073372234, -1273625233, -823628301, 835224507, 991807863, 382, 0]
), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0
, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000
, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0
x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00
000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000
000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000
000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000
000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000
000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000
000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000
000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000
000000000000], 1002, 0)
single_dose_regimen = DosageRegimen(100)
first(single_dose_regimen.data)

DataFrameRow (9 columns)

timecmtamtevidiiaddlratedurationss
Float64Int64Float64Int8Float64Int64Float64Float64Int8
10.01100.010.000.00.00

to build a single subject

s1 = Subject(id=1, events=single_dose_regimen, covariates=(Wt=70,))
Subject
  ID: 1
  Events: 1
  Covariates: Wt

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 it to generate a population of subjects with a random weight generated from the covariate function above

pop = Population(map(i -> Subject(id = i, events=single_dose_regimen, covariates=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
Pumas.ConstantCovar{NamedTuple{(:Wt,),Tuple{Int64}}}((Wt = 74,))

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])
    Vc  = tvv * (Wt/70) * exp(η[2])
    Ka = tvka * exp(η[3])
  end
  @covariates Wt

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

  @derived begin
      cp = @. 1000*(Central / Vc)
      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, tvka = 1, Ω = PDMats.PDiagMat{Float64,Array{Float64,
1}}(3, [0.09, 0.09, 0.09], [11.11111111111111, 11.11111111111111, 11.111111
11111111]), σ_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
StringFloat64Float64?Float64?Float64Int8Int64?Float64Int64
110.0missingmissing100.0110.074
210.00.00.00.00missing0.074
311.02136.141928.230.00missing0.074
412.02803.413018.930.00missing0.074
513.02977.333186.010.00missing0.074
614.02986.352872.350.00missing0.074

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"
1776-element Array{String,1}:
 "ev"
 "ev"
 "ev"
 "ev"
 "ev"
 "ev"
 "ev"
 "ev"
 "ev"
 "ev"
 ⋮
 "ev"
 "ev"
 "ev"
 "ev"
 "ev"
 "ev"
 "ev"
 "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⁻¹
    time:          hr
    auc:           mg hr L⁻¹
    aumc:          mg hr² L⁻¹
    λz:            hr⁻¹
    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
StringQuantit…
11132198.0 mg hr L⁻¹
22102079.0 mg hr L⁻¹
33106989.0 mg hr L⁻¹
4457458.7 mg hr L⁻¹
5571282.8 mg hr L⁻¹
6696739.3 mg hr L⁻¹
77161621.0 mg hr L⁻¹
8872576.2 mg hr L⁻¹
9961923.5 mg hr L⁻¹
1010132283.0 mg hr L⁻¹
1111157841.0 mg hr L⁻¹
1212124295.0 mg hr L⁻¹
1313109570.0 mg hr L⁻¹
1414110130.0 mg hr L⁻¹
151583332.2 mg hr L⁻¹
1616118001.0 mg hr L⁻¹
171781171.9 mg hr L⁻¹
1818119144.0 mg hr L⁻¹
1919210948.0 mg hr L⁻¹
202090731.3 mg hr L⁻¹
2121169402.0 mg hr L⁻¹
2222139547.0 mg hr L⁻¹
2323124435.0 mg hr L⁻¹
2424108008.0 mg hr L⁻¹
cmax = NCA.cmax(ncadf)

24 rows × 2 columns

idcmax
StringQuantit…
112986.35 mg L⁻¹
224845.71 mg L⁻¹
333300.26 mg L⁻¹
443188.57 mg L⁻¹
554587.95 mg L⁻¹
663826.88 mg L⁻¹
773786.89 mg L⁻¹
883206.47 mg L⁻¹
994795.85 mg L⁻¹
10105447.36 mg L⁻¹
11112948.35 mg L⁻¹
12127299.77 mg L⁻¹
13134210.48 mg L⁻¹
14144776.76 mg L⁻¹
15156813.03 mg L⁻¹
16165285.48 mg L⁻¹
17173288.68 mg L⁻¹
18184481.94 mg L⁻¹
19198022.22 mg L⁻¹
20204503.52 mg L⁻¹
21213520.06 mg L⁻¹
22225590.86 mg L⁻¹
23236493.77 mg L⁻¹
24245435.34 mg L⁻¹

Or generate the entire NCA report using

report = NCAReport(ncadf)
first(report,6)

6 rows × 39 columns (omitted printing of 32 columns)

iddoseamttlagtmaxcmaxtlastclast
StringQuantit…Quantit…Quantit…Quantit…Quantit…Quantit…
11100.0 mg0.0 hr4.0 hr2986.35 mg L⁻¹72.0 hr566.15 mg L⁻¹
22100.0 mg0.0 hr4.0 hr4845.71 mg L⁻¹72.0 hr93.876 mg L⁻¹
33100.0 mg0.0 hr3.0 hr3300.26 mg L⁻¹72.0 hr326.361 mg L⁻¹
44100.0 mg0.0 hr4.0 hr3188.57 mg L⁻¹72.0 hr26.0065 mg L⁻¹
55100.0 mg0.0 hr3.0 hr4587.95 mg L⁻¹72.0 hr20.3659 mg L⁻¹
66100.0 mg0.0 hr4.0 hr3826.88 mg L⁻¹72.0 hr166.316 mg L⁻¹

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
StringFloat64Float64?Float64?Float64Int8Int64Float64Int64String
110.0missingmissing100.0110.074ev
211.02136.141928.230.0020.074ev
312.02803.413018.930.0020.074ev
413.02977.333186.010.0020.074ev
514.02986.352872.350.0020.074ev
615.02941.062521.280.0020.074ev

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, covariates=[:Wt], observations=[:dv])
Population
  Subjects: 24
  Covariates: Wt
  Observables: dv

where

  • covariates is a symbol that maps to the columns of covariates

  • observations is a symbol that maps to the columns of 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())
Iter     Function value   Gradient norm 
     0     1.134964e+04     1.928833e+01
 * time: 7.700920104980469e-5
     1     1.134928e+04     3.577934e+01
 * time: 0.14695405960083008
     2     1.134867e+04     2.016275e+01
 * time: 0.2311711311340332
     3     1.134846e+04     1.467966e+01
 * time: 0.3096771240234375
     4     1.134813e+04     2.651939e+00
 * time: 0.4403500556945801
     5     1.134812e+04     2.339406e+00
 * time: 0.5061051845550537
     6     1.134809e+04     1.934372e+00
 * time: 0.5611600875854492
     7     1.134805e+04     3.034282e-01
 * time: 0.6209461688995361
     8     1.134805e+04     2.558618e-01
 * time: 0.7131052017211914
     9     1.134804e+04     8.600182e-02
 * time: 0.7634930610656738
    10     1.134804e+04     2.772043e-02
 * time: 0.8129270076751709
    11     1.134804e+04     9.824259e-03
 * time: 0.8591082096099854
    12     1.134804e+04     4.457464e-03
 * time: 0.9450790882110596
    13     1.134804e+04     1.025930e-03
 * time: 0.988976001739502
    14     1.134804e+04     8.285445e-05
 * time: 1.0336220264434814
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:              Pumas.FOCEI
Log-likelihood value:                   -11348.038
Number of subjects:                             24
Number of parameters:         Fixed      Optimized
                                  0              5
Observation records:         Active        Missing
    dv:                        1728              0
    Total:                     1728              0

---------------------
           Estimate
---------------------
tvcl        0.96606
tvv        20.124
tvka        0.89507
Ω₁,₁        0.10253
Ω₂,₂        0.077902
Ω₃,₃        0.074778
σ_prop      0.038981
---------------------

Part 3: Infer the results

infer provides the model inference

infer(res)
Calculating: variance-covariance matrix. Done.
Asymptotic inference results

Successful minimization:                      true

Likelihood approximation:              Pumas.FOCEI
Log-likelihood value:                   -11348.038
Number of subjects:                             24
Number of parameters:         Fixed      Optimized
                                  0              5
Observation records:         Active        Missing
    dv:                        1728              0
    Total:                     1728              0

--------------------------------------------------------------------
          Estimate           SE                      95.0% C.I.
--------------------------------------------------------------------
tvcl       0.96606         0.063144         [ 0.8423   ;  1.0898  ]
tvv       20.124           1.1628           [17.845    ; 22.403   ]
tvka       0.89507         0.066005         [ 0.7657   ;  1.0244  ]
Ω₁,₁       0.10253         0.026178         [ 0.05122  ;  0.15384 ]
Ω₂,₂       0.077902        0.020065         [ 0.038576 ;  0.11723 ]
Ω₃,₃       0.074778        0.04165          [-0.0068549;  0.15641 ]
σ_prop     0.038981        0.0015704        [ 0.035903 ;  0.042059]
--------------------------------------------------------------------

Part 4: Inspect the results

inspect gives you the

  • model predictions

  • residuals

  • Empirical Bayes estimates

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

6 rows × 7 columns

idtimeeviddvdv_preddv_ipredWt
StringFloat64Int64Float64Float64Float64Int64
111.001928.232705.811956.574
212.003018.933686.232724.7574
313.003186.013967.482997.5374
414.002872.353968.573064.4174
515.002521.283860.473046.4174
616.003474.893712.772994.174
resids = DataFrame(wresiduals(res))
first(resids, 6)

6 rows × 6 columns

idtimedv_wresdv_iwreswres_approxWt
StringFloat64Float64Float64FOCEI…Int64
111.0-0.411365-0.073188FOCEI()74
212.00.1879170.546827FOCEI()74
313.0-0.02812960.318467FOCEI()74
414.0-0.642131-0.317451FOCEI()74
515.0-1.17404-0.87309FOCEI()74
616.00.5355750.813335FOCEI()74
ebes = empirical_bayes(res)
24-element Array{NamedTuple{(:η,),Tuple{Array{Float64,1}}},1}:
 (η = [-0.2904324910014926, 0.32051585069031013, -0.025910521490302388],)
 (η = [0.15199534709748194, 0.006694677876328507, -0.08505621551491206],)
 (η = [0.031515603554869896, 0.4075062260692816, -0.08656419062156831],)
 (η = [0.6880578364465235, 0.30321844824353766, -0.3622518118272251],)
 (η = [0.33469295296182006, -0.1832827147147013, 0.16122027365885316],)
 (η = [0.08295344196269373, 0.12360173157291095, -0.1729395453249503],)
 (η = [-0.511131927491953, 0.11194536524752642, -0.41615521162861957],)
 (η = [0.2718347359378102, 0.16701157563252736, 0.14320197042334237],)
 (η = [0.5927403091418617, -0.10695202878979747, 0.04688600142116455],)
 (η = [-0.17845655253482828, -0.15779252808464944, -0.06272480097206558],)
 ⋮
 (η = [-0.01110822809004068, -0.10290139041393323, -0.26598794317569],)
 (η = [0.2700672576227783, 0.3799927856330493, 0.31484327694740916],)
 (η = [-0.12203376166661534, 0.023553558052422405, 0.28187878876154926],)
 (η = [-0.5312831786274118, -0.37389475635153313, 0.2347824182097107],)
 (η = [0.028003309206508557, -0.23279073374554443, 0.0009664479429166229],)
 (η = [-0.43588668749431664, 0.3543966457784553, 0.053406350455314465],)
 (η = [-0.26436530654479284, -0.2302933903372198, 0.0378864488734128],)
 (η = [0.019386224183180666, -0.19523044653084168, 0.03783817560327249],)
 (η = [-0.025519566815914455, -0.2973288415844167, 0.01380599186212406],)

The inspect function provides all the results at once

resout = DataFrame(inspect(res))
first(resout, 6)
Calculating: predictions, weighted residuals, empirical bayes. Done.

6 rows × 13 columns (omitted printing of 4 columns)

idtimeeviddvdv_preddv_ipreddv_wresdv_iwreswres_approx
StringFloat64Int64Float64Float64Float64Float64Float64FOCEI…
111.001928.232705.811956.5-0.411365-0.073188FOCEI()
212.003018.933686.232724.750.1879170.546827FOCEI()
313.003186.013967.482997.53-0.02812960.318467FOCEI()
414.002872.353968.573064.41-0.642131-0.317451FOCEI()
515.002521.283860.473046.41-1.17404-0.87309FOCEI()
616.003474.893712.772994.10.5355750.813335FOCEI()