Introduction to Pumas

Chris Rackauckas
July 19, 2019

Introduction

This is an introduction to Pumas, a software for pharmacometric modeling and simulation.

The basic workflow of Pumas is:

  1. Build a model.

  2. Define subjects or populations to simulate or estimate.

  3. Analyze the results with post-processing and plots.

We will show how to build a multiple-response PK/PD model via the @model macro, define a subject with multiple doses, and analyze the results of the simulation. This tutorial is made to be a broad overview of the workflow and more in-depth treatment of each section can be found in the subsequent tutorials and documentation.

Installation

To install Pumas, first install JuliaPro. Then use the command:

using Pkg
Pkg.add("Pumas")

Doing it this way, Pumas and its dependencies will install automatically.

Getting Started

To load the package, use

using Pumas

Using the Model Macro

Now let's define a model. A model is defined in an @model block. Inside of this block we have a few subsections. The first of which is @param. In here we define what kind of parameters we have. For this model we will define a vector parameter θ of size 12:

@param begin
  θ  VectorDomain(12)
end

Next we define our random effects. The random effects are defined by a distribution from Distributions.jl. For more information on defining distributions, please see the Distributions.jl documentation. For this tutorial, we wish to have a multivariate normal of 11 uncorrelated random effects, so we utilize the syntax:

using LinearAlgebra
@random begin
  η ~ MvNormal(Matrix{Float64}(I, 11, 11))
end

Notice that here we imported I from LinearAlgebra and and said that our Normal distribution's covariance is said I, the identity matrix.

Now we define our pre-processing step in @pre. This is where we choose how the parameters, random effects, and the covariates collate. We define the values and give them a name as follows:

@pre begin
    Ka1     = θ[1]
    CL      = θ[2]*exp(η[1])
    Vc      = θ[3]*exp(η[2])
    Q       = θ[4]*exp(η[3])
    Vp      = θ[5]*exp(η[4])
    Kin     = θ[6]*exp(η[5])
    Kout    = θ[7]*exp(η[6])
    IC50    = θ[8]*exp(η[7])
    IMAX    = θ[9]*exp(η[8])
    γ       = θ[10]*exp(η[9])
    Vmax    = θ[11]*exp(η[10])
    Km      = θ[12]*exp(η[11])
end

Next we define the @init block which gives the inital values for our differential equations. Any variable not mentioned in this block is assumed to have a zero for its starting value. We wish to only set the starting value for Resp, and thus we use:

@init begin
    Resp = θ[6]/θ[7]
end

Now we define our dynamics. We do this via the @dynamics block. Differential variables are declared by having a line defining their derivative. For our model, we use:

@dynamics begin
    Ev1'    = -Ka1*Ev1
    Cent'   =  Ka1*Ev1 - (CL+Vmax/(Km+(Cent/Vc))+Q)*(Cent/Vc)  + Q*(Periph/Vp)
    Periph' =  Q*(Cent/Vc)  - Q*(Periph/Vp)
    Resp'   =  Kin*(1-(IMAX*(Cent/Vc)^γ/(IC50^γ+(Cent/Vc)^γ)))  - Kout*Resp
end

Lastly we utilize the @derived macro to define our post-processing. We can output values using the following:

@derived begin
    ev1    = Ev1
    cp     = Cent / θ[3]
    periph = Periph
    resp   = Resp
end

The @model block is all of these together, giving us the following model:

using LinearAlgebra
model = @model begin

    @param begin
      θ  VectorDomain(12)
    end

    @random begin
      η ~ MvNormal(Matrix{Float64}(I, 11, 11))
    end

    @pre begin
        Ka1     = θ[1]
        CL      = θ[2]*exp(η[1])
        Vc      = θ[3]*exp(η[2])
        Q       = θ[4]*exp(η[3])
        Vp      = θ[5]*exp(η[4])
        Kin     = θ[6]*exp(η[5])
        Kout    = θ[7]*exp(η[6])
        IC50    = θ[8]*exp(η[7])
        IMAX    = θ[9]*exp(η[8])
        γ       = θ[10]*exp(η[9])
        Vmax    = θ[11]*exp(η[10])
        Km      = θ[12]*exp(η[11])
    end

    @init begin
        Resp = θ[6]/θ[7]
    end

    @dynamics begin
        Ev1'    = -Ka1*Ev1
        Cent'   =  Ka1*Ev1 - (CL+Vmax/(Km+(Cent/Vc))+Q)*(Cent/Vc)  + Q*(Periph/Vp)
        Periph' =  Q*(Cent/Vc)  - Q*(Periph/Vp)
        Resp'   =  Kin*(1-(IMAX*(Cent/Vc)^γ/(IC50^γ+(Cent/Vc)^γ)))  - Kout*Resp
    end

    @derived begin
        ev1    = Ev1
        cp     = Cent / θ[3]
        periph = Periph
        resp   = Resp
    end
