There are a bunch of new features out now in the mrgsolve 1.0.0 release. Most of these are related to model syntax (what / how you write your model in the .mod or .cpp file) rather than model functionality.

library(mrgsolve)
library(tidyverse)

1 Use THETA(n) to refer to THETAn

All models will be able to use THETA(n) to refer to THETAn. For example:

[ nmxml ] 
run = 100
root = "cppfile"

[ main ] 
double CL = THETA(1) * pow(WT/70, 0.75);

You’ll notice that we referred to a nonmem run with [ nmxml ] so this syntax is meant be one small step toward making it easier to translate your model from the nonmem control stream to mrgsolve format. This is pretty simple, but will mean less time removing parentheses when copying over model code.

2 New plugin called `autodec`

Plugins are extra functionality that you can invoke when coding your model. For example, plugins exist for calculating time after dose, adding Rcpp functionality to your model etc and more.

One new plugin is called autodec. This plugin will automatically declare C++ variables for you. For example, I can code the snipped from the previous block using autodec plugin and mrgsolve will discover that CL needs to be declared:

[ plugin ] autodec

[ nmxml ] 
run = 100
root = "cppfile"

[ main ] 
CL = THETA(1) * pow(WT/70, 0.75);

The motivation behind this feature is to make it easier to code models with lots of variables when all your variables are simple double precision floating point numbers.

3 New plugin called `nm-vars`

This plugin provided a nonmem-like syntax for certain model elements, including

Fn, Dn, Rn, and ALAGn
A(n), A_0(n) and DADT(n)
T (used for SOLVERTIME in $ODE)

So your mrgsolve model translated from nonmem might look like this

see(modlib("nm-like"))

Building nm-like ... done.

$PROB Model written with some nonmem-like syntax features

$PLUGIN nm-vars autodec

$PARAM
THETA1 = 1, THETA2 = 21, THETA3 = 1.3, WT = 70, F1I = 0.5, D2I = 2
KIN = 100, KOUT = 0.1, IC50 = 10, IMAX = 0.9

$CMT @number 3

$PK
CL = THETA(1) * pow(WT/70, 0.75); 
V  = THETA(2); 
KA = THETA(3);

F1 = F1I;
D2 = D2I;
A_0(3) = KIN / KOUT;

$DES 
CP = A(2)/V;
INH = IMAX*CP/(IC50 + CP);
  
DADT(1) = -KA*A(1);
DADT(2) =  KA*A(1) - (CL/V)*A(2);
DADT(3) =  KIN * (1-INH) - KOUT * A(3);

$ERROR
CP = A(2)/V;

Notice that not all nonmem syntax is supported; just select data structures that use A and DADT to refer to compartments and differential equations. Also notice that we still require ; at the end of each line and we still need to use pow(base, exponent) to calculate exponent of some number.

4 Compartment

The $CMT block has two new options: @number and @prefix that let you quickly specify a numbered series of compartments in the model similar to nonmem. The default @prefix is A so that this code

$CMT @number 3

will put compartments A1, A2, A3 into your model.

5 New model in `modlib`: `nm-like`

You can see a model coded with these new syntax features in the internal model library; the model is called nm-like

mod <- modlib("nm-like")
see(mod)


Model file:  nm-like.cpp 
$PROB Model written with some nonmem-like syntax features

$PLUGIN nm-vars autodec

$PARAM
THETA1 = 1, THETA2 = 21, THETA3 = 1.3, WT = 70, F1I = 0.5, D2I = 2
KIN = 100, KOUT = 0.1, IC50 = 10, IMAX = 0.9

$CMT @number 3

$PK
CL = THETA(1) * pow(WT/70, 0.75); 
V  = THETA(2); 
KA = THETA(3);

F1 = F1I;
D2 = D2I;
A_0(3) = KIN / KOUT;

$DES 
CP = A(2)/V;
INH = IMAX*CP/(IC50 + CP);
  
DADT(1) = -KA*A(1);
DADT(2) =  KA*A(1) - (CL/V)*A(2);
DADT(3) =  KIN * (1-INH) - KOUT * A(3);

$ERROR
CP = A(2)/V;

Here, I’ve coded the THETAs in a parameter block; these would ordinarily come into the model via $NMXML or $NMEXT, which automatically import these parameters and estimates. But note that we can still refer to THETA(1) and others … this resolves to THETA1 regardless.

6 Audit

When mrgsolve loads your model, it checks the $ODE block to make sure you have differential equations for every compartment in your model. A long time ago, I started calling this check audit and it was something that was on by default but could be turned off via the call to mread().

Starting with version 1.0.0, audit can be controlled at the $ODE block level like this

[ ode ] @audit

dxdt_A1 = ...
dxdt_A2 = ...

to explicitly tell mrgsolve to audit the equations or

[ ode ] @!audit

dxdt_A1 = ...
dxdt_A2 = ...

