pmxNODE is a package to facilitate the application of neural ordinary differential equations and the integration of neural networks in pharmacometric modeling software, including Monolix, NONMEM, and nlmixr2. It allows to utilize NN functions that are not commonly available in said pharmacometric software. The way pmxNODE works is that it translates NN functions into explicit equations that describe the calculations within a neural network, described in the publication Low-dimensional Neural ODEs accounting for inter-individual variability implemented in Monolix and NONMEM by Bräm, Steiert, Steffens, Pfister, and Koch (https://doi.org/10.1002/psp4.13265). It further allows to automatically initialize the neural network parameters, following the initialization approach presented in the publication above.
The general workflow for Monolix with pmxNODE consists of few steps:
Write a Monolix model file, as you would do normaly. However, you can utilize NN functions in your model code, e.g., for complex or unknown model parts.
Convert the written model file with the nn_converter_mlx function from the pmxNODE
package. This function can also directly generate a .mlxtran
file with automatically initialized model parameters.
We suggest to fit the model to the data first without
inter-individual variability on neural networks parameters (argument
pop_only = True in the nn_converter_mlx function) and add the
random effects in a second run, where parameters were initialized with
last estimates.
We first need to load the pmxNODE package. Since the automated generation of .mlxtran files is done with the lixoftConnectors package, this also needs to be loaded and initialized. The installation instructions are provided on the MonolixSuite-website.
Some examples are available in the pmxNODE package. To see all
example files, you can use the get_example_list function. To copy an
example to a folder of your choice, you can use the copy_example function. After calling the
copy_examples function, two files
should be in the target folder, a data file and a Monolix model
file.
get_example_list()
copy_examples(
target_folder = "~/pmxNODE",
example_nr = 1,
example_software = "Monolix"
)Let’s have a look at the model file:
DESCRIPTION:
A NODE for iv drug administration with concentration and time-dependent NN.
Volume of distribution is estimated; initial value for V is 2.
[LONGITUDINAL]
input = {V=2}
PK:
; PK model definition
depot(target=C)
EQUATION:
ddt_C = NNc(state=C,min_init=0.5,max_init=10) + amtDose * NNct(state=t,min_init=0.5,max_init=10,time_nn=TRUE)
Cc = C/V
OUTPUT:
output = Cc Before fitting the model, it needs to be converted with the nn_converter_mlx function. In addition to
the path/file name of the unconverted Monolix model and the argument on
including inter-individual variability for the neural network
parameters, a Monolix .mlxtran file can be automatically
generated with the gen_mlx_file = TRUE argument. In order to do
so, a data file and the header types must be provided. If no file name
for the new Monolix file is given through the mlx_name
argument, the Monolix file name is automatically generated based on the
Monolix model file. Note that a suffix is added to the file name, either
_pop if pop_only = TRUE or _ind if
pop_only = FALSE.
nn_converter_mlx(mlx_path = "~/pmxNODE/mlx_example1_model.txt",
pop_only = TRUE,
gen_mlx_file = TRUE,
mlx_name = "~/pmxNODE/mlx_example1",
data_file = "~/pmxNODE/data_example1_mlx.csv",
header_types = c("id","time","amount","observation"))Now the converted model file has included all the code needed for the NODE:
DESCRIPTION:
A NODE for iv drug administration with concentration and time-dependent NN.
Volume of distribution is estimated; initial value for V is 2.
[LONGITUDINAL]
input = {V,Wc_11,Wc_12,Wc_13,Wc_14,Wc_15,bc_11,bc_12,bc_13,bc_14,bc_15,Wc_21,Wc_22,Wc_23,Wc_24,Wc_25,bc_21,Wct_11,Wct_12,Wct_13,Wct_14,Wct_15,bct_11,bct_12,bct_13,bct_14,bct_15,Wct_21,Wct_22,Wct_23,Wct_24,Wct_25,bct_21}
PK:
; PK model definition
depot(target=C)
EQUATION:
hc_1 = max(0,Wc_11 * C + bc_11)
hc_2 = max(0,Wc_12 * C + bc_12)
hc_3 = max(0,Wc_13 * C + bc_13)
hc_4 = max(0,Wc_14 * C + bc_14)
hc_5 = max(0,Wc_15 * C + bc_15)
NNc = Wc_21 * hc_1 + Wc_22 * hc_2 + Wc_23 * hc_3 + Wc_24 * hc_4 + Wc_25 * hc_5 + bc_21
hct_1 = max(0,-(Wct_11^2) * t + bct_11)
hct_2 = max(0,-(Wct_12^2) * t + bct_12)
hct_3 = max(0,-(Wct_13^2) * t + bct_13)
hct_4 = max(0,-(Wct_14^2) * t + bct_14)
hct_5 = max(0,-(Wct_15^2) * t + bct_15)
NNct = Wct_21 * hct_1 + Wct_22 * hct_2 + Wct_23 * hct_3 + Wct_24 * hct_4 + Wct_25 * hct_5 + bct_21
ddt_C = NNc + amtDose * NNct
Cc = C/V
OUTPUT:
output = Cc The model can be automatically run from R with the function run_mlx from the pmxNODE package.
In order to get the parameter estimations from the population fit
(without inter-individual variability), the pre_fixef_extractor_mlx function can be
utilized.
These parameter estimates can be given as additional argument
pre_fixef to the nn_converter_mlx function. To include
inter-individual variability, the population argument is set to false
(pop_only = FALSE) in the nn_converter_mlx function.
The final model with inter-individual variability can then be fitted
again with the run_mlx function.
est_parms <- pre_fixef_extractor_mlx("~/pmxNODE/mlx_example1_pop.mlxtran")
nn_converter_mlx(mlx_path = "~/pmxNODE/mlx_example1_model.txt",
pop_only = FALSE,
pre_fixef = est_parms,
gen_mlx_file = TRUE,
mlx_name = "~/pmxNODE/mlx_example1",
data_file = "~/pmxNODE/data_example1_mlx.csv",
header_types = c("id","time","amount","observation"))
run_mlx("~/pmxNODE/mlx_example1_ind.mlxtran")We can check now the predictions from the NODE model.
ggplot(predictions) +
geom_point(aes(x = time, y = DV)) +
geom_line(aes(x = time, y = indivPred_mode), color = "blue") +
geom_line(aes(x = time, y = popPred), color = "red") +
facet_wrap(~id)
ggplot(predictions) +
geom_point(aes(x = DV,y = indivPred_mode)) +
geom_abline(slope = 1, intercept = 0)
ggplot(predictions) +
geom_point(aes(x = DV, y = indWRes_mode)) +
geom_abline(slope = 0, intercept = 0)
ggplot(predictions) +
geom_point(aes(x = time, y = indWRes_mode)) +
geom_abline(slope = 0, intercept = 0)
Now if we want to investigate what the NNs in the NODE have learned, we can plot the derivatives versus the states. This visualizes the dynamics on ODE-level identified by the NODE. To functions are available for this, either the der_state_plot_mlx or the rhs_plot_mlx functions. The first one generates the derivative vs. state plot for a single NN. The name of the NN, the minimal and maximal input to the NN, and either the estimated parameters or directly the name of the fitted Monolix file must be given. Additionally, it needs to be specified if the NN is a time-dependent NN.
der_state_plot_mlx("c", min_state = 0, max_state = 10, mlx_file = "~/pmxNODE/mlx_example1_ind.mlxtran", plot_type = "ggplot")
der_state_plot_mlx("ct", min_state = 0, max_state = 24, mlx_file = "~/pmxNODE/mlx_example1_ind.mlxtran", time_nn = TRUE, plot_type = "ggplot")
With the second one, the entire right-hand side of a differential equation can be plotted, e.g., the combination of multiple NNs or NNs combined with mechanistic parts. The right-hand side equation must be given as a string, the variable for the x-axis needs to be defined, and the inputs must be given as a dataframe with columns for each variable in the right-hand side equation. Additionally, a vector with information concerning time-dependency of the NNs in the right-hand side equation mus be provided or else all NNs are assumed to be non-time-depenedent. Note that for NNc inputs, the predictions need to be multiplied with the volume of distribution since the inputs to NNc in the model is amount and not concentration.
est_parms <- pre_fixef_extractor_mlx("~/pmxNODE/mlx_example1_ind.mlxtran")
rhs_inputs <- data.frame(id = predictions$id,
NNc = predictions$popPred*est_parms["V_pop"],
NNct = predictions$time,
dose = 10)
rhs_plot_mlx("NNc + dose * NNct", x_var = "NNc", inputs = rhs_inputs, mlx_file = "~/pmxNODE/mlx_example1_ind.mlxtran", time_nn = c(FALSE, TRUE))
rhs_plot_mlx("NNc + dose * NNct", x_var = "NNct", inputs = rhs_inputs, mlx_file = "~/pmxNODE/mlx_example1_ind.mlxtran", time_nn = c(FALSE, TRUE))
Similar plots can also be generated on individual level.
ind_der_state_plot_mlx("c", min_state = 0, max_state = 10, mlx_file = "~/pmxNODE/mlx_example1_ind.mlxtran", plot_type = "ggplot")
ind_der_state_plot_mlx("ct", min_state = 0, max_state = 24, mlx_file = "~/pmxNODE/mlx_example1_ind.mlxtran", time_nn = TRUE, plot_type = "ggplot")
ind_rhs_plot_mlx("NNc + dose * NNct", x_var = "NNc", group = "id", inputs = rhs_inputs, mlx_file = "~/pmxNODE/mlx_example1_ind.mlxtran", time_nn = c(FALSE, TRUE))
ind_rhs_plot_mlx("NNc + dose * NNct", x_var = "NNct", group = "id", inputs = rhs_inputs, mlx_file = "~/pmxNODE/mlx_example1_ind.mlxtran", time_nn = c(FALSE, TRUE))