end
PumasModel
  Parameters: θ
  Random effects: η
  Covariates: 
  Dynamical variables: Ev1, Cent, Periph, Resp
  Derived: ev1, cp, periph, resp
  Observed: ev1, cp, periph, resp

Building a Subject

Now let's build a subject to simulate the model with. A subject defines three components:

  1. The dosage regimen

  2. The covariates of the indvidual

  3. Observations associated with the individual.

Our model did not make use of covariates so we will ignore (2) for now, and (3) is only necessary for fitting parameters to data which will not be covered in this tutorial. Thus our subject will be defined simply by its dosage regimen.

To do this, we use the DosageRegimen constructor. It uses terms from the NMTRAN format to specify its dose schedule. The first value is always the dosing amount. Then there are optional arguments, the most important of which is time which specifies the time that the dosing occurs. For example,

DosageRegimen(15, time=0)
Pumas.DosageRegimen(1×8 DataFrames.DataFrame
│ Row │ time    │ cmt   │ amt     │ evid │ ii      │ addl  │ rate    │ ss  
 │
│     │ Float64 │ Int64 │ Float64 │ Int8 │ Float64 │ Int64 │ Float64 │ Int8
 │
├─────┼─────────┼───────┼─────────┼──────┼─────────┼───────┼─────────┼─────
─┤
│ 1   │ 0.0     │ 1     │ 15.0    │ 1    │ 0.0     │ 0     │ 0.0     │ 0   
 │)

is a dosage regimen which simply does a single dose at time t=0 of amount 15. If we use arrays, then the dosage regimen will be the grouping of the values. For example, let's define a dose of amount 15 at times t=0,4,8, and 12:

regimen = DosageRegimen([15,15,15,15], time=[0,4,8,12])
Pumas.DosageRegimen(4×8 DataFrames.DataFrame
│ Row │ time    │ cmt   │ amt     │ evid │ ii      │ addl  │ rate    │ ss  
 │
│     │ Float64 │ Int64 │ Float64 │ Int8 │ Float64 │ Int64 │ Float64 │ Int8
 │
├─────┼─────────┼───────┼─────────┼──────┼─────────┼───────┼─────────┼─────
─┤
│ 1   │ 0.0     │ 1     │ 15.0    │ 1    │ 0.0     │ 0     │ 0.0     │ 0   
 │
│ 2   │ 4.0     │ 1     │ 15.0    │ 1    │ 0.0     │ 0     │ 0.0     │ 0   
 │
│ 3   │ 8.0     │ 1     │ 15.0    │ 1    │ 0.0     │ 0     │ 0.0     │ 0   
 │
│ 4   │ 12.0    │ 1     │ 15.0    │ 1    │ 0.0     │ 0     │ 0.0     │ 0   
 │)

Let's define our subject to have id=1 and this multiple dosing regimen:

subject = Subject(id=1,evs=regimen)
Subject
  ID: 1
  Events: 4

Running a Simulation

The main function for running a simulation is simobs. simobs on a population simulates all of the population (in parallel), while simobs on a subject simulates just that subject. If we wish to change the parameters from the initialized values, then we pass them in. Let's simulate subject 1 with a set of chosen parameters:

fixeffs = (θ = [
          1, # Ka1  Absorption rate constant 1 (1/time)
          1, # CL   Clearance (volume/time)
          20, # Vc   Central volume (volume)
          2, # Q    Inter-compartmental clearance (volume/time)
          10, # Vp   Peripheral volume of distribution (volume)
          10, # Kin  Response in rate constant (1/time)
          2, # Kout Response out rate constant (1/time)
          2, # IC50 Concentration for 50% of max inhibition (mass/volume)
          1, # IMAX Maximum inhibition
          1, # γ    Emax model sigmoidicity
          0, # Vmax Maximum reaction velocity (mass/time)
          2  # Km   Michaelis constant (mass/volume)
          ],)

