Exercise 1.2: PK Model Development

Task

Develop the non compartment model and compartment model in R and MCSim

Non compartment model

\[C(t) = \frac{F\cdot D \cdot ka}{V \cdot BW \cdot (k_a - k_e)}\cdot(\mathrm{e}^{-k_et}-\mathrm{e}^{-k_a{t}})\]

Compartment model

\[\frac{dA_{gut}}{dt}=-k_aA_{gut}\]

\[\frac{dA_{cpt}}{dt}=k_aA_{gut}-k_eA_{cpt}\] \[C_{cpt}=\frac{A_{cpt}}{V \cdot BW}\]

Parameters

\(t\): Time (hr)

\(C\): Concentration (mg/L)

\(F\): Absorption fraction (-)

\(D\): Ingestion dose (mg)

\(k_a\): Absorption rate (/hr)

\(k_e\): Elimination rate (/hr)

\(V\): Distribution volume (L/kg BW)

\(BW\): Body weight (kg)

\(A_{gut}\): Amount in gut lumen (mg)

Introduction of model development in MCSim

# Model description file (this is a comment)

<Global variable specifications>

States = { 
  <state variables for the model, such as quantity> 
}
Outputs = { 
  <output variables, such as concentration> 
}
Inputs = { 
  <input variables, such as exposure dose>
}
Initialize {
  <Equations for initializing or scaling model parameters>
}
Dynamics {
  <Equations for computing derivatives of the state variables>
}
CalcOutputs {
  <Equations for computing output variables>
}
End. # mandatory ending keyword

General syntax

Variable assignments

<variable-name> = <constant-value-or-expression> ;

Colon conditional assignments

<variable-name> = (<test> ? <value-if-true> : <value-if-false>);

For example

Adjusted_param = (Input_var > 0.0 ? Param * 1.1 : Param);

Syntax of (simulation) input-file

For the basic simulation

# Input-file (text after # are comments)
<Global assignments and specifications>
Simulation {
  <Specifications for first simulation>
}
Simulation {
  <Specifications for second simulation>
}
# Unlimited number of simulation specifications
End. # Mandatory End keyword. Everything after this line is ignored

1.1 Define the state variables States{}

States are variables for which a first-order differential equation is defined in the Dynamics{} section. Remember, when naming the variables, be sure to let all variable names begin with a capital letter followed by meaningful lower case subscripts. In States{} section, each variable needs to separate by comma ,. In addition, remember to provide the unit that is used for the state variables.

1.2 Define the outputs variables Outputs{}

Outputs are dependent model variables (obtainable at any time as analytical functions of the states, inputs or parameters) that do not have dynamics. They must receive assignments in either the Dynamics or CalcOutputs sections. When constructing the model, we need to have an output variable for checking mass balance. Then, assign the variable that we are interested in or have an experiment to make the comparison, such as the blood concentration.

1.3 Define the input variable(s) Inputs{}

Inputs are variables independent of the other variables, but eventually varying with time (for example an exposure concentration to a chemical). GNU MCSim has various functions to define the different exposure types as:

• PerDose(): # specifies a periodic input of constant

  PerDose(<magnitude>, <period>, <initial-time>, <exposure-time>);

• PerExp(): # specifies a periodic exponential input.

  PerExp(<magnitude>, <period>, <initial-time>, <decay-constant>);

• PerTransit(): #models a delayed input mechanism (e.g., gut absorption)

  PerTransit(<magnitude>, <period>, <initial-time-in-period>, <decay-constant>, <number-of-input-compartments>);

• NDoses(): # specifies a number of stepwise inputs of variable magnitude and their starting times

  NDoses(<n>, <list-of-magnitudes>, <list-of-initial-times>);

• Spikes(): # specifies a number of instantaneous inputs of variable magnitude and their exact times of occurrence.

  Spikes(<n>, <list-of-magnitudes>, <list-of-times>);

