# Introduction to Pumas

This is an introduction to Pumas, a software for pharmaceutical 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 single doses, and analyze the results of the simulation. Fitting of data using any of the methods available in Pumas is showcased in the tutorials. This tutorial is designed as a broad overview of the workflow and a more in-depth treatment of each section can be found in the tutorials.

## Working Example

Let's start by showing a complete simulation code, and then break down how it works.

using Random
using Pumas
using PumasUtilities
using CairoMakie
######### Turnover model

inf_2cmt_lin_turnover = @model begin
@param   begin
tvcl ∈ RealDomain(lower=0)
tvvc ∈ RealDomain(lower=0)
tvq ∈ RealDomain(lower=0)
tvvp ∈ RealDomain(lower=0)
Ω_pk ∈ PDiagDomain(4)
σ_prop_pk ∈ RealDomain(lower=0)
# PD parameters
tvturn ∈ RealDomain(lower=0)
tvebase ∈ RealDomain(lower=0)
tvec50 ∈ RealDomain(lower=0)
Ω_pd ∈ PDiagDomain(1)
end

@random begin
ηpk ~ MvNormal(Ω_pk)
ηpd ~ MvNormal(Ω_pd)
end

@pre begin
CL    = tvcl    * exp(ηpk[1])
Vc    = tvvc    * exp(ηpk[2])
Q     = tvq     * exp(ηpk[3])
Vp    = tvvp    * exp(ηpk[4])

ebase = tvebase * exp(ηpd[1])
ec50  = tvec50
emax  = 1
turn  = tvturn
kout  = 1/turn
kin0  = ebase   * kout
end

@init begin
Resp = ebase
end

@vars begin
conc  := Central/Vc
edrug := emax*conc/(ec50 + conc)
kin   := kin0*(1-edrug)
end

@dynamics begin
Central'    =  - (CL/Vc)*Central + (Q/Vp)*Peripheral - (Q/Vc)*Central
Peripheral' = (Q/Vc)*Central - (Q/Vp)*Peripheral
Resp'       = kin - kout*Resp
end

@derived begin
dv   ~ @. Normal(conc, conc*σ_prop_pk)
end
end

turnover_params = (tvcl = 1.5,
tvvc = 25.0,
tvq = 5.0,
tvvp = 150.0,
tvturn = 10,
tvebase = 10,
tvec50 = 0.3,
Ω_pk = Diagonal([0.05,0.05,0.05,0.05]),
Ω_pd = Diagonal([0.05]),
σ_prop_pk = 0.15,
σ_add_pd = 0.2)
regimen = DosageRegimen(150, rate = 10, cmt = 1)
pop     = map(i -> Subject(id = i,events = regimen), 1:10)
sd_obstimes = [0.25, 0.5, 0.75, 1, 2, 4, 8,
12, 16, 20, 24, 36, 48, 60, 71.9] # single dose observation times
Random.seed!(123)
sims = simobs(inf_2cmt_lin_turnover,
pop,
turnover_params,
obstimes = sd_obstimes)
sim_plot(inf_2cmt_lin_turnover,
sims, observations =[:dv],
figure = (fontsize = 18, ),
axis = (xlabel = "Time (hr)",
ylabel = "Concentration (ng/mL)", ))

sim_plot(inf_2cmt_lin_turnover,
sims, observations =[:resp],
figure = (fontsize = 18, ),
axis = (xlabel = "Time (hr)",
ylabel = "Response", ))

In this code, we defined a nonlinear mixed effects model by describing the parameters, the random effects, the dynamical model, and the derived (result) values. Then we generated a population of 10 subjects who received a single dose of 150mg, specified parameter values, simulated the model, and generated a plot of the results. Now let's walk through this process!

## Using the Model Macro

First we define the model. The simplest way to do is via the @model DSL. 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 structural model parameters of PK and PD and their corresponding variances where applicable:

