RxODE behaves strangely with IOV when the number of threads is higher than 1

[luyckxn]

Hi,

I have noticed something strange when updating my RxODE version from 1.0.5 to the latest available versions (first 1.0.8 then 1.0.9, the latest version available on CRAN).

In brief, when there is some treatment IOV in the model (e.g. on KA, D1 or ALAG1) and the number of threads is higher than 1, there seems to be some randomness in the results, even if there is no variability at all in the simulation.

I could reproduce my problem with a minimalist example. I have a 2-cpt PK model with IOV on the lag time and a dataset of 10 subjects. There is no variability in the simulation in the sense all ETA's and the IOV_ALAG1 have been exported into the dataset (here attached: dataset.zip ).

My code is as follows:

library(RxODE)
library(ggplot2)

model <- RxODE({
  KA=THETA_KA*exp(ETA_KA)
  CL=THETA_CL*exp(ETA_CL)
  V2=THETA_V2*exp(ETA_V2)
  V3=THETA_V3*exp(ETA_V3)
  Q=THETA_Q*exp(ETA_Q)
  S2=V2
  ALAG1=THETA_ALAG1*exp(ETA_ALAG1 + IOV_ALAG1)
  d/dt(A_DEPOT)=-KA*A_DEPOT
  d/dt(A_CENTRAL)=KA*A_DEPOT + Q*A_PERIPHERAL/V3 + (-CL/V2 - Q/V2)*A_CENTRAL
  d/dt(A_PERIPHERAL)=-Q*A_PERIPHERAL/V3 + Q*A_CENTRAL/V2
  d/dt(A_OUTPUT)=CL*A_CENTRAL/V2
  lag(A_DEPOT)=ALAG1
  CP=A_CENTRAL/S2
})

theta <- c(THETA_KA=1, THETA_CL=5, THETA_V2=80, THETA_V3=20, THETA_Q=4, THETA_ALAG1=5)

dataset <- read.csv("dataset.csv") # This dataset contains 3 boluses given at time 0, 24 & 48

RxODE::setRxThreads(6)  # Change 6 to 1 is a working workaround

results <- RxODE::rxSolve(object=model, params=theta, omega=NULL, sigma=NULL, events=dataset)

plot(results, CP)

On my laptop, when I execute this code, it gives different (and incorrect) results every time with RxODE v1.0.9, except if the number of threads is set to 1. With RxODE v1.0.5, the simulation works just fine and I get the exact same results every time.

Could you please tell me if you are able to reproduce the issue? The CSV dataset is attached to this issue.

Thanks a lot!

[mattfidler]

Thank you @luyckxn,

I will see if I can reproduce your issue.

[mattfidler]

Here is what I get with the development version:

setwd("~/src/RxODE/tests/testthat/");source("test-issue-430.R", echo=TRUE)
#> 
#> > library(RxODE)
#> RxODE 1.1.0 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> 
#> > library(ggplot2)
#> 
#> > model <- RxODE({
#> +   KA=THETA_KA*exp(ETA_KA)
#> +   CL=THETA_CL*exp(ETA_CL)
#> +   V2=THETA_V2*exp(ETA_V2)
#> +   V3=THETA_V3*exp(ETA_V3)
#> +   Q=THETA_Q*exp(E .... [TRUNCATED] 
#> 
#> > theta <- c(THETA_KA=1, THETA_CL=5, THETA_V2=80, THETA_V3=20, THETA_Q=4, THETA_ALAG1=5)
#> 
#> > dataset <- qs::qread("test-issue-430.qs") # This dataset contains 3 boluses given at time 0, 24 & 48
#> 
#> > print(head(dataset))
#>   ID TIME EVID MDV  AMT CMT RATE DOSENO  IOV_ALAG1    ETA_KA     ETA_CL
#> 1  1  0.0    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#> 2  1  0.0    1   1 1000   1    0      1 0.09510191 0.3797291 0.06294607
#> 3  1  0.5    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#> 4  1  1.0    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#> 5  1  1.5    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#> 6  1  2.0    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#>        ETA_V2    ETA_V3     ETA_Q ETA_ALAG1
#> 1 -0.02601346 0.2148261 0.1453031 0.3023562
#> 2 -0.02601346 0.2148261 0.1453031 0.3023562
#> 3 -0.02601346 0.2148261 0.1453031 0.3023562
#> 4 -0.02601346 0.2148261 0.1453031 0.3023562
#> 5 -0.02601346 0.2148261 0.1453031 0.3023562
#> 6 -0.02601346 0.2148261 0.1453031 0.3023562
#> 
#> > RxODE::setRxThreads(6)  # Change 6 to 1 is a working workaround
#> 
#> > results <- RxODE::rxSolve(object=model, params=theta, omega=NULL, sigma=NULL, events=dataset)
#> 
#> > plot(results, CP) |> print()