1.4 Parameters of the model

Here we need to define the “default” value of input parameters. Unlike previous variables that need to use curly brackets to group these variables, the input parameters can be put in any places without the restriction. In addition, don’t forget to put a semicolon (;) after the definition of each parameter. If the input parameter doesn’t have a default value, just put a semicolon after the name of the parameter or used an arbitrary number.

The PBPK parameters always include the parameters that can describe absorption, distribution (e.g., partition coefficient), metabolism (e.g., Michaelis-Menten constant), and elimination (e.g., rate Constant). These parameters can be classified to chemical-specific (e.g., molecular weight) and anatomy (e.g., body weight, blood flow, tissue weight/volume) parameter. In addition, the unit conversion factor can define in this part.

The parameters of the PBPK model might include:

• Unit conversion factor

• Exposure modeling parameters

• Chemical-specific parameters

• Physiological parameters

• Pharmacokinetic parameters

• Scale parameters (no default value)

1.5 Define the parameter initialization and scaling Initialize{}

This section is used to scale specific model parameters and resolve dependencies between parameters. Generally, the scaling involves a change of units or conversion from percentage to actual units.

1.6 Define the dynamics of the simulation Dynamics{}

The equations given in this section will be called by the integration routines at each integration step. It includes specification of differential equations.

The derivative of a state variable is defined using the dt() operator, as shown here:

dt(state-variable) '=' constant-value-or-expression ';'

1.7 Define the output calculations CalcOutputs{}

In this section, the equations are given in this section will be called by the simulation program at each output time specified by a Print() or PrintStep() statement. In this way, output computations are done efficiently, only when values are to be saved. Here is the example to assign the variables in CalcOutputs{} section. Again, don’t forget the semicolon and put the keyword End. at the end of model.

