Dosage Regimens, Subjects, and Populations

In Pumas, subjects are represented by the Subject type and collections of subjects are represented as Vectors of Subjects and aliased as a Population. Subjects are defined by their identifier, observations, covariates, and events. In this section we will specify the methods used for defining Subjects programmatically or using the read_pumas function that reads in data that follows the Pumas NLME Data Format (PumasNDF) data format. Before we look at Subjects, we will take a look at how to define events as represented by the DosageRegimen type.

Dosage Regimen Terminology

When subjects are subjected to treatment it is represented by an event in Pumas. Administration of a drug is represented by a DosageRegimen that describes the amount, type, frequency and route. DosageRegimens can either be constructed programmatically using the DosageRegimen constructor or from a data source in the PumasNDF format using read_pumas. The names of the inputs are the same independent of how the DosageRegimen is constructed. The definition of the values are as follows:

• amt: the amount of the dose. This is the only required value.
• time: the time at which the dose is given. Defaults to 0.
• evid: the event id. 1 specifies a normal event. 3 means it's a reset event, meaning that the value of the dynamical variable is reset to the amt at the dosing event. If 4, then the dynamical value and time are reset, and then a final dose is given. Defaults to 1.
• ii: the inter-dose interval. For steady state events, this is the length of time between successive doses. When addl is specified, this is the length of time between dose events. Defaults to 0.
• addl: the number of additional events of the same type, spaced by ii. Defaults to 0.
• rate: the rate of administration. If 0, then the dose is instantaneous. Otherwise the dose is administrated at a constant rate for a duration equal to amt/rate.
• ss: an indicator for whether the dose is a steady state dose. A steady state dose is defined as the result of having applied the dose with the interval ii infinitely many successive times. 0 indicates that the dose is not a steady state dose. 1 indicates that the dose is a steady state dose. 2 indicates that it is a steady state dose that is added to the previous amount. The default is 0.
• route: route of administration to be used in NCA analysis if it is carried out with the integrated interface inside @model. Defaults to NullRoute which is basically no route specified.

This specification leads to the following default constructor for the DosageRegimen type:

DosageRegimen(
  amt::Numeric;
  time::Numeric = 0,
  cmt::Union{Numeric,Symbol} = 1,
  evid::Numeric = 1,
  ii::Numeric = zero.(time),
  addl::Numeric = 0,
  rate::Numeric = zero.(amt)./oneunit.(time),
  duration::Numeric = zero(amt)./oneunit.(time),
  ss::Numeric = 0,
  route::NCA.Route = NCA.NullRoute,
)

julia> DosageRegimen(100; ii = 24, addl = 6)
DosageRegimen
 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

julia> DosageRegimen(50; ii = 12, addl = 13)
DosageRegimen
 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     50.0     1     12.0     13      0.0       0.0     0  NullRoute

julia> DosageRegimen(200; ii = 24, addl = 2)
DosageRegimen
 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    200.0     1     24.0      2      0.0       0.0     0  NullRoute

julia> DosageRegimen(200; ii = 24, addl = 2, route = NCA.IVBolus)
DosageRegimen
 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    200.0     1     24.0      2      0.0       0.0     0  IVBolus

evs = DosageRegimen(
  regimen1::DosageRegimen,
  regimen2::DosageRegimen;
  offset = nothing
)

julia> e1 = DosageRegimen(100; ii = 24, addl = 6)
DosageRegimen
 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

julia> e2 = DosageRegimen(50; ii = 12, addl = 13)
DosageRegimen
 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     50.0     1     12.0     13      0.0       0.0     0  NullRoute

julia> evs = DosageRegimen(e1, e2)
DosageRegimen
 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
   2 │     0.0      1     50.0     1     12.0     13      0.0       0.0     0  NullRoute

julia> DosageRegimen(e1, e2; offset = 10)
DosageRegimen
 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
   2 │   178.0      1     50.0     1     12.0     13      0.0       0.0     0  NullRoute

Steady-State Dosing

In addition to applying the specified dose, steady-state dosing resets the compartments to a steady state (if ss = 1) or adds the the steady-state values to the current compartment values (if ss = 2).

There are three classes of steady-state doses:

The steady state of the system can be viewed as the result of having applied the steady-state dose from the infinite past, depending on the type of the dose as a constant infusion or as a sequence of doses. These so-called implied past doses do not affect the values of the system prior to the time point (time) of the dose record.

Constant Infusion

Steady-state events with constant infusion are characterized by

1. amt = 0: zero dose amount

1. rate > 0: positive infusion rate

The time point (time) of a steady-state dose record with constant infusion denotes the time point at which the infusion from the infinite past ends.

We illustrate this abstract concept with an example, using the following one-compartment model:

model = @model begin
    @pre begin
        Ka = 1.0
        CL = 1.0
        Vc = 10.0
    end

    @covariates BIOAV LAG

    @dosecontrol begin
        bioav = (; Central = BIOAV)
        lags = (; Central = LAG)
    end

    @dynamics Depots1Central1

    @observed begin
        conc = @. Central / Vc
    end
end

PumasModel
  Parameters: 
  Random effects: 
  Covariates: BIOAV, LAG
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: 
  Observed: conc

In this model, absorption rate and clearance are fixed to 1, volume of distribution is fixed to 10, and bioavailability and absorption lag time of the central compartment are set to covariates BIOAV and LAG, respectively. We study a subject

df = DataFrame(;
    id = 1,
    time = 0.0,
    amt = 0.0,
    cmt = :Depot,
    rate = 10.0,
    ii = 0.0,
    evid = 1,
    ss = 1,
    conc = missing,
    BIOAV = 0.5,
    LAG = 4.0,
)
sub = only(read_pumas(df; observations = [:conc], covariates = [:BIOAV, :LAG]))

Subject
  ID: 1
  Events: 1
  Observations: conc: (n=0)
  Covariates: BIOAV, LAG

with dosage regimen

sub.events

and covariates

sub.covariates(0.0)

(BIOAV = 0.5,
 LAG = 4.0,)

at the time point of the dose record. We see that this dose record is a steady-state event (evid = 1 and ss = 1) at time = 0.0 with constant infusion (amt = 0.0 and rate = 10.0 > 0) into the Depot compartment. Since ss = 1, the compartment amounts are reset to steady-state values at time = 0.0. The main question is: What are the steady-state values of the compartments at time = 0.0?

The steady state can be viewed as having applied an infusion of rate = 10.0 constantly from the infinite past up until the time point time = 0.0 of the steady-state event. We approximate the constant infusion with an infusion over a sufficiently long time span $[t_0, 0]$.

# Extract the compartment (`cmt`) of the steady-state dose
ss_dose = only(sub.events)
(; cmt) = ss_dose

# Start of the long infusion
t₀ = -500.0

sim_long_infusion = simobs(
    model,
    Subject(
        sub;
        covariates = sub.covariates(0.0),
        events = DosageRegimen(ss_dose.rate * (-t₀); time = t₀, cmt, ss_dose.rate),
    ),
    (;);
    obstimes = range(t₀, 0.0; length = 1_001),
)

SimulatedObservations
  Simulated variables: 
  Time: -500.0:0.5:0.0

The concentration and in particular also the compartment values converge to a steady state at time = 0.0. This visual observation can be confirmed by checking that the left-sided derivatives of the compartment values at time = 0.0 are zero:

Pumas.FiniteDiff.finite_difference_derivative(
    t -> sim_long_infusion.sol(-t),
    0.0,
    Val(:forward),
)

2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 0.0
 0.0

Thus, at time = 0.0, the time point of the steady-state event, the compartment values are set to these (approximations of the) steady-state values.

Depot_ss, Central_ss = sim_long_infusion.sol(0.0)

2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
  10.0
 100.0

model_approx = @model begin
    @pre begin
        Ka = 1.0
        CL = 1.0
        Vc = 10.0
    end

    @covariates BIOAV LAG

    @dosecontrol begin
        bioav = (; Central = BIOAV)
        lags = (; Central = LAG)
    end

    @init begin
        Depot = 10.0
        Central = 100.0
    end

    @dynamics Depots1Central1

    @observed begin
        conc = @. Central / Vc
    end
end

PumasModel
  Parameters: 
  Random effects: 
  Covariates: BIOAV, LAG
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: 
  Observed: conc

The only effect of a steady-state event with constant infusion is this change in the compartment values at the time point of the steady-state event. After this time point, simulation continues as usual.

# Default simulation
obstimes = range(0, 50; length = 1_001)
sim = simobs(model, sub, (;); obstimes)

# Simulation based on approximation of the steady state
sim_approx = simobs(
    model_approx,
    # Same subject, but without the steady-state event
    Subject(sub.id, sub.observations, sub.covariates, nothing, sub.time),
    (;);
    obstimes,
)

sim_approx.observed.conc ≈ sim.observed.conc

true

Multiple Bolus Doses

Steady-state events with multiple bolus doses are characterized by

1. amt > 0: positive dose amount

2. rate = 0: zero infusion rate

The time point (time) of a steady-state event with multiple bolus doses denotes the time point at which the bolus is injected.

Again we illustrate this type of steady-state events with an example, using the same one-compartment model as above:

model = @model begin
    @pre begin
        Ka = 1.0
        CL = 1.0
        Vc = 10.0
    end

    @covariates BIOAV LAG

    @dosecontrol begin
        bioav = (; Central = BIOAV)
        lags = (; Central = LAG)
    end

    @dynamics Depots1Central1

    @observed begin
        conc = @. Central / Vc
    end
