using Dates
using Random
using Pumas
using PumasUtilities
Random.seed!(1234)
Generating and Simulating Populations
1 Introduction
In this tutorial, we will cover the fundamentals of generating populations to simulate with Pumas. We will demonstrate how to specify dosage regimens and covariates, and then how to piece these together to form a population to simulate.
2 The model
Below is a Pumas model that specifies a 1-compartment oral absorption system with between-subject variability on all the parameters. Details of the model specification are provided in the documentation.
= @model begin
model @param begin
∈ RealDomain(; lower = 0)
tvcl ∈ RealDomain(; lower = 0)
tvvc ∈ RealDomain(; lower = 0)
tvka ∈ PSDDomain(3)
Ω ∈ RealDomain(; lower = 0)
σ_prop end
@random begin
~ MvNormal(Ω)
η end
@covariates Wt
@pre begin
= tvcl * (Wt / 70)^0.75 * exp(η[1])
CL = tvvc * (Wt / 70) * exp(η[2])
Vc = tvka * exp(η[3])
Ka end
@dynamics begin
' = -Ka * Depot
Depot' = Ka * Depot - Central * CL / Vc
Centralend
@derived begin
= @. Central / Vc
conc ~ @. Normal(conc, conc * σ_prop)
dv end
end
PumasModel
Parameters: tvcl, tvvc, tvka, Ω, σ_prop
Random effects: η
Covariates: Wt
Dynamical system variables: Depot, Central
Dynamical system type: Matrix exponential
Derived: conc, dv
Observed: conc, dv
3 Setting up parameters
Next we provide the initial estimates of the parameters to simulate from. The fixed effects are provided in the typical values tv
parameters (tvcl
, tvvc
, tvka
) and the between-subject variability parameters are provided in the Ω matrix as a covariance matrix. So, 0.04 variance on Ω₁,₁
suggests a 20% coefficient of variation (\(\sqrt{0.04} = 0.2\)). Similarly, σ_prop
has a 20% proportional residual error.
=
fixeffs = 0.4, tvvc = 20, tvka = 1.1, Ω = Diagonal([0.04, 0.04, 0.04]), σ_prop = 0.2) (; tvcl
(tvcl = 0.4,
tvvc = 20,
tvka = 1.1,
Ω = [0.04 0.0 0.0; 0.0 0.04 0.0; 0.0 0.0 0.04],
σ_prop = 0.2,)
4 Single dose example
DosageRegimen()
is the function that lets you construct a dosing regimen. The first argument of the DosageRegimen
is amt
and is a positional argument and not a keyword argument. All subsequent arguments need to be named, since they are keyword arguments. Let’s try a simple example where you provide a 100 mg dose at time = 0
.
= DosageRegimen(100; time = 0) ev
Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
---|---|---|---|---|---|---|---|---|---|---|
Float64 | Int64 | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
1 | 0.0 | 1 | 100.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
As you can see above, we provided a single 100 mg dose. DosageRegimen
provides some defaults when it creates the dosage regimen dataset, time = 0
, evid = 1
, cmt = 1
, rate = 0
, ii = 0
and addl = 0
.
Note that ev
is of type DosageRegimen
. Specified like above, DosageRegimen
is one of the four fundamental building block of a Subject
(more on Subject
below).
4.1 Building Subjects
Let’s create a single subject
= Subject(;
s1 = 1,
id = ev,
events = (; Wt = 70),
covariates = (; dv = nothing),
observations )
Subject
ID: 1
Events: 1
Observations: dv: (n=0)
Covariates: Wt
Note that the s1
Subject
created above is composed of:
id
: an unique identifierobservations
: observations, represented by aNamedTuple
covariates
: covariates, represented by aNamedTuple
events
: events, represented by aDosageRegimen
DataFrame(s1.events)
Row | time | amt | evid | cmt | rate | duration | ss | ii | route |
---|---|---|---|---|---|---|---|---|---|
Float64 | Float64 | Int8 | Int64 | Float64 | Float64 | Int8 | Float64 | NCA.Route | |
1 | 0.0 | 100.0 | 1 | 1 | 0.0 | 0.0 | 0 | 0.0 | NullRoute |
The events are presented by basic information such as the dose of drug, the time of dose administration, the compartment number for administration, and whether the dose is an instantaneous input or an infusion.
Below is how the covariates are represented, as a NamedTuple
.
s1.covariates
Pumas.ConstantCovar{NamedTuple{(:Wt,), Tuple{Int64}}}((Wt = 70,))
(Note: defining distributions for covariates will be discussed in detail later.)
Using this one subject, s1
, let us simulate a simple concentration time profile using the model above, and plot the result:
= simobs(model, s1, fixeffs; obstimes = 0:0.1:120)
sims1
sim_plot(
model,
sims1;= [:conc],
observations = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "Predicted Concentration (ng/mL)"),
axis )
4.2 Building Populations
Now, lets create one more subject, s2
.
= Subject(;
s2 = 2,
id = ev,
events = (; Wt = 70),
covariates = (; dv = nothing),
observations )
Subject
ID: 2
Events: 1
Observations: dv: (n=0)
Covariates: Wt
If we want to simulate both s1
and s2
together, we need to bring these subjects together to form a Population
. A Population
is essentially a collection (vector) of subjects.
= [s1, s2] twosubjs
Population
Subjects: 2
Covariates: Wt
Observations: dv
Let’s see the details of the first and the second subject:
1] twosubjs[
Subject
ID: 1
Events: 1
Observations: dv: (n=0)
Covariates: Wt
2] twosubjs[
Subject
ID: 2
Events: 1
Observations: dv: (n=0)
Covariates: Wt
Now, we can simulate this Population
of 2 subjects, and plot the result:
= simobs(model, twosubjs, fixeffs; obstimes = 0:0.1:120)
sims2
sim_plot(
model,
sims2;= [:conc],
observations = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "Predicted Concentration (ng/mL)"),
axis )
Similarly, we can build a population of any number of subjects. But before we do that, let’s dive into covariate generation.
4.3 Covariates
As was discussed earlier, a Subject
can also be provided details regarding covariates. In the model above, a continuous covariate body weight Wt
impacts both CL
and Vc
. Let us now specify covariates to a population of 10 subjects:
choose_covariates() = (; Wt = rand(55:80))
choose_covariates (generic function with 1 method)
choose_covariates
will randomly choose a Wt
between 55-80 kgs
We can make a vector with covariates for ten subjects through an array comprehension:
= [choose_covariates() for i = 1:10] cvs
Now, we add these covariates to the population as below. The map(f, xs)
will return the result of f
on each element of the collection xs
. Let’s map a function that build’s a subject with the randomly chosen covariates in order to build a Population
:
= map(
pop_with_covariates -> Subject(;
i = i,
id = ev,
events = choose_covariates(),
covariates = (; dv = nothing),
observations
),1:10,
)
Population
Subjects: 10
Covariates: Wt
Observations: dv
Simulate into the population, and visualize the output:
= simobs(model, pop_with_covariates, fixeffs; obstimes = 0:0.1:120)
sims3
sim_plot(
model,
sims3;= [:conc],
observations = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "Predicted Concentration (ng/mL)"),
axis )
5 Multiple dose example
DosageRegimen
provides additional controls over the specification of the dosage regimen. For example, ii
defines the “interdose interval”, or the time difference between two doses, while addl
defines how many additional times to repeat a dose. Thus, let’s define a dose of 100
that’s repeated 7
times at 24
hour intervals:
= DosageRegimen(100; ii = 24, addl = 6) md
Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
---|---|---|---|---|---|---|---|---|---|---|
Float64 | Int64 | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
1 | 0.0 | 1 | 100.0 | 1 | 24.0 | 6 | 0.0 | 0.0 | 0 | NullRoute |
Let’s create a new subject, s3
with this dosage regimen, and see the results:
= Subject(;
s3 = 3,
id = md,
events = (; Wt = 70),
covariates = (; dv = nothing),
observations
)
= simobs(model, s3, fixeffs; obstimes = 0:0.1:240)
sims4
sim_plot(
model,
[sims4];= [:conc],
observations = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "Predicted Concentration (ng/mL)"),
axis )
6 Combining dosage regimens
We can also combine dosage regimens to build a more complex regimen. Recall from the introduction that using arrays will build the element-wise combinations. Thus, let’s build a dose of 500
into compartment 1
at time 0
, and 7
doses into compartment 1
of 100
spaced by 24
hours:
= DosageRegimen([500, 100]; time = [0, 24], cmt = 1, addl = [0, 6], ii = [0, 24]) ldmd
Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
---|---|---|---|---|---|---|---|---|---|---|
Float64 | Int64 | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
1 | 0.0 | 1 | 500.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
2 | 24.0 | 1 | 100.0 | 1 | 24.0 | 6 | 0.0 | 0.0 | 0 | NullRoute |
Let’s see if this result matches our intuition:
= Subject(;
s4 = 4,
id = ldmd,
events = (; Wt = 70),
covariates = (; dv = nothing),
observations
)
= simobs(model, s4, fixeffs; obstimes = 0:0.1:120)
sims5
sim_plot(
model,
[sims5];= [:conc],
observations = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "Predicted Concentration (ng/mL)"),
axis )
Another way to build complex dosage regimens is to combine previously constructed regimens into a single regimen. For example:
= DosageRegimen(500; cmt = 1, time = 0, addl = 0, ii = 0)
e1 = DosageRegimen(100; cmt = 1, time = 24, addl = 6, ii = 24)
e2 = DosageRegimen(e1, e2) evs
Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
---|---|---|---|---|---|---|---|---|---|---|
Float64 | Int64 | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
1 | 0.0 | 1 | 500.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
2 | 24.0 | 1 | 100.0 | 1 | 24.0 | 6 | 0.0 | 0.0 | 0 | NullRoute |
= Subject(;
s5 = 5,
id = evs,
events = (; Wt = 70),
covariates = (; dv = nothing),
observations
)
= simobs(model, s5, fixeffs; obstimes = 0:0.1:120)
sims6
sim_plot(
model,
[sims6];= [:conc],
observations = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "Predicted Concentration (ng/mL)"),
axis )
is the same regimen as before.
Putting these ideas together, we can define a Population
where individuals with different covariates undergo different regimens, and simulate them all together:
= DosageRegimen(100; ii = 24, addl = 6)
e3 = DosageRegimen(50; ii = 12, addl = 13)
e4 = DosageRegimen(200; ii = 24, addl = 2)
e5
= map(
pop1 -> Subject(;
i = i,
id = e3,
events = choose_covariates(),
covariates = (; dv = nothing),
observations
),1:5,
)
= map(
pop2 -> Subject(;
i = i,
id = e4,
events = choose_covariates(),
covariates = (; dv = nothing),
observations
),6:8,
)
= map(
pop3 -> Subject(;
i = i,
id = e5,
events = choose_covariates(),
covariates = (; dv = nothing),
observations
),9:10,
)
= vcat(pop1, pop2, pop3)
pop
= simobs(model, pop, fixeffs; obstimes = 0:0.1:120)
sims7
sim_plot(
model,
sims7;= [:conc],
observations = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "Predicted Concentration (ng/mL)"),
axis )
7 Defining Infusions
An infusion is a dosage which is defined as having a non-zero positive rate at which the drug enters the system. Let’s define a single infusion dose of total amount 100
with a rate of 3
unit amount per time into the compartment 2
:
= DosageRegimen(100; rate = 3, cmt = 2) inf
Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
---|---|---|---|---|---|---|---|---|---|---|
Float64 | Int64 | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
1 | 0.0 | 2 | 100.0 | 1 | 0.0 | 0 | 3.0 | 33.3333 | 0 | NullRoute |
Now let’s simulate a subject undergoing this treatment strategy:
= Subject(;
s6 = 5,
id = inf,
events = (; Wt = 70),
covariates = (; dv = nothing),
observations
)
= simobs(model, s6, fixeffs; obstimes = 0:0.1:120)
sims8
sim_plot(
model,
[sims8];= [:conc],
observations = (; fontsize = 18),
figure = (; xlabel = "Time (hr)", ylabel = "Predicted Concentration (ng/mL)"),
axis )
8 Final Note on Julia Programming
Note that all of these functions are standard Julia functions, and thus standard Julia programming constructions can be utilized to simplify the construction of large populations. We already demonstrated the use of map
and array comprehension, but we can also make use of constructs like for
loops and if-else
statements.
9 Conclusion
This tutorial shows the tools for generating populations of infinite complexity, defining covariates and dosage regimens on the fly and simulating the results of the model.