1 Model components

This chapter details the different components of a model in mrgsolve. Each component listed here is maintained within a “model object”. This is an updatable S4 object in R that contains all of the basic information required to properly configure and simulate from the model.

1.1 Parameter list

The parameter list is an updatable set of name-value pairs. Referencing the name of an item in the parameter list will substitute the current value associated with that name. While the name “parameter” may have a certain connotation in the modeling world, in mrgsolve a “parameter” could be any category of numeric data: covariates (e.g. WT, AGE, SEX), flags, other numeric data that we commonly call “parameter” (e.g. CL or VC).

The parameter list is declared in the code block $PARAM. While there may be multiple $PARAM blocks in a model, these are condensed to a single parameter list stored in the model object. The names and numbers of all parameters in the model must be declared at the time that the model is compiled. Also, a default value for each parameter must be declared at model compile time, but the value of each parameter may be updated in one of several ways.

The parameters in a model object can be queried or updated with the param() function.

See also: 2.2.4, ?param in the R help system after loading mrgsolve.

1.1.1 Central role of parameters in planning simulations

The data items in the parameter list are more than just values associated with a name. When an name is added to the parameter list, that name becomes a key word that mrgsolve will start to recognize in input data sets or when manipulating the model object.

For example, when you want to include a covariate in the model, say weight (WT), you’ll include a column in the data set called WT that will indicate the weight of this or that patient. It is crucial that you also list WT in $PARAM with some default value. It helps if that value is sensible too. When mrgsolve receives the data set prior to simulating, the WT column is matched up with the WT parameter name. As mrgsolve works its way through the input data set (from person to person or from time to time), the value of WT is updated so that the symbol WT in $MAIN or $ODE or $TABLE always points to the value of WT. If the WT name is not in the parameter list, it won’t matter if it is in the data set or not. Only listing a name in $PARAM gets it “into the game”.

Understanding the parameter update mechanism is very important for planning complicated simulations with mrgsolve. Please see the information in 3.1 and in 11.3.

1.2 Compartment list

Like the parameter list, the compartment list is a series of name-value pairs. The compartment list defines the number, names, and initial values of each compartment in the model. The names, numbers, and order of the compartment in a model is established at the time of model compile and changes to the compartment list require re-compilation of the model.

Compartments are declared in one of two code blocks: $INIT and $CMT. Nominal initial values must be supplied for each compartment. The main difference between $INIT and $CMT is that $CMT assumes a default initial value of 0 for each compartment; thus only compartment names are entered. When using $INIT, both names and values must be explicitly stated for each compartment.

The initial values for each compartment can be queried with the init() function. There are several different ways to set the initial conditions in a model; section 11.2 illustrates several of these.

See also: section 11.2 and ?init in the R help system after loading mrgsolve.

1.3 Simulation time grid

The mrgsolve model object stores the parameters for the series of time points to be output for a a simulation. This is the default output time grid that will be used if not over-ridden by another mechanism.

The elements of the simulation time grid are: start, end, delta and add. start, end, delta are passed to seq() as from, to, and by, respectively. add is any arbitrary vector of additional times to simulate.

The simulation time grid in a model object may be queried with the stime() function or by printing the model object to the R console.

See also section 3.2 for discussion of the simulation time grid and input data sets and 1.3.1 and 11.4 for using time grid objects .

1.3.1 tgrid objects

A tgrid object has start, end, delta and add attributes. This object is independent of the model object. tgrid objects may be created and combined to create complex sampling designs.

See section 11.4 for examples and usage.

1.4 Solver settings

mrgsolve uses the DLSODA solver from ODEPACK. Several of the settings for that solver are stored in the model object and passed to the solver when the problem is started. Settings include: atol, rtol, maxsteps, hmax, hmin, ixpr, mxhnil.

1.4.1 atol

Absolute tolerance parameter. Adjust this value lower when you see state variables (compartments) that are becoming very small and possibly turning negative. For example:

mod <- modlib("viral1", end = 144) 
. Loading model from cache.
out <- mrgsim_e(mod, ev(amt = 1000)) %>% filter(V < 0)

out
. # A tibble: 6 × 8
.      ID  time expos        T         I         V  logV logChange
.   <dbl> <dbl> <dbl>    <dbl>     <dbl>     <dbl> <dbl>     <dbl>
. 1     1  95    1000 5187630. -2.74e-10 -1.79e-12   NaN       NaN
. 2     1  97.5  1000 5298189. -3.18e-10 -2.09e-12   NaN       NaN
. 3     1 120    1000 6252828. -1.55e-10 -1.01e-12   NaN       NaN
. 4     1 122.   1000 6354547. -2.90e-10 -1.90e-12   NaN       NaN
. 5     1 125    1000 6455422. -1.97e-10 -1.30e-12   NaN       NaN
. 6     1 128.   1000 6555459. -7.81e-11 -5.21e-13   NaN       NaN

Adjusting atol to 1E-20 or 1E-30 will prevent this.

mrgsim_e(mod, ev(amt = 1000), atol = 1E-20)  %>% filter(time %in% out$time)
. # A tibble: 6 × 8
.      ID  time expos        T        I        V  logV logChange
.   <dbl> <dbl> <dbl>    <dbl>    <dbl>    <dbl> <dbl>     <dbl>
. 1     1  95    1000 5187630. 1.33e-11 8.72e-14 -13.1     -18.7
. 2     1  97.5  1000 5298189. 5.13e-12 3.36e-14 -13.5     -19.2
. 3     1 120    1000 6252828. 1.01e-15 6.62e-18 -17.2     -22.9
. 4     1 122.   1000 6354547. 3.94e-16 2.58e-18 -17.6     -23.3
. 5     1 125    1000 6455422. 1.54e-16 1.01e-18 -18.0     -23.7
. 6     1 128.   1000 6555459. 5.99e-17 3.92e-19 -18.4     -24.1