end

PumasModel
  Parameters: 
  Random effects: 
  Covariates: BIOAV, LAG
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: 
  Observed: conc

Absorption rate and clearance are fixed to 1, volume of distribution is fixed to 10, and bioavailability and absorption lag time of the central compartment are set to covariates BIOAV and LAG, respectively. We study a subject

df = DataFrame(;
    id = 1,
    time = 0.0,
    amt = 10.0,
    cmt = :Central,
    rate = 0.0,
    ii = 5.0,
    addl = 2,
    evid = 1,
    ss = 1,
    conc = missing,
    BIOAV = 0.5,
    LAG = 4.0,
)
sub = only(read_pumas(df; observations = [:conc], covariates = [:BIOAV, :LAG]))

Subject
  ID: 1
  Events: 3
  Observations: conc: (n=0)
  Covariates: BIOAV, LAG

with dosage regimen

sub.events

and covariates

sub.covariates(0.0)

(BIOAV = 0.5,
 LAG = 4.0,)

at the time point of the dose record. We see that this dose record is a steady-state event (evid = 1 and ss = 1) at time = 0.0 with multiple bolus doses (amt = 10.0 > 0 and rate = 0.0) that are applied to the Central compartment every ii = 5.0 time units and repeated an additional addl = 2 times after the time point time = 0.0 of the dose record. Since ss = 1, the compartment amounts are reset to steady-state values at time = 0.0. The main question is: What are the steady-state values of the compartments at time = 0.0?

The steady state can be viewed as the compartment values resulting from having applied the sequence of implied doses from the infinite past up until the time point time = 0.0 of the steady-state event. We approximate this infinite sequence of doses with a sufficiently long sequence of doses.

# Extract dose amount (`amt`), interdose interval (`ii`), compartment (`cmt`),
# and time point (`time`) of the steady-state dose
ss_dose = only(sub.events)
(; amt, ii, cmt, time) = ss_dose

# Number of implied doses that are applied
n_implied = 40

sim_many_bolus = simobs(
    model,
    Subject(
        sub;
        covariates = sub.covariates(time),
        events = DosageRegimen(
            amt;
            time = time - n_implied * ii,
            ii,
            addl = n_implied - 1,
            cmt,
        ),
    ),
    (;);
    obstimes = range(time - n_implied * ii, time; length = 50 * n_implied + 1),
)

SimulatedObservations
  Simulated variables: 
  Time: -200.0:0.1:0.0

The concentration and importantly also the compartment values converge to a cycle of ii time units, i.e., values at time = 0.0 are (approximately) equal to values at time - ii = -5.0.

sim_many_bolus.sol(0.0) ≈ sim_many_bolus.sol(-5.0)

true

Thus, the compartment values

Depot_ss, Central_ss = sim_many_bolus.sol(0.0)

2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
  0.0
 11.498194694281786

at time = 0.0 are an approximation of the steady-state values that we are interested in. The simulation of the steady-state event can be performed by setting the compartement values at time = 0.0 to these steady-state values and continuing the simulation as usual.

model_approx = @model begin
    @pre begin
        Ka = 1.0
        CL = 1.0
        Vc = 10.0
    end

    @covariates BIOAV LAG

    @dosecontrol begin
        bioav = (; Central = BIOAV)
        lags = (; Central = LAG)
    end

    @init begin
        Depot = 0.0
        Central = 11.498194694
    end

    @dynamics Depots1Central1

    @observed begin
        conc = @. Central / Vc
    end
end

PumasModel
  Parameters: 
  Random effects: 
  Covariates: BIOAV, LAG
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: 
  Observed: conc

# Default simulation
obstimes = range(0, 50; length = 1_001)
sim = simobs(model, sub, (;); obstimes)

# Simulation based on approximation of the steady state
sim_approx = simobs(
    model_approx,
    # Same subject, but with a regular dose instead of the steady-state event
    Subject(
        sub;
        events = DosageRegimen(
            ss_dose.amt;
            ss_dose.time,
            ss_dose.ii,
            ss_dose.cmt,
            ss_dose.addl,
            ss_dose.evid,
            ss_dose.rate,
            ss_dose.duration,
            ss = 0,
        ),
    ),
    (;);
    obstimes,
)

sim_approx.observed.conc ≈ sim.observed.conc

true

Multiple Infusions

Steady-state events with multiple infusions are characterized by

1. amt > 0: positive dose amount

2. rate > 0: positive infusion rate

The time point (time) of such a steady-state event denotes the time point at which the infusion is started.

Also for this steady-state type we show an example, using the same one-compartment model as above:

model = @model begin
    @pre begin
        Ka = 1.0
        CL = 1.0
        Vc = 10.0
    end

    @covariates BIOAV LAG

    @dosecontrol begin
        bioav = (; Central = BIOAV)
        lags = (; Central = LAG)
    end

    @dynamics Depots1Central1

    @observed begin
        conc = @. Central / Vc
    end
end

PumasModel
  Parameters: 
  Random effects: 
  Covariates: BIOAV, LAG
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: 
  Observed: conc

Absorption rate and clearance are fixed to 1, volume of distribution is fixed to 10, and bioavailability and absorption lag time of the central compartment are set to covariates BIOAV and LAG, respectively. We inspect a subject

df = DataFrame(;
    id = 1,
    time = 0.0,
    amt = 5.0,
    cmt = :Central,
    rate = 10.0,
    ii = 5.0,
    addl = 3,
    evid = 1,
    ss = 1,
    conc = missing,
    BIOAV = 0.5,
    LAG = 4.0,
)
sub = only(read_pumas(df; observations = [:conc], covariates = [:BIOAV, :LAG]))

Subject
  ID: 1
  Events: 4
  Observations: conc: (n=0)
  Covariates: BIOAV, LAG

with dosage regimen

sub.events

and covariates

sub.covariates(0.0)

(BIOAV = 0.5,
 LAG = 4.0,)

at the time point of the dose record. We see that this dose record is a steady-state event (evid = 1 and ss = 1) at time = 0.0 with multiple infusion doses (amt = 5.0 > 0 and rate = 10.0 > 0) that are applied to the Central compartment every ii = 5.0 time units and repeated an additional addl = 3 times after the time point time = 0.0 of the dose record. Since ss = 1, the compartment amounts are reset to steady-state values at time = 0.0. Again the main question is: What are the steady-state values of the compartments at time = 0.0?

As for steady-state events with multiple bolus doses, the steady state can be viewed as the compartment values resulting from having applied the sequence of implied doses from the infinite past up until the time point time = 0.0 of the steady-state event. We approximate this infinite sequence of infusions with a sufficiently long sequence of infusions.

# Extract dose amount (`amt`), interdose interval (`ii`), compartment (`cmt`),
# rate (`rate`), and time point (`time`) of the steady-state dose
ss_dose = only(sub.events)
(; amt, ii, cmt, time) = ss_dose

# Number of implied doses that are applied
n_implied = 40

sim_many_infusion = simobs(
    model,
    Subject(
        sub;
        covariates = sub.covariates(time),
        events = DosageRegimen(
            amt;
            time = time - n_implied * ii,
            rate = ss_dose.rate,
            ii,
            addl = n_implied - 1,
            cmt,
        ),
    ),
    (;);
    obstimes = range(time - n_implied * ii, time; length = 50 * n_implied + 1),
)

SimulatedObservations
  Simulated variables: 
  Time: -200.0:0.1:0.0

The concentration and importantly also the compartment values converge to a cycle of ii = 5.0 time units, i.e., values at time = 0.0 are (approximately) equal to values at time - ii = -5.0.

sim_many_infusion.sol(0.0) ≈ sim_many_infusion.sol(-5.0)

true

Thus, the compartment values

Depot_ss, Central_ss = sim_many_infusion.sol(0.0)

2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
 0.0
 5.821563689981831

at time = 0.0 are an approximation of the steady-state values that we are interested in. The simulation of the steady-state event can be performed by setting the compartement values at time = 0.0 to these steady-state values and continuing the simulation as usual.

model_approx = @model begin
    @pre begin
        Ka = 1.0
        CL = 1.0
        Vc = 10.0
    end

    @covariates BIOAV LAG

    @dosecontrol begin
        bioav = (; Central = BIOAV)
        lags = (; Central = LAG)
    end

    @init begin
        Depot = 0.0
        Central = 5.821563689
    end

    @dynamics Depots1Central1

    @observed begin
        conc = @. Central / Vc
    end
end

PumasModel
  Parameters: 
  Random effects: 
  Covariates: BIOAV, LAG
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: 
  Observed: conc

# Default simulation
obstimes = range(0, 50; length = 1_001)
sim = simobs(model, sub, (;); obstimes)

# Simulation based on approximation of the steady state
sim_approx = simobs(
    model_approx,
    # Same subject, but with a regular dose instead of the steady-state event
    Subject(
        sub;
        events = DosageRegimen(
            ss_dose.amt;
            ss_dose.time,
            ss_dose.ii,
            ss_dose.cmt,
            ss_dose.addl,
            ss_dose.evid,
            ss_dose.rate,
            ss_dose.duration,
            ss = 0,
        ),
    ),
    (;);
    obstimes,
)

sim_approx.observed.conc ≈ sim.observed.conc

true

Comparison With NONMEM