@param   begin
tvcl ∈ RealDomain(lower=0)
tvvc ∈ RealDomain(lower=0)
tvq ∈ RealDomain(lower=0)
tvvp ∈ RealDomain(lower=0)
Ω_pk ∈ PDiagDomain(4)
σ_prop_pk ∈ RealDomain(lower=0)
# PD parameters
tvturn ∈ RealDomain(lower=0)
tvebase ∈ RealDomain(lower=0)
tvec50 ∈ RealDomain(lower=0)
Ω_pd ∈ PDiagDomain(1)
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 uncorrelated random effects, one for PK and another for PD so we utilize the syntax:

ηpk ~ MvNormal(Ω_pk)
ηpd ~ MvNormal(Ω_pd)

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
CL = tvcl * exp(ηpk[1])
Vc = tvvc * exp(ηpk[2])
Q = tvq * exp(ηpk[3])
Vp = tvvp * exp(ηpk[4])

ebase = tvebase*exp(ηpd[1])
ec50 = tvec50
emax = 1
turn = tvturn
kout = 1/turn
kin0 = ebase*kout
end

Next we define the @init block which gives the initial 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 = ebase
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
Central' =  - (CL/Vc)*Central + (Q/Vp)*Peripheral - (Q/Vc)*Central
Peripheral' = (Q/Vc)*Central - (Q/Vp)*Peripheral
Resp' = kin - kout*Resp
end

Next we setup alias variables that can be used later in the code. Such alias code can be setup in the @vars block

@vars begin
conc := Central/Vc
edrug := emax*conc/(ec50 + conc)
kin := kin0*(1-edrug)
end

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

@derived begin
dv ~ @. Normal(conc, sqrt(conc^2*σ_prop_pk))
end

## Building a Subject

Now let's build a subject to simulate the model with. A subject is defined by the following components:

1. An identifier - id
2. The dosage regimen - events
3. The covariates of the individual - covariates
4. Observations associated with the individual - observations
5. The timing of the sampling - time
6. A vector of times if the covariates are time-varying - covariate_time
7. Interpolation direction for covariates - covariate_direction

Our model did not make use of covariates so we will ignore (3, 6 and 7) for now, and (4) 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. 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(150, time=0)

is a dosage regimen which simply does a single dose at time t=0 of amount 150. Let's assume the dose is an infusion administered at the rate of 10 mg/hour into the first compartment

regimen = DosageRegimen(150, rate=10, cmt=1)

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

subject = Subject(id = 1, events = regimen)

You can also create a collection of subjects which becomes a Population.

pop = Population(map(i -> Subject(id= i, events = regimen), 1:10))

## Running a Simulation

The main function for running a simulation is simobs. simobs on a Population simulates all of the population, 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:

turnover_params = (tvcl = 1.5,
tvvc = 25.0,
tvq = 5.0,
tvvp = 150.0,
tvturn = 10,
tvebase = 10,
tvec50 = 0.3,
Ω_pk = Diagonal([0.05,0.05,0.05,0.05]),
Ω_pd = Diagonal([0.05]),
σ_prop_pk = 0.15,

Random.seed!(123)
sims = simobs(inf_2cmt_lin_turnover,
pop,
turnover_params,
obstimes = sd_obstimes)

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

using PlottingUtilities
sim_plot(inf_2cmt_lin_turnover, sims, observations = :dv)

Note that we can use the Customization from PumasUtilities to further modify the plot. For example,

sim_plot(inf_2cmt_lin_turnover,
sims, observations =[:dv],
figure = (fontsize = 18, ),
axis = (xlabel = "Time (hr)",
ylabel = "Concentration (ng/mL)", ))

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 = (ηpk = rand(4), ηpd = rand(1),)

If a population simulation is required with no random effects, then the values of the η's can be set to zero that will result in a simulation only at the mean level:

etas = zero_randeffs(inf_2cmt_lin_turnover,
turnover_params,
pop)
Random.seed!(123)
sims = simobs(inf_2cmt_lin_turnover,
pop,
turnover_params,
etas;
obstimes = sd_obstimes)
fig =Figure(resolution=(1000,600))
f1 = sim_plot(fig[1,1], inf_2cmt_lin_turnover,
sims, observations =[:dv],
figure = (fontsize = 18, ),
axis = (xlabel = "Time (hr)",
ylabel = "Concentration (ng/mL)", ))
f2 = sim_plot(fig[1,2], inf_2cmt_lin_turnover,
sims, observations =[:resp],
figure = (fontsize = 18, ),
axis = (xlabel = "Time (hr)",
ylabel = "Response", ))
fig                           

You still see variability in the plot above, mainly due to the residual variability components in the model. It is quite trivial to change the parameter estimates of only a subset of parameters as below

turnover_params_wo_sigma  = (turnover_params..., σ_prop_pk = 0.0, σ_add_pd = 0.0)

and now if we perform the simulation again without random effects to get only the population mean,

etas = zero_randeffs(inf_2cmt_lin_turnover,
turnover_params,
pop)
Random.seed!(123)
sims = simobs(inf_2cmt_lin_turnover,
pop,
turnover_params_wo_sigma,
etas;
obstimes = sd_obstimes)
fig =Figure(resolution=(1000,600))
f1 = sim_plot(fig[1,1], inf_2cmt_lin_turnover,
sims, observations =[:dv],
figure = (fontsize = 18, ),
axis = (xlabel = "Time (hr)",
ylabel = "Concentration (ng/mL)", ))
f2 = sim_plot(fig[1,2], inf_2cmt_lin_turnover,
sims, observations =[:resp],
figure = (fontsize = 18, ),
axis = (xlabel = "Time (hr)",
ylabel = "Response", ))
fig 

## Handling the SimulatedObservations

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

sims[1].observations[:dv]

is the array of dv values at the associated time points for subject 1. We can convert this to a normal Pumas Subjects using the constructor:

julia> sims_subjects = Subject.(sims)
Population
Subjects: 10
Observations: dv, resp

and sims_subjects can then be used fit the model parameters based on a simulated dataset. We can also turn the simulated population into a DataFrame via using the DataFrame command:

julia> DataFrame(sims)
160×8 DataFrame
Row │ id      time     dv              resp           amt        evid  cmt      rate
│ String  Float64  Float64?        Float64?       Float64?   Int8  Int64?   Float64
─────┼───────────────────────────────────────────────────────────────────────────────────
1 │ 1          0.0   missing         missing            150.0     1        1     10.0
2 │ 1          0.25        0.096826        9.96663  missing       0  missing      0.0
3 │ 1          0.5         0.187597        9.88839  missing       0  missing      0.0
4 │ 1          0.75        0.272731        9.78409  missing       0  missing      0.0
5 │ 1          1.0         0.352616        9.66347  missing       0  missing      0.0
6 │ 1          2.0         0.626531        9.10546  missing       0  missing      0.0
7 │ 1          4.0         1.01152         7.93649  missing       0  missing      0.0
8 │ 1          8.0         1.43067         5.95748  missing       0  missing      0.0
⋮  │   ⋮        ⋮           ⋮               ⋮            ⋮       ⋮       ⋮        ⋮
154 │ 10        16.0         1.47042         3.52807  missing       0  missing      0.0
155 │ 10        20.0         0.81701         3.10522  missing       0  missing      0.0
156 │ 10        24.0         0.599389        3.09228  missing       0  missing      0.0
157 │ 10        36.0         0.467029        3.55706  missing       0  missing      0.0
158 │ 10        48.0         0.426772        3.89715  missing       0  missing      0.0
159 │ 10        60.0         0.392367        4.14396  missing       0  missing      0.0
160 │ 10        71.9         0.361073        4.36054  missing       0  missing      0.0
145 rows omitted

From there, any Julia tools can be used to analyze these arrays and DataFrames.

## Conclusion

This tutorial covered basic workflow for how to build a model and simulate results from it. The subsequent tutorials 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 standard input data.
3. Fitting models