to disable the audit. The audit is on by default and should only be turned off when you have written some special construct in [ ode ] to code the differential equations (e.g. odes are written programmatically).

7 New time spacer for event sequences

The new feature is an ii spacer for sequences of event objects.

Using event objects, we can create complicated dosing regimens from simpler ones. For example, we might have a week of 50 mg dosing followed by a week of 100 mg dosing and then up to 200 mg dosing. We can accomplish that with an event sequence

First create the individual pieces

a <- ev(amt = 50,  ii = 24, addl = 6)
b <- ev(amt = 100, ii = 24, addl = 6)
c <- ev(amt = 200, ii = 24, addl = 6)

Then put them in a sequence

data <- ev_seq(a, b, c)
data

Events:
  time amt ii addl cmt evid
1    0  50 24    6   1    1
2  168 100 24    6   1    1
3  336 200 24    6   1    1

When they are sequenced, you’ll see that the second piece (b) starts one dosing interval after the last dose in the first piece (a).

We can put a 24 hour spacer between a and b

seq(a, wait = 24, b)

Events:
  time amt ii addl cmt evid
1    0  50 24    6   1    1
2  192 100 24    6   1    1

Here, the last dose in a is given, we wait one dosing interval (24 hours) then wait another 24 hours (via wait) and then start b.

Rather than using wait, we can use ii to specify the amount of time from the last dose in a to the first dose in b. So if we want to wait 3 days between the last dose in a and the first dose in b

ev_seq(a, ii = 3*24, b)

Events:
  time amt ii addl cmt evid
1    0  50 24    6   1    1
2  216 100 24    6   1    1

Notice that this same behavior can be achieved with wait but using ii might be easier to interpret in some cases.

8 Collapse matrices

Sometimes in your model, you have multiple OMEGA or SIGMA blocks like this

code <- '
$OMEGA  @name first
1 2 3

$OMEGA @name second
4
'
mod <- mcode("collapse", code, compile = FALSE)

omat(mod)

$first
    [,1] [,2] [,3]
1:     1    0    0
2:     0    2    0
3:     0    0    3

$second
    [,1]
4:     4

If I want to update these values but I only have a 4x4 matrix

mat

     [,1] [,2] [,3] [,4]
[1,]  0.1  0.0  0.0  0.0
[2,]  0.0  0.2  0.0  0.0
[3,]  0.0  0.0  0.3  0.0
[4,]  0.0  0.0  0.0  0.4

I can’t do it; mrgsolve wants a 3x3 matrix for the first slot and 1x1 matrix for the second.

A new function will collapse the model matrix into larger matrices

mod <- collapse_omega(mod)
omat(mod)

$...
    [,1] [,2] [,3] [,4]
1:     1    0    0    0
2:     0    2    0    0
3:     0    0    3    0
4:     0    0    0    4

mod <- update(mod, omat = mat)

Warning: The following argument was passed to `mrgsolve::update()`, but
it is either an invalid name (check your spelling?) or not an
eligible attribute for update:
✖ omat

omat(mod)

$...
    [,1] [,2] [,3] [,4]
1:     1    0    0    0
2:     0    2    0    0
3:     0    0    3    0
4:     0    0    0    4

You can’t split matrices up into smaller chunks, but if there are more than two blocks, you can select which blocks to join

code <- '
$OMEGA  @name first
1 2 3

$OMEGA @name second
4

$OMEGA @name third
5 6 7 8
'
mod <- mcode("collapse2", code, compile = FALSE)

omat(mod)

$first
    [,1] [,2] [,3]
1:     1    0    0
2:     0    2    0
3:     0    0    3

$second
    [,1]
4:     4

$third
    [,1] [,2] [,3] [,4]
5:     5    0    0    0
6:     0    6    0    0
7:     0    0    7    0
8:     0    0    0    8

mod <- collapse_omega(mod, range = c(2, NA), name = "remainder")
omat(mod)

$first
    [,1] [,2] [,3]
1:     1    0    0
2:     0    2    0
3:     0    0    3

$remainder
    [,1] [,2] [,3] [,4] [,5]
4:     4    0    0    0    0
5:     0    5    0    0    0
6:     0    0    6    0    0
7:     0    0    0    7    0
8:     0    0    0    0    8

9 Deprecated simeta(n) and simeps(n)

We recently rolled out a feature where the user could call simeta() and pass an integer that would indicate a single ETA() to update, leaving the others alone. This proved to be more than a little dangerous when $OMEGA had off-diagonal elements. We will be deprecating this feature starting with version 1.0.0. When you load a model that contains simeta(n) or simeps(n), there will be a warning about the syntax. The warning includes some code you can write to silence the warning, but this will be temporary. In a future release, we will always warn when this syntax us used and eventually generate an error. It is recommended to always use simeta() or simeps() when re-simulating random effects.