Steady-state events in Pumas are designed to behave in the same way as in NONMEM. To illustrate this consistency, we compare the simulations produced by Pumas and NONMEM for different types of steady-state events.

For this comparison we use the examples and NONMEM simulation results in the nmtests repository, available under the GPL-2.0 license.

Simulations are performed with three different Pumas implementations of the one-compartment model used in NONMEM:

• A model with a pre-defined dynamical system with closed-form solutions of its dynamics and steady states ("Closed Form").

• A model with a linear ordinary differential equation (ODE) as dynamical system whose dynamics and steady states are computed with matrix exponentials ("Matrix Exponential").

• A model with a dynamical system whose dynamics and steady states are computed with numerical solvers ("Differential Equation Solver").

model_closedform = @model begin
    @pre begin
        Ka = 1.5
        CL = 1.1
        Vc = 20.0
    end

    @covariates LAGT DUR2 BIOAV

    @dosecontrol begin
        bioav = (Central = BIOAV,)
        lags = (Central = LAGT,)
        duration = (Central = DUR2,)
    end

    @dynamics Depots1Central1

    @observed begin
        CP = @. 1000 * Central / Vc
    end
end

PumasModel
  Parameters: 
  Random effects: 
  Covariates: LAGT, DUR2, BIOAV
  Dynamical system variables: Depot, Central
  Dynamical system type: Closed form
  Derived: 
  Observed: CP

model_matrixexp = @model begin
    @pre begin
        Ka = 1.5
        CL = 1.1
        Vc = 20.0
    end

    @covariates LAGT DUR2 BIOAV

    @dosecontrol begin
        bioav = (Central = BIOAV,)
        lags = (Central = LAGT,)
        duration = (Central = DUR2,)
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - CL / Vc * Central
    end

    @observed begin
        CP = @. 1000 * Central / Vc
    end
end

PumasModel
  Parameters: 
  Random effects: 
  Covariates: LAGT, DUR2, BIOAV
  Dynamical system variables: Depot, Central
  Dynamical system type: Matrix exponential
  Derived: 
  Observed: CP

model_diffeq = @model begin
    @options begin
        checklinear = false
    end

    @pre begin
        Ka = 1.5
        CL = 1.1
        Vc = 20.0
    end

    @covariates LAGT DUR2 BIOAV

    @dosecontrol begin
        bioav = (Central = BIOAV,)
        lags = (Central = LAGT,)
        duration = (Central = DUR2,)
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - CL / Vc * Central
    end

    @observed begin
        CP = @. 1000 * Central / Vc
    end
end

PumasModel
  Parameters: 
  Random effects: 
  Covariates: LAGT, DUR2, BIOAV
  Dynamical system variables: Depot, Central
  Dynamical system type: Nonlinear ODE
  Derived: 
  Observed: CP

For all three Pumas models, simulations coincide with the NONMEM output:

The Subject Constructor

The dosage regimen is only a subset of what we need to fully specify a subject. As mentioned above, we use the Subject type to represent individuals in Pumas. As seen below, then can either be constructed from primitives (time, events covariates, and observations) or read from a tabular format that will be mapped to these concepts. They can be constructed using the Subject constructor programmatically:

Subject

Subject(;id = "1",
         observations = NamedTuple(),
         events = nothing,
         time = observations isa AbstractDataFrame ? observations.time : nothing,
         covariates::Union{Nothing, NamedTuple} = nothing,
         covariates_time = observations isa AbstractDataFrame ? observations.time : nothing,
         covariates_direction = :left)

julia> Subject()
Subject
  ID: 1

julia> Subject(
         id = 20,
         events = DosageRegimen(200, ii = 24, addl = 2),
         covariates = (WT = 14.2, HT = 5.2),
       )
Subject
  ID: 20
  Events: 3
  Covariates: WT, HT

julia> Subject(covariates = (WT = [14.2, 14.7], HT = fill(5.2, 2)), covariates_time = [0, 3])
Subject
  ID: 1
  Covariates: WT, HT

or using read_pumas from tabular data:

read_pumas(filepath::AbstractString; missingstring = ["", ".", "NA"], kwargs...)
read_pumas(df::AbstractDataFrame; kwargs...)

You can also create Subjects from an existing Subject:

Subject(subject::Subject; id::AbstractString, events)

julia> s1 = Subject(; id = "AKJ491", events = DosageRegimen(1.0; time = 1.0))
Subject
  ID: AKJ491
  Events: 1

julia> Subject(s1; id = "AKJ492")
Subject
  ID: AKJ492
  Events: 1

The simulated output of a PumasModel as a result of simobs is Pumas.SimulatedObservations. Passing this simobs output into a vectorized Subject will result in a Population (an alias for a vector of Subjects) that is equivalent to the output of read_pumas. This is a convenient feature that allows one to simulate data and turn it back into a Population that can then be passed into a fit function:

Subject(simsubject::SimulatedObservations)

sims = simobs(model, pop, params)

Subject.(sims)

Understanding the Subject

The Subject is a central concept when it comes to simulation and fitting in Pumas. To do either of the two we need a model and data. The model is built using the PumasModel or PumasEMModel but carries no information about events, covariate values, or observed values. These all come from the Subject. In Pumas material you will also see Population mentioned and that is nothing but a collection of Subjects. Specifically, the collection type is a Vector such that if sub1, sub2, and sub3 are Subjects then [sub1, sub2, sub3] is a Population.

As mentioned, a model alone does not carry all the information needed to perform a simulation. In theory, we could divide the information into two categories:

1. Observations
2. Covariates

Both sets of information is always associated with a specific time point that refers to when the information was collected. In the case of dosing information the information is typically not collected but taken from protocol even if we would ideally always want to know the actual dosing time rather than the protocol time in a fitting context.

Covariates can further be divided into to categories: events and actual covariates. Events are typically dosing records including information such as additional doses, whether to treat the dose as if it has been given from an infinite past (steady-state dosing), and so on. Actual covariates would be information about the subject that does not directly affect events and that are not observations of derived variables in the statistical model. This could be demographics, formulation information that does not fit the dosing regimen, occasion indices, and so on. In Pumas, these two sets of covariate information fall into the two sub-categories events and covariates.

Observations

The data used for fitting model parameters and the data that will be simulated in simobs is specified in what we call observations

  Subject(id = "1", time = [...], observations = [...])

If we have data on record we will often use the read_pumas function that is designed to map tabular data to Subjects. This will be explained below. It is also possible to input data directly in the Subject constructor as follows:

Subject(
    id = "1",
    time = [0.0, 2.0, 4.0, 24.0],
    observations = (; dv = [3.0, 1.2, 0.8, 0.001]),
)

This requires the time vector to match up with the vectors inside the NamedTuple passed to observations. If we had two dependent variables we would simply add the second one to the NamedTuple

Subject(
    id = "1",
    time = [0.0, 2.0, 4.0, 24.0],
    observations = (; dv = [3.0, 1.2, 0.8, 0.001], dv2 = [0.1, 0.6, 0.9, 1.05]),
)

Where each element of each vector in the observations input matches the time vector. Some times there will not be a complete overlap between the two sets of sampling time points and in that case we need to add missing values. Let us say that dv2 is not sampled for the first two hours but it is additionally sampled at 48 hours then it might look like the following:

Subject(
    id = "1",
    time = [0.0, 2.0, 4.0, 24.0, 48.0],
    observations = (;
        dv = [3.0, 1.2, 0.8, 0.001, missing],
        dv2 = [missing, missing, 0.9, 1.05, 0.2],
    ),
)

Notice, that one dependent variable cannot have two observations at the same time point. Sometimes this can occur in data if pre-dose plasma samples are labelled as having the same time value as the dose and the post-dose sample. This is an example in the category mentioned above where the actual sample time is not reflected in the data set as you cannot realistically take a pre-dose sample, dose, and take a post-dose sample with a nanosecond. If this occurs, you need to shift the time point of the pre-dose sample slightly. Below, we add a pre-dose sample 3 minutes prior to dosing by setting time to -3/60 = -0.05:

Subject(
    id = "1",
    time = [-0.05, 0.0, 2.0, 4.0, 24.0, 48.0],
    observations = (;
        dv = [0.0, 3.0, 1.2, 0.8, 0.001, missing],
        dv2 = [missing, missing, missing, 0.9, 1.05, 0.2],
    ),
)

The same can be done in a DataFrame if using read_pumas.

As mentioned above, if we have data on record it is more often the case that we use read_pumas. The direct construction of Subjects is often related to simulation exercises where we want to simulate data instead. When simulating data, we can either choose to specify what variables to simulate and store in the simulated subject through the observations input or we can choose to not input anything and in that case all variables returned from @derived will be treated as relevant. Consider the following @derived block:

  @derived begin
    cp1 := @. Central / Vc
    dv1 ~ @. Normal(cp1, abs(cp1) * σ)
    cp2 := @. Metabolite / Vm
    dv2 ~ @. Normal(cp2, abs(cp2) * σ)
  end

This is a parent-metabolite model with proportional error models. cp1 and cp2 are suppressed from the output through := so the relevant variables to consider for data simulation is parent concentration dv1 and metabolite concentration dv2. If we want to simulate data from this model we need to prepare our subjects in one of the following two ways:

Subject(id = "1-both-1", time = [0.0, 1.0, 3.0, 4.0])

or

Subject(id = "1-both-2", time = [0.0, 1.0, 3.0, 4.0], observations = (:dv1, :dv2))