sim = simobs(model, subject, fixeffs)
Pumas.SimulatedObservations{Pumas.Subject{Nothing,Nothing,Array{Pumas.Event
,1},Nothing},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePr
ecision{Float64}},NamedTuple{(:ev1, :cp, :periph, :resp),NTuple{4,Array{Flo
at64,1}}}}(Subject
  ID: 1
  Events: 4
, 0.0:1.0:36.0, (ev1 = [15.0, 5.51819, 2.03003, 0.746809, 15.2747, 5.61926,
 2.06718, 0.760473, 15.2798, 5.62111  …  4.68203e-6, 1.72262e-6, 6.35957e-7
, 2.09712e-7, 6.93428e-8, 2.26354e-8, 7.32975e-9, 2.52184e-9, 8.61695e-10, 
2.87507e-10], cp = [0.0, 0.310673, 0.378429, 0.385297, 0.371218, 0.661221, 
0.706952, 0.692245, 0.657675, 0.928418  …  0.485284, 0.452455, 0.421846, 0.
393506, 0.366909, 0.342143, 0.319052, 0.297519, 0.277441, 0.258717], periph
 = [0.0, 2.82347, 4.17095, 4.45998, 4.36897, 6.97508, 8.07065, 8.10708, 7.7
7384, 10.1519  …  5.76705, 5.37969, 5.01839, 4.67705, 4.36222, 4.06782, 3.7
9317, 3.53721, 3.29851, 3.07588], resp = [5.0, 0.726534, -2.38598, -3.76747
, -4.23693, -5.60638, -6.65935, -7.07949, -7.1205, -7.82036  …  -5.97321, -
5.65524, -5.33876, -5.02546, -4.71238, -4.40052, -4.09085, -3.78344, -3.478
44, -3.17551]))

We can then plot the simulated observations by using the plot command:

using Plots
plot(sim)

Note that we can use the attributes from Plots.jl to further modify the plot. For example,

plot(sim,
     color=2,thickness_scaling=1.5,
     legend=false, lw=2)

Notice that in our model we said that there was a single parameter θ so our input parameter is a named tuple with just the name θ. When we only give the parameters, the random effects are automatically sampled from their distributions. If we wish to prescribe a value for the random effects, we pass initial values similarly:

randeffs = (η = rand(11),)
sim = simobs(model, subject, fixeffs, randeffs)
plot(sim)

The points which are saved are by default at once every hour until one day after the last event. If you wish to change the saving time points, pass the keyword argument obstimes. For example, let's save at every 0.1 hours and run the simulation for 19 hours:

sim = simobs(model, subject, fixeffs, randeffs, obstimes = 0:0.1:19)
plot(sim)

Handling the SimulatedObservations

The resulting SimulatedObservations type has two fields. sim.times is an array of time points for which the data was saved. sim.derived is the result of the derived variables. From there, the derived variabes are accessed by name. For example,

sim[:cp]
191-element Array{Float64,1}:
 0.0                
 0.07056857729487821
 0.13288275278255193
 0.1877877868729826 
 0.23604600778739995
 0.2783447990408738 
 0.3153039218915116 
 0.3474818003312944 
 0.3753819050898568 
 0.3994573778583226 
 ⋮                  
 1.1275116014140318 
 1.1218178171278705 
 1.1162123942312068 
 1.1106930077397865 
 1.1052574096171284 
 1.09990338452835   
 1.0946287488974582 
 1.089431350907351  
 1.084309070499816

is the array of cp values at the associated time points. We can turn this into a DataFrame via using the DataFrame command:

df = DataFrame(sim)
first(df,6) # Print only the first 6: the DataFrame is huge!

6 rows × 9 columns

timeev1cpperiphrespamtevidcmtrate
Float64Float64Float64Float64Float64Float64Int8Int64⍰Float64
10.015.00.00.05.015.0110.0
20.015.00.00.05.00.00missing0.0
30.113.57260.07056860.01168824.469150.00missing0.0
40.212.2810.1328830.0445554.118530.00missing0.0
50.311.11230.1877880.09557373.88640.00missing0.0
60.410.05480.2360460.1620483.732130.00missing0.0