1.4.2 rtol

Relative tolerance parameter. Adjust this value lower when you want more precision around the calculation of state variables as the system advances.

1.4.3 maxsteps

This is the maximum number of steps the solver will take when advancing from one time to the next. If the solver can’t make it in maxsteps it will stop and give an error message like this:

DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
      taken on this call before reaching TOUT     
In above message, I = 
[1] 2000
In above message, R = 
[1] 0.0004049985
DLSODA-  ISTATE (=I1) illegal.
In above message, I = 
[1] -1
DLSODA-  Run aborted.. apparent infinite loop.    
Error in (function (x, data, idata = null_idata, carry.out = character(0),  : 
  error from XERRWD

You might see this when you have to integrate along time between records in a data set. There isn’t necessarily a problem, but the solver might have to advance over many doses to get to the next record and it only has a limited number of steps it can take between those records before it stops with this error.

When you see this, increase maxsteps to 50000 or larger.

But keep in mind that sometimes the solver can’t make it to the next record because there are issues with the model. It might take thousands of steps to make it 24 hours down the road. In that case, go back to the model code and look for problems in how it is coded.

1.4.4 hmax

The maximum step size. By default, the solver will take steps of different sizes based on what is happening in the simulation. Setting hmax tells the solver not to take a step larger than that value. So in a model where time is in hours, reducing hmax to 0.1 will prevent the solver from taking a step larger than 0.1 hours as it tries to advance to the next time. The will slow down the simulation a bit. But sometimes helpful when the solver starts taking large steps. We don’t recommend using this routinely; for most applications, it should be reserved for troubleshooting situations. If your model doesn’t give the results that you want without setting hmax, we’d recommend a new setup where this isn’t needed.

1.4.5 hmin

The minimum step size. Only set this if you know what you’re doing.

1.4.6 ixpr

A flag to enable printing messages to the R console when the solver switches between non-stiff and stiff solving modes. Rarely used.

1.4.7 mxhnil

The maximum number of messages printed when the model is solving. If you have a lot of messages, keep working on your model code.

1.5 Functions

There are four C++ functions that mrgsolve creates and manages: PREAMBLE, MAIN, ODE, TABLE. Each function is created from an entire code block in the model specification file. The user is responsible for writing correct C++ code in each of these blocks. mrgsolve will parse these blocks and augment this code with the necessary elements to create the C++ function.

These functions may be specified in any order in the model specification file, but there is a specific calling order for these functions. Recognizing and understanding this calling order will help understand how the different pieces of the model specification fit together.

Just prior to starting the problem, mrgsolve calls $PREAMBLE. Then, during advance from time T1 to T2, first $MAIN is called, then $ODE is called repeatedly as the solver finds the values of state variables at T2, and, once the solution is found, $TABLE is called to calculate derived quantities at T2 and to specify variables that should be included in the model output. So, it is helpful to write model specification files in the order:

  1. $PREAMBLE - called only once just prior to processing the first record of the data set
  2. $MAIN - before advancing the system
  3. $ODE - the system advances to T2
  4. $TABLE - after advancing the system

But the order in which they are coded will not affect model compilation or the simulation result.

1.5.1 The $PREAMBLE function

The PREAMBLE function gets called only once, just prior to processing the first record of the data set. This function is composed of C++ code and is used to initialize variables and get them set up prior to starting on the problem.

See 2.2.13 for details.

1.5.2 The $MAIN function

The MAIN function gets called at least once before the the solver advances from the current time (T1) to the next time (T2). In the MAIN function, the user may:

  • Set initial conditions for any compartment
  • Derive new variables to be used in the model
  • Write covariate models
  • Add between-subject variability to quantities to structural model parameters (e.g. CL or VC).

In addition to getting called once per record, the MAIN function may be called several times prior to starting the simulation run. The MAIN function is also called whenever the user queries the compartment list.

See 2.2.7 for details.

1.5.3 The $ODE function

The ODE function is where the user writes the model differential equations.
Any derived quantity that depends on a state variable and is used to advanced the system must be calculated inside $ODE. But, this function is called repeatedly during the simulation run, so any calculation that can be moved out of $ODE (for example: to $MAIN) should be.

See 2.2.9 for details.

1.5.4 The $TABLE function

The TABLE function is called after the solver advances in time. The purpose of TABLE is to allow the user to interact with the values of the state variables after advancing, potentially derive new variables, and to insert different outputs into the table of simulated results.

See 2.2.11 for details.

1.6 Random effect variances

The mrgsolve model object keeps track of a arbitrary number of block matrices that are used to simulate variates from multivariate normal distributions. Users can specify OMEGA matrices for simulating between-subject random effects (one draw per individual) or SIGMA matrices for simulating within-subject random effects (one draw per observation).

The user may use the revar() function to query both OMEGA and SIGMA.

1.6.1 OMEGA

The matrices are specified in $OMEGA blocks in the model specification file.

OMEGA may be queried or updated with the omat() function.

1.6.2 SIGMA

The matrices are specified in $SIGMA blocks in the model specification file.

SIGMA may be queried or updated by the smat() function.