They will both contain the same information in this case because the omission of observations means "simulate everything" and in the second case we include observations and we mention everything. The list of variables to simulate should be given as a Tuple of Symbols (the names of the dv with a : in front).

However, if we wanted to simulate only parent or only metabolite we would have to do something like the following:

Subject(id = "1-parent", time = [0.0, 1.0, 3.0, 4.0], observations = (:dv1,))
Subject(id = "1-metabolite", time = [0.0, 1.0, 3.0, 4.0], observations = (:dv2,))

and then we could create a Subject from the output of simobs called with each subject by itself and we would get two different subjects back: one with only parent concentrations, and one with only metabolite concentrations.

PumasNDF

The PumasNDF is a specification for building a Population (an alias for a vector of Subjects) from tabular data. Generally this tabular data is given by a database like a CSV. The CSV has columns described as follows:

• id: the ID of the individual. Each individual should have a unique integer, or string.
• time: the time corresponding to the row. Should be unique per id, i.e. no duplicate time values for a given subject.
• evid: the event id. 1 specifies a normal dose event. 3 means it's a reset event, meaning that the value of the dynamical variable is reset to the amt at the dosing event. If 4, then the dynamical value and time are reset, and then a final dose is given. Defaults to 0 if amt is 0 or missing, and 1 otherwise.
• amt: the amount of the dose. If the evid column exists and is non-zero, this value should be non-zero. Defaults to 0.
• ii: the inter-dose interval. When addl is specified, this is the length of time to the next dose. For steady state events, this is the length of time between successive doses. Defaults to 0, and is required to be non-zero on rows where a steady-state event is specified.
• addl: the number of additional doses of the same amt to give. Defaults to 0.
• rate: the rate of administration. If 0, then the dose is instantaneous. Otherwise, the dose is administered at a constant rate for a duration equal to amt/rate. A rate=-2 allows the rate to be determined by Dose Control Parameters (DCP). Defaults to 0.
• ss: an indicator for whether the dose is a steady state dose. A steady state dose is defined as the result of having applied the dose with the interval ii infinitely many successive times. 0 indicates that the dose is not a steady state dose. 1 indicates that the dose is a steady state dose. 2 indicates that it is a steady state dose that is added to the previous amount. The default is 0.
• cmt: the compartment being dosed. Defaults to 1.
• duration: the duration of administration. If 0, then the dose is instantaneous. Otherwise, the dose is administered at a constant rate equal to amt/duration. Defaults to 0.
• Observation and covariate columns should be given as a time series of values of matching type. Constant covariates should be constant through the full column. Time points without a measurement should be denoted by a dot (.).

If a column does not exist, its values are imputed to be the defaults. Special notes:

• If rate and duration exist, then it is enforced that amt=rate*duration.
• All values and header names are interpreted as lower case.

using CSV
data = CSV.read(
    joinpath("pathtomyfile", "mydata.csv"),
    DataFrame;
    missingstring = ["", ".", "NA", "BQL"],
)

For more information check out the CSV documentation.

PumasNDF Checks

The read_pumas function performs some general checks on the provided data and informs the user about inconsistency in the data. It throws an error in case of invalid data, reporting the row number and column name causing the problem such that the user can resolve the issue.

Following is the list of checks applied by the read_pumas function.

Necessary columns in case of event and non-event data

In case of event_data=true (the default), the dataset must contain id, time, amt, and observations columns.

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    cmt = [1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

[ Info: The input has keys: [:id, :time, :cmt, :dv, :age, :sex, :evid]
ERROR: PumasDataError: The input must have: `id, time, amt, and observations` when `event_data` is `true`
[...]

In case of event_data=false, only the id column is required.

df = DataFrame(
    time = [0, 1],
    cmt = [1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex], event_data = false)

# output

ERROR: ArgumentError: column name "id" not found in the data frame; existing most similar names are: "dv" and "evid"
[...]

Event data without evid column

In case of event_data=true (the default), a warning is displayed if the provided dataset doesn't have an evid column.

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    cmt = [1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

┌ Warning: Your dataset has dose event but it hasn't an evid column. We are adding 1 for dosing rows and 0 for others in evid column. If this is not the case, please add your evid column.
│
└ @ Pumas
Population
  Subjects: 1
  Covariates: age, sex
  Observations: dv

Non-numeric/string entries in an observation column

If there are non-numeric or string entries in an observation column, read_pumas throws an error and reports row(s) and column(s) having this issue.

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    cmt = [1, 2],
    dv = [missing, "k@"],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: [1], row = [2], col = dv]  We expect the dv column to be of numeric type.
These are the unique non-numeric values present in the column dv: ("k@",)
[...]

Non-numeric/string entries in the amt column

Similarly, if there are non-numeric or string entries in an observation column, read_pumas throws an error and reports row(s) and column(s) having this issue.

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = ["k8", 0],
    cmt = [1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: [1], row = [1], col = amt]  We expect the amt column to be of numeric type.
These are the unique non-numeric values present in the column amt: ("k8",)
[...]

cmt must be a positive integer or valid string/symbol for non-zero evid data record

The cmt column should contain positive numbers or string/symbol identifiers of the compartment being dosed.

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    cmt = [-1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = cmt] column had invalid value: -1. Pumas does not currently support negative compartment values to turn off compartments.
[...]

amt can only be missing or zero when evid is zero

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 5],
    cmt = [1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 2, col = evid] amt can only be missing or 0 when evid is 0
[...]

amt can only be positive or zero when evid is 1

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [-10, 0],
    cmt = [1, 2],
    evid = [1, 0],
    dv = [10, 8],
    age = [45, 45],
    sex = ["M", "M"],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = evid] amt can only be positive or zero when evid is 1