Same as model file, we highly recommend to use comments (#) to annotate your code in the beginning. This comment might include:

• Brief description of this input file
• What model file that is used to conduct the simulation
• Others…

1.8 Setting integrator Integrate()

GNU MCSim provides three types of integrator that include Lsodes, Cvodes, and Euler. Here, we used Lsodes (Livermore Solver for Ordinary Differential Equation) as a primary solver, which was originated from the SLAC Fortran library. The Lsodes function uses adaptive numerical methods to advance a solution to a system of ordinary differential equations one time-step, given values for the variables Y and X. It solves the initial value problem for stiff or non-stiff systems of ordinary differential equations (ODE).

The syntax for Lsodes is: Integrate(Lsodes, <rtol>, <atol>, <method>);

where <rtol> is a scalar specifying the relative error tolerance for each integration step. The scalar <atol> specifies the absolute error tolerance parameter. Those tolerances are used for all state variables. The estimated local error for a state variable y is controlled so as to be roughly less (in magnitude) than rtol*|y| + atol. Thus the local error test passes if, for each state variable, either the absolute error is less than , or the relative error is less than <rtol>. Set <rtol> to zero for pure absolute error control, and set <atol> to zero for pure relative error control. Caution: actual (global) errors may exceed these local tolerances, so choose them conservatively. Decreasing the tolerance leads to more accurate results but at the cost of significant increase in time taken.

The <method> flag should be 0 (zero) for non-stiff differential systems and 1 or 2 for stiff systems. You should try flag 0 or 1 and select the fastest for equal accuracy of output unless insight from your system leads you to choose one of them a priori. If you specify 2 you should provide the Jacobian of your differential system. The good starting point for and is about 1e-6. The default integration method in Lsodes with parameters 1e-5, 1e-7 and 1.

1.9 Setting paramerer(s)

For the simple simulation, we can assign the given parameter values in the input file to replace the default values in the model file. The parameter can be set in the simulation{} section as a local setting or outside the section as a global setting. Same as the parameter setting in the model file, use = to define the given value of parameters and put ; after the definition.

Input functions

These functions can use to different exposure types

- PerDose(): # specifies a periodic input of constant

    PerDose(<magnitude>, <period>, <initial-time>, <exposure-time>);

  
- PerExp(): # specifies a periodic exponential input.

    PerExp(<magnitude>, <period>, <initial-time>, <decay-constant>);  

  
- PerTransit(): models a delayed input mechanism  

    PerTransit(<magnitude>, <period>, <initial-time-in-period>, 
              <decay-constant>, <number-of-input-compartments>);  
    
              
- NDoses(): specifies a number of stepwise inputs of variable magnitude and their starting times
    
    NDoses(<n>, <list-of-magnitudes>, <list-of-initial-times>);


- Spikes(): specifies a number of instantaneous inputs of variable magnitude and their exact times of occurrence.
    
    Spikes(<n>, <list-of-magnitudes>, <list-of-times>);

1.10 Setting simulation simulation{}

After the global specifications, we can define the local parameters. In addition, we need to define the output through Print() or PrintStep() functions. The arguments of Print() and PrintStep() are the comma-separated list of variable names (at least one and up to MAX_PRINT_VARS, which is set to 10 by default) and followed by given time-points (limited by the available memory at run time) as:

Print(<identifier1>, <identifier2>, ..., <time1>, <time2>, ...);

PrintStep(<identifier1>, <identifier2>, ..., <start-time>, <end-time>, <time-step>);

Note: The output variables cannot receive assignments in simulation input files.

Solution

Non compartment model

• Non compartment model (R)

ffpk <- function(IngDose, Fgutabs, kgutabs, kelim, Vdist, BW, t){
    A <- (IngDose * Fgutabs * kgutabs)/(Vdist * BW * (kgutabs - kelim))
    Conc <- A * exp(-kelim * t) - A * exp(-kgutabs * t)
    return(Conc)
}

• Non compartment model (MCSim)

Outputs = {
  A,
  Conc,
}

# ingested input             (mg)
IngDose = 1.0; 

# Absorption fraction        (-)
Fgutabs = 1.0; #

# Absorption rate constant   (/h)
kgutabs = 1.0; #

# Elimination rate constants (/h)
kelim = 1.0;

# Distribution volumes       (L/kg BW)
Vdist = 1.0;

# Body weight                (kg)
BW = 70;

CalcOutputs { 
  A = (IngDose * Fgutabs * kgutabs)/(Vdist * BW * (kgutabs - kelim));
  Conc = A * exp(-kelim * t) - A * exp(-kgutabs * t);
}

End.

Compartment model

• Compartment model (R)

pbtk1cpt <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    dAgut = - kgutabs * Agut
    dAcpt = kgutabs * Agut - ke * Acpt

    Ccpt = Acpt / Vdist * BW
    
    list(c(dAgut, dAcpt), "Ccpt" = Ccpt) 
  })
}

• Compartment model (MCSim)

States  = { 
  Acpt      # Quantity in central compartment (mg)
};  

Inputs = { 
  Oral_input  # Drug input (mg)
}; 

Outputs = {
  Ccpt,     # Concentration in central compartment (mg/l)
};

# Oral input modeling
Absdose;
IngDose    = 1.0; # ingested input (mg)
Fgutabs    = 1.0; #
Period     = 0.0; # period of the exposure/no exposure cycle (h)
Tlag       = 0.0; # Absorption lagtime (h)
kgutabs    = 0.1; # Intestinal absortion rate (/h)
Oral_input = PerExp (Absdose, Period, Tlag, kgutabs);

# Distribution volumes (L/kg BW)
Vdist = 1;

# Body weight (kg)
BW = 70;

# Elimination rate constants (/h)
kelim = 1;

Initialize{ Absdose = IngDose * Fgutabs; }

Dynamics { 
  dt (Acpt) = kgutabs * Oral_input - kelim * Acpt;
}

CalcOutputs { 
  Ccpt = Acpt  / Vdist * BW;
}

End.