^{Created on 2021-06-08 by the reprex package (v2.0.0)}

[mattfidler]

For me, I get the same results with 1 thread:

setwd("~/src/RxODE/tests/testthat/");source("test-issue-430.R", echo=TRUE)
#> 
#> > library(RxODE)
#> RxODE 1.1.0 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> 
#> > library(ggplot2)
#> 
#> > model <- RxODE({
#> +   KA=THETA_KA*exp(ETA_KA)
#> +   CL=THETA_CL*exp(ETA_CL)
#> +   V2=THETA_V2*exp(ETA_V2)
#> +   V3=THETA_V3*exp(ETA_V3)
#> +   Q=THETA_Q*exp(E .... [TRUNCATED] 
#> 
#> > theta <- c(THETA_KA=1, THETA_CL=5, THETA_V2=80, THETA_V3=20, THETA_Q=4, THETA_ALAG1=5)
#> 
#> > dataset <- qs::qread("test-issue-430.qs") # This dataset contains 3 boluses given at time 0, 24 & 48
#> 
#> > print(head(dataset))
#>   ID TIME EVID MDV  AMT CMT RATE DOSENO  IOV_ALAG1    ETA_KA     ETA_CL
#> 1  1  0.0    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#> 2  1  0.0    1   1 1000   1    0      1 0.09510191 0.3797291 0.06294607
#> 3  1  0.5    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#> 4  1  1.0    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#> 5  1  1.5    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#> 6  1  2.0    0   0   NA   1    0     NA 0.09510191 0.3797291 0.06294607
#>        ETA_V2    ETA_V3     ETA_Q ETA_ALAG1
#> 1 -0.02601346 0.2148261 0.1453031 0.3023562
#> 2 -0.02601346 0.2148261 0.1453031 0.3023562
#> 3 -0.02601346 0.2148261 0.1453031 0.3023562
#> 4 -0.02601346 0.2148261 0.1453031 0.3023562
#> 5 -0.02601346 0.2148261 0.1453031 0.3023562
#> 6 -0.02601346 0.2148261 0.1453031 0.3023562
#> 
#> > RxODE::setRxThreads(1)  # Change 6 to 1 is a working workaround
#> 
#> > results <- RxODE::rxSolve(object=model, params=theta, omega=NULL, sigma=NULL, events=dataset)
#> 
#> > plot(results, CP) |> print()

^{Created on 2021-06-08 by the reprex package (v2.0.0)}

[mattfidler]

It seems like you may have some platform dependent changes based on the fix from #399

[mattfidler]

Hi @luyckxn ,

Although I couldn't reproduce your issue, it seems like the update I just pushed may fix your problem.

If you could try to confirm by updating RxODE to the development version by:

devtools::install_github("nlmixrdevelopment/RxODE")

If this does not fix the issue, I'm unsure what needs to be done to fix it. Regardless, I think the implemented lock should be done per individual so I will be keeping the fix.

[luyckxn]

Hi matt,

Thank you very much for your quick reply.

Look carefully at the 2 plots you provided (with respectively 4 threads and 1 thread).

They differ slightly... This is exactly the issue I have.
Only the second one is correct, with a single thread (I qualified the output against mrgsolve).

I will try with the development version of RxODE.

Thanks.

[luyckxn]

I installed the development version of RxODE.

The problem disappears.

Thank you very much for resolving this issue.

[mattfidler]

Ok. I can test for these differences then.

Thank you for confirming the fix.