[...]

Observations at the time of dose

Observations should be missing at the time of dose (or when amt > 0).

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    cmt = [1, 2],
    evid = [1, 0],
    dv = [10, 8],
    age = [45, 45],
    sex = ["M", "M"],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = dv] an observation is present at the time of dose in column dv. A blank record (`missing`) is required at time of dosing, i.e. when `amt` is positive.
[...]

Steady-state dosing requires ii column

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    ss = [1, 0],
    cmt = [1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: your dataset does not have ii which is a required column for steady state dosing.
[...]

Steady-state dosing with multiple infusions or bolus doses requires positive ii

In case of a steady-state event with multiple infusions or bolus doses (ss = 1 or ss = 2 and amt > 0), the value of the dose interval column ii must be non-zero.

If the rate column is not provided it is assumed to be zero.

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    ss = [1, 0],
    ii = [0, 0],
    cmt = [1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = ii] for steady-state dosing the value of the interval column ii must be non-zero but was 0
[...]

Steady-state dosing with constant infusion requires zero ii

In case of a steady-state event with constant infusion (ss = 1 or ss = 2, amt = 0 and rate > 0), the value of the dose interval column ii must be zero.

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [0, 0],
    ss = [1, 0],
    rate = [2, 0],
    ii = [1, 0],
    cmt = [1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = ii] for steady-state infusion the value of the interval column ii must be zero but was 1
[...]

Steady-state dosing with constant infusion requires zero addl

In case of a steady-state event with constant infusion (ss = 1 or ss = 2, amt = 0 and rate > 0), the value of the additional dose column addl must be zero.

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [0, 0],
    ss = [1, 0],
    rate = [2, 0],
    ii = [0, 0],
    addl = [5, 0],
    cmt = [1, 2],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
    evid = [1, 0],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = addl] for steady-state infusion the value of the additional dose column addl must be zero but was 5
[...]

Non-zero addl requires ii column

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    addl = [5, 0],
    cmt = [1, 2],
    evid = [1, 0],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: The dataset has an `addl` column specified with non-zero entries but does not have an `ii` column specified. It is not possible to specify the additional doses without providing the interval of time between them. Please specify which column the `ii` information is present in in your `read_pumas` call.
[...]

Non-zero addl requires positive ii

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    addl = [5, 0],
    ii = [0, 0],
    cmt = ["Depot", "Central"],
    evid = [1, 0],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = ii]  ii must be positive for addl > 0
[...]

Positive ii requires non-zero addl

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    addl = [0, 0],
    ii = [12, 0],
    cmt = ["Depot", "Central"],
    evid = [1, 0],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = addl]  addl must be positive for ii > 0
[...]

ii can only be missing or zero when evid is zero

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    addl = [5, 2],
    ii = [12, 4],
    cmt = ["Depot", "Central"],
    evid = [1, 0],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 2, col = evid]  ii can only be missing or zero when evid is 0
[...]

addl can only be positive or zero when evid is 1

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    addl = [-10, 0],
    ii = [12, 0],
    cmt = ["Depot", "Central"],
    evid = [1, 0],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = evid]  addl can only be positive or zero when evid is 1
[...]

evid must be nonzero when amt is positive or addl and ii are positive

When amt is positive, evid must be non-zero as evid=0 indicates an observation record.

df = DataFrame(
    id = [1, 1],
    time = [0, 1],
    amt = [10, 0],
    addl = [5, 0],
    ii = [12, 0],
    cmt = ["Depot", "Central"],
    evid = [0, 0],
    dv = [missing, 8],
    age = [45, 45],
    sex = ["M", "M"],
)
read_pumas(df; observations = [:dv], covariates = [:age, :sex])

# output

ERROR: PumasDataError: [Subject id: 1, row = 1, col = evid] amt can only be missing or 0 when evid is 0
[...]