From there, any Julia tools can be used to analyze these arrays and DataFrames. For example, if we wish the plot the result of ev1 over time, we'd use the following:

plot(sim.times,sim[:ev1])

Using these commands, a Julia program can be written to post-process the program however you like!

Conclusion

This tutorial covered basic workflow for how to build a model and simulate results from it. The subsequent models will go into more detail in the components, such as:

  1. More detailed treatment of specifying populations, dosage regimens, and covariates.

  2. Reading in dosage regimens and observations from NMTRAN data.

using PumasTutorials
PumasTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

Appendix

These tutorials are part of the PumasTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

To locally run this tutorial, do the following commands:

using PumasTutorials
PumasTutorials.weave_file("introduction","introduction.jmd")

Computer Information:

Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-8550U CPU @ 1.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = "C:\Users\accou\AppData\Local\atom\app-1.38.2\atom.exe"  -a
  JULIA_NUM_THREADS = 4

Package Information:

Status `C:\Users\accou\.julia\environments\v1.1\Project.toml`
[621f4979-c628-5d54-868e-fcf4e3e8185c] AbstractFFTs 0.4.1
[c52e3926-4ff0-5f6e-af25-54175e0327b1] Atom 0.8.8
[f0abef60-9ec0-11e9-27de-db6506a91768] AutoOffload 0.1.0
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.4.2
[4ece37e6-a012-11e8-38cd-91247efc2c34] Bioequivalence 0.1.0
[336ed68f-0bac-5ca0-87d4-7b16caf5d00b] CSV 0.5.9
[c5f51814-7f29-56b8-a69c-e4d8f6be1fde] CUDAdrv 3.0.1
[be33ccc6-a3ff-5ff2-a52e-74243cff1e17] CUDAnative 2.2.1
[49dc2e85-a5d0-5ad3-a950-438e2897f1b9] Calculus 0.5.0
[7057c7e9-c182-5462-911a-8362d720325c] Cassette 0.2.5
[34da2185-b29b-5c13-b0c7-acf172513d20] Compat 2.1.0
[3a865a2d-5b23-5a0f-bc46-62713ec82fae] CuArrays 1.1.0
[667455a9-e2ce-5579-9412-b964f529a492] Cubature 1.4.0
[a93c6f00-e57d-5684-b7b6-d8193f3e46c0] DataFrames 0.18.4
[82cc6244-b520-54b8-b5a6-8a565e85f1d0] DataInterpolations 0.2.0
[31a5f54b-26ea-5ae9-a837-f05ce5417438] Debugger 0.5.0
[bcd4f6db-9728-5f36-b5f7-82caef46ccdb] DelayDiffEq 5.9.1
[2b5f629d-d688-5b77-993f-72d75c75574e] DiffEqBase 5.16.3
[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 1.2.0
[31c91b34-3c75-11e9-0341-95557aab0344] DiffEqBenchmarks 0.1.0
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.5.2+
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.13.0
[01453d9d-ee7c-5054-8395-0335cb756afa] DiffEqDiffTools 0.14.0
[aae7a2af-3d4f-5e19-a356-7da93b79d9d0] DiffEqFlux 0.5.2
[071ae1c0-96b5-11e9-1965-c90190d839ea] DiffEqGPU 0.1.0
[c894b116-72e5-5b58-be3c-e6d8d4ac2b12] DiffEqJump 6.1.1+
[8f2b45d5-b17b-5532-9e92-98ae0077e2e3] DiffEqMachineLearning 0.1.0
[78ddff82-25fc-5f2b-89aa-309469cbf16f] DiffEqMonteCarlo 0.15.1
[77a26b50-5914-5dd7-bc55-306e6241c503] DiffEqNoiseProcess 3.3.1
[9fdde737-9c7f-55bf-ade8-46b3f136cc48] DiffEqOperators 3.5.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.2.0
[a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9] DiffEqProblemLibrary 4.3.0
[41bf760c-e81c-5289-8e54-58b1f1f8abe2] DiffEqSensitivity 3.3.0
[6d1b261a-3be8-11e9-3f2f-0b112a9a8436] DiffEqTutorials 0.1.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.6.0
[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.20.0
[e30172f5-a6a5-5a46-863b-614d45cd2de4] Documenter 0.23.0
[587475ba-b771-5e3f-ad9e-33799f191a9c] Flux 0.8.3
[f6369f11-7733-5829-9624-2563aa707210] ForwardDiff 0.10.3+
[ba82f77b-6841-5d2e-bd9f-4daf811aec27] GPUifyLoops 0.2.5
[c91e804a-d5a3-530f-b6f0-dfbca275c004] Gadfly 1.1.0
[bc5e4493-9b4d-5f90-b8aa-2b2bcaad7a26] GitHub 5.1.1
[7073ff75-c697-5162-941a-fcdaad2a7d2a] IJulia 1.18.1
[42fd0dbc-a981-5370-80f2-aaf504508153] IterativeSolvers 0.8.1
[033835bb-8acc-5ee8-8aae-3f567f8a3819] JLD2 0.1.2
[e5e0dc1b-0480-54bc-9374-aad01c23163d] Juno 0.7.0
[2d691ee1-e668-5016-a719-b2531b85e0f5] LIBLINEAR 0.5.1
[7f56f5a3-f504-529b-bc02-0b1fe5e64312] LSODA 0.4.0
[6f1fad26-d15e-5dc8-ae53-837a1d7b8c9f] Libtask 0.3.0
[c7f686f2-ff18-58e9-bc7b-31028e88f75d] MCMCChains 0.3.10
[33e6dc65-8f57-5167-99aa-e5a354878fb2] MKL 0.0.0
[961ee093-0014-501f-94e3-6117800e7a78] ModelingToolkit 0.5.0
[4886b29c-78c9-11e9-0a6e-41e1f4161f7b] MonteCarloIntegration 0.0.1
[2774e3e8-f4cf-5e23-947b-6d7e65073b56] NLsolve 4.0.0
[8faf48c0-8b73-11e9-0e63-2155955bfa4d] NeuralNetDiffEq 0.1.0
[09606e27-ecf5-54fc-bb29-004bd9f985bf] ODEInterfaceDiffEq 3.3.1
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.12.0
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.2.0
[14b8a8f1-9102-5b29-a752-f990bacb7fe1] PkgTemplates 0.6.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.25.3
[92933f4c-e287-5a05-a399-4b506db050ca] ProgressMeter 1.0.0
[d7b8c89e-ad89-52e0-b9fd-d0ed321fa021] Pumas 0.1.0
[b7b41870-aa11-11e9-048a-09266ec4a62f] PumasTutorials 0.0.1
[438e738f-606a-5dbb-bf0a-cddfbfd45ab0] PyCall 1.91.2
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.8.1
[1fd47b50-473d-5c70-9696-f719f8f3bcdc] QuadGK 2.1.0
[612083be-0b0f-5412-89c1-4e7c75506a58] Queryverse 0.3.1
[6f49c342-dc21-5d91-9882-a32aef131414] RCall 0.13.3
[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 0.20.0
[37e2e3b7-166d-5795-8a7a-e32c996b4267] ReverseDiff 0.3.1
[295af30f-e4ad-537b-8983-00126c2a3abe] Revise 2.1.6
[2b6d1eac-7baa-5078-8adc-e6a3e659f14f] SingleFloats 0.1.3
[47a9eef4-7e08-11e9-0b38-333d64bd3804] SparseDiffTools 0.4.0
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.11.0
[4c63d2b9-4356-54db-8cca-17b64c39e42c] StatsFuns 0.8.0
[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.11.0
[9672c7b4-1e72-59bd-8a11-6ac3964bc41f] SteadyStateDiffEq 1.5.0
[789caeaf-c7a9-5a7d-9973-96adeb23e2a0] StochasticDiffEq 6.6.0
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 3.6.1
[fd094767-a336-5f1f-9728-57cf17d0bbfb] Suppressor 0.1.1
[6fc51010-71bc-11e9-0e15-a3fcc6593c49] Surrogates 0.1.0
[9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c] Tracker 0.2.2
[fce5fe82-541a-59a6-adf8-730c64b5f9a0] Turing 0.6.18
[1986cc42-f94f-5a68-af5c-568840ba703d] Unitful 0.16.0
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.9.1
[e88e6eb3-aa80-5325-afca-941959d7151f] Zygote 0.3.2