cannot simulate large model

Hi Thank you for your work!

I am trying to simulate large models with mrgsolve. Some models are simulated OK, the another ones result in errors and stop simulation at initial step. For example I am trying to start the model in attachment with 103 ODEs. default.cpp.txt The code

mod_th1 <- mread('default', "Y:\\some_folder")

mod_th1 %>%
  mrgsim(delta=0.1, maxsteps=5e4, hmax=0.001, hmin=0) %>%
  plot()

result in

DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step  
     (H = step size). Solver will continue anyway.
In above message, R1 =
[1] 5.444999e-243 2.469979e+225
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step  
     (H = step size). Solver will continue anyway.
In above message, R1 =
[1] 5.444999e-243 2.469979e+225
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step  
     (H = step size). Solver will continue anyway.
In above message, R1 =
[1] 5.444999e-243 2.469979e+225
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step  
     (H = step size). Solver will continue anyway.
In above message, R1 =
[1] 5.444999e-243 2.469979e+225
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step  
     (H = step size). Solver will continue anyway.
In above message, R1 =
[1] 5.444999e-243 2.469979e+225
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step  
     (H = step size). Solver will continue anyway.
In above message, R1 =
[1] 5.444999e-243 2.469979e+225
DLSODA-  Warning..Internal T (=R1) and H (=R2) are
      such that in the machine, T + H = T on the next step  
     (H = step size). Solver will continue anyway.
In above message, R1 =
[1] 5.444999e-243 2.469979e+225
...

I tried to simulate the same model using RxODE in R and it was successful. BTW RxODE uses LSODA which is (as I understand) the same algorithm as in mrgsolve. I'v used the same options for simulations

# RxODE
solution <- model_default$solve(
    theta_default, 
    events=qd_default_task, 
    init=init_default, 
    maxsteps = 50000L, 
    hmax=0.001,
    hmin=0, 
    returnType="data.frame"
)

Can you please help me with that? Are there any hidden settings to make it work in mrgsolve? Are there some limitations for model complexity and size?

The similar behavior can be seen for another models from my collection.

Hi @metelkin

Thanks for posting the model. I don't believe there is any theoretical maximum size for the model. I do think the solvers are the same and you've identified the different settings that I would have tried as well.

When I run the model I get the following. It looks like some compartments are going to some extreme values.

I'm wondering if you have the SBML for this? Assuming that you did this programmatically. I could try to run it through a translator that I have. I can't say anything is wrong exactly, but that's the place where I would start.

I would be happy to work on this with you either here or in a private repository if you would like.

Best Regards, Kyle


library(mrgsolve)

mod <- mread("~/demo/default_cpp.txt")
#> Compiling default_cpp_txt ... done.

mrgsim(mod, end = 4, delta = 1, mxhnil = 1)
#> DLSODA-  Warning..Internal T (=R1) and H (=R2) are
#>       such that in the machine, T + H = T on the next step  
#>      (H = step size). Solver will continue anyway.
#> In above message, R1 =
#> [1]  8.603459e-74 7.141320e-218
#> DLSODA-  Above warning has been issued I1 times.  
#>      It will not be issued again for this problem.
#> In above message, I = 
#> [1] 1
#> Model:  default_cpp_txt 
#> Dim:    5 x 105 
#> Time:   0 to 4 
#> ID:     1 
#>      ID time Th1_pl_amt Th1_ln_amt Th1_nld_amt Th1_ld_amt Th1_nle_amt
#> [1,]  1    0     351675      0.000           0          0           0
#> [2,]  1    1     315212      0.000       17612      17612           0
#> [3,]  1    2     281301      4.756       32762      32762        1191
#> [4,]  1    3     251468     16.697       44082      44082        4182
#> [5,]  1    4     224320     32.711       53043      53043        8194
#>      Th1_le_amt IL2_ln_amt IL2_pl_amt IL2_nld_amt IL2_ld_amt IL2_nle_amt
#> [1,]          0  0.000e+00      1.361    0.000000   0.000000   0.0000000
#> [2,]          0  0.000e+00      1.501    0.000000   0.000000   0.0000000
#> [3,]       1191  0.000e+00      1.632    0.004756   0.004756   0.0000000
#> [4,]       4182  1.542e-06      1.747    0.016712   0.016712   0.0003861
#> [5,]       8194  6.562e-06      1.852    0.032772   0.032772   0.0016434
#>      IL2_le_amt IL3_ln_amt IL3_pl_amt IL3_nld_amt IL3_ld_amt IL3_nle_amt
#> [1,]  0.0000000  0.000e+00     0.0000    0.000000   0.000000   0.0000000
#> [2,]  0.0000000  0.000e+00     0.1407    0.000000   0.000000   0.0000000
#> [3,]  0.0000000  0.000e+00     0.2715    0.004756   0.004756   0.0000000
#> [4,]  0.0003861  1.542e-06     0.3866    0.016712   0.016712   0.0003861
#> [5,]  0.0016434  6.562e-06     0.4913    0.032772   0.032772   0.0016434
#>      IL3_le_amt IL10_ln_amt IL10_pl_amt IL10_nld_amt IL10_ld_amt
#> [1,]  0.0000000   0.000e+00       3.585     0.000000    0.000000
#> [2,]  0.0000000   0.000e+00       3.725     0.000000    0.000000
#> [3,]  0.0000000   0.000e+00       3.856     0.004756    0.004756
#> [4,]  0.0003861   1.542e-06       3.971     0.016712    0.016712
#> [5,]  0.0016434   6.562e-06       4.076     0.032772    0.032772
#>      IL10_nle_amt IL10_le_amt IL21_ln_amt IL21_pl_amt IL21_nld_amt
#> [1,]    0.0000000   0.0000000   0.000e+00      0.0000     0.000000
#> [2,]    0.0000000   0.0000000   0.000e+00      0.1407     0.000000
#> [3,]    0.0000000   0.0000000   0.000e+00      0.2715     0.004756
#> [4,]    0.0003861   0.0003861   1.542e-06      0.3866     0.016712
#> [5,]    0.0016434   0.0016434   6.562e-06      0.4913     0.032772
#>      IL21_ld_amt IL21_nle_amt IL21_le_amt IFNg_ln_amt IFNg_pl_amt
#> [1,]    0.000000    0.0000000   0.0000000    0.000000        0.00
#> [2,]    0.000000    0.0000000   0.0000000    0.000000       45.44
#> [3,]    0.004756    0.0000000   0.0000000    0.000000       87.70
#> [4,]    0.016712    0.0003861   0.0003861    0.000498      124.88
#> [5,]    0.032772    0.0016434   0.0016434    0.002120      158.71
#>      IFNg_nld_amt IFNg_ld_amt IFNg_nle_amt IFNg_le_amt GMCSF_ln_amt
#> [1,]        0.000       0.000       0.0000      0.0000    0.000e+00
#> [2,]        0.000       0.000       0.0000      0.0000    0.000e+00
#> [3,]        1.536       1.536       0.0000      0.0000    0.000e+00
#> [4,]        5.398       5.398       0.1247      0.1247    1.927e-07
#> [5,]       10.586      10.586       0.5309      0.5309    8.203e-07
#>      GMCSF_pl_amt GMCSF_nld_amt GMCSF_ld_amt GMCSF_nle_amt GMCSF_le_amt
#> [1,]        1.033     0.0000000    0.0000000     0.000e+00    0.000e+00
#> [2,]        1.051     0.0000000    0.0000000     0.000e+00    0.000e+00
#> [3,]        1.067     0.0005946    0.0005946     0.000e+00    0.000e+00
#> [4,]        1.081     0.0020890    0.0020890     4.826e-05    4.826e-05
#> [5,]        1.094     0.0040965    0.0040965     2.054e-04    2.054e-04
#>      TNFa_ln_amt TNFa_pl_amt TNFa_nld_amt TNFa_ld_amt TNFa_nle_amt
#> [1,]   0.000e+00       1.425       0.0000      0.0000      0.00000
#> [2,]   0.000e+00       8.853       0.0000      0.0000      0.00000
#> [3,]   0.000e+00      15.762       0.2512      0.2512      0.00000
#> [4,]   8.141e-05      21.839       0.8825      0.8825      0.02039
#> [5,]   3.465e-04      27.370       1.7306      1.7306      0.08678
#>      TNFa_le_amt MHC2_th0_lca_ln_is_amt TCRc_MHC2_th0_lca_ln_is_amt
#> [1,]     0.00000              3.994e-07                   0.000e+00
#> [2,]     0.00000             -1.534e-03                   1.534e-03
#> [3,]     0.00000             -1.004e+05                   1.004e+05
#> [4,]     0.02039             -5.168e+20                   5.168e+20
#> [5,]     0.08678             -1.102e+52                   1.102e+52
#>      TCRc_th0_lca_ln_is_amt CD80_th0_lca_ln_is_amt
#> [1,]              6.075e-08              4.622e-08
#> [2,]             -1.534e-03             -6.256e-03
#> [3,]             -1.004e+05             -4.918e+08
#> [4,]             -5.168e+20             -5.998e+30
#> [5,]             -1.102e+52             -1.374e+75
#>      CD28_CD80_th0_lca_ln_is_amt CD28_th0_lca_ln_is_amt
#> [1,]                   0.000e+00              1.417e-08
#> [2,]                   6.256e-03             -1.219e-02
#> [3,]                   4.918e+08             -1.575e+09
#> [4,]                   5.998e+30             -3.673e+31
#> [5,]                   1.374e+75             -1.774e+76
#>      CD86_th0_lca_ln_is_amt CD28_CD86_th0_lca_ln_is_amt
#> [1,]              1.884e-08                   0.000e+00
#> [2,]             -5.929e-03                   5.929e-03
#> [3,]             -1.084e+09                   1.084e+09
#> [4,]             -3.073e+31                   3.073e+31
#> [5,]             -1.636e+76                   1.636e+76
#>      OX40L_th0_lca_ln_is_amt OX40_OX40L_th0_lca_ln_is_amt
#> [1,]               2.763e-09                    0.000e+00
#> [2,]              -3.663e-05                    3.663e-05
#> [3,]              -2.466e+03                    2.466e+03
#> [4,]              -1.342e+19                    1.342e+19
#> [5,]              -3.201e+50                    3.201e+50
#>      OX40_th0_lca_ln_is_amt MHC2_th0_ideca_ln_is_amt
#> [1,]              4.868e-09                1.884e-07
#> [2,]             -3.662e-05               -7.235e-04
#> [3,]             -2.466e+03               -2.235e+04
#> [4,]             -1.342e+19               -2.558e+19
#> [5,]             -3.201e+50               -2.701e+49
#>      TCRc_MHC2_th0_ideca_ln_is_amt TCRc_th0_ideca_ln_is_amt
#> [1,]                     0.000e+00                6.075e-08
#> [2,]                     7.236e-04               -7.236e-04
#> [3,]                     2.235e+04               -2.235e+04
#> [4,]                     2.558e+19               -2.558e+19
#> [5,]                     2.701e+49               -2.701e+49
#>      CD80_th0_ideca_ln_is_amt CD28_CD80_th0_ideca_ln_is_amt
#> [1,]                4.672e-09                     0.000e+00
#> [2,]               -6.324e-04                     6.324e-04
#> [3,]               -3.332e+06                     3.332e+06
#> [4,]               -1.441e+26                     1.441e+26
#> [5,]               -3.335e+65                     3.335e+65
#>      CD28_th0_ideca_ln_is_amt CD86_th0_ideca_ln_is_amt
#> [1,]                1.417e-08                5.853e-10
#> [2,]               -8.166e-04               -1.842e-04
#> [3,]               -5.588e+06               -2.256e+06
#> [4,]               -3.710e+26               -2.269e+26
#> [5,]               -1.554e+66               -1.221e+66
#>      CD28_CD86_th0_ideca_ln_is_amt OX40L_th0_ideca_ln_is_amt
#> [1,]                     0.000e+00                 2.763e-09
#> [2,]                     1.842e-04                -3.663e-05
#> [3,]                     2.256e+06                -2.466e+03
#> [4,]                     2.269e+26                -1.342e+19
#> [5,]                     1.221e+66                -3.201e+50
#>      OX40_OX40L_th0_ideca_ln_is_amt OX40_th0_ideca_ln_is_amt
#> [1,]                      0.000e+00                4.868e-09
#> [2,]                      3.663e-05               -3.662e-05
#> [3,]                      2.466e+03               -2.466e+03
#> [4,]                      1.342e+19               -1.342e+19
#> [5,]                      3.201e+50               -3.201e+50
#>      MHC2_th0_bi_ln_is_amt TCRc_MHC2_th0_bi_ln_is_amt
#> [1,]             3.165e-07                  0.000e+00
#> [2,]            -1.215e-03                  1.216e-03
#> [3,]            -6.307e+04                  6.307e+04
#> [4,]            -2.038e+20                  2.038e+20
#> [5,]            -1.714e+51                  1.714e+51
#>      TCRc_th0_bi_ln_is_amt CD80_th0_bi_ln_is_amt
#> [1,]             6.075e-08             1.751e-10
#> [2,]            -1.216e-03            -2.370e-05
#> [3,]            -6.307e+04            -5.970e+04
#> [4,]            -2.038e+20            -1.021e+24
#> [5,]            -1.714e+51            -5.500e+62
#>      CD28_CD80_th0_bi_ln_is_amt CD28_th0_bi_ln_is_amt
#> [1,]                  0.000e+00             1.417e-08
#> [2,]                  2.370e-05            -3.905e-04
#> [3,]                  5.970e+04            -2.208e+06
#> [4,]                  1.021e+24            -8.640e+25
#> [5,]                  5.500e+62            -1.075e+65
#>      CD86_th0_bi_ln_is_amt CD28_CD86_th0_bi_ln_is_amt
#> [1,]             1.166e-09                  0.000e+00
#> [2,]            -3.668e-04                  3.668e-04
#> [3,]            -2.148e+06                  2.148e+06
#> [4,]            -8.538e+25                  8.538e+25
#> [5,]            -1.070e+65                  1.070e+65
#>      OX40L_th0_bi_ln_is_amt OX40_OX40L_th0_bi_ln_is_amt
#> [1,]              3.894e-09                   0.000e+00
#> [2,]             -5.161e-05                   5.161e-05
#> [3,]             -4.897e+03                   4.897e+03
#> [4,]             -5.290e+19                   5.290e+19
#> [5,]             -4.975e+51                   4.975e+51
#>      OX40_th0_bi_ln_is_amt MHC2_th0_m1_ld_is_amt
#> [1,]             4.868e-09             1.417e-07
#> [2,]            -5.161e-05            -5.440e-04
#> [3,]            -4.897e+03            -1.264e+04
#> [4,]            -5.290e+19            -8.181e+18
#> [5,]            -4.975e+51            -2.762e+48
#>      TCRc_MHC2_th0_m1_ld_is_amt TCRc_th0_m1_ld_is_amt
#> [1,]                  0.000e+00             6.075e-08
#> [2,]                  5.442e-04            -5.441e-04
#> [3,]                  1.264e+04            -1.264e+04
#> [4,]                  8.181e+18            -8.181e+18
#> [5,]                  2.762e+48            -2.762e+48
#>      CD80_th0_m1_ld_is_amt CD28_CD80_th0_m1_ld_is_amt
#> [1,]             9.144e-09                  0.000e+00
#> [2,]            -1.238e-03                  1.238e-03
#> [3,]            -5.596e+07                  5.596e+07
#> [4,]            -2.870e+29                  2.870e+29
#> [5,]            -1.347e+73                  1.347e+73
#>      CD28_th0_m1_ld_is_amt CD86_th0_m1_ld_is_amt
#> [1,]             1.417e-08             1.834e-08
#> [2,]            -7.009e-03            -5.771e-03
#> [3,]            -6.626e+08            -6.066e+08
#> [4,]            -7.522e+30            -7.235e+30
#> [5,]            -8.026e+74            -7.892e+74
#>      CD28_CD86_th0_m1_ld_is_amt OX40L_th0_m1_ld_is_amt
#> [1,]                  0.000e+00              2.763e-09
#> [2,]                  5.771e-03             -3.663e-05
#> [3,]                  6.066e+08             -2.466e+03
#> [4,]                  7.235e+30             -1.342e+19
#> [5,]                  7.892e+74             -3.201e+50
#>      OX40_OX40L_th0_m1_ld_is_amt OX40_th0_m1_ld_is_amt
#> [1,]                   0.000e+00             4.868e-09
#> [2,]                   3.663e-05            -3.662e-05
#> [3,]                   2.466e+03            -2.466e+03
#> [4,]                   1.342e+19            -1.342e+19
#> [5,]                   3.201e+50            -3.201e+50
#>      MHC2_th0_m1_nld_is_amt TCRc_MHC2_th0_m1_nld_is_amt
#> [1,]              1.417e-07                   0.000e+00
#> [2,]             -5.440e-04                   5.442e-04
#> [3,]             -1.264e+04                   1.264e+04
#> [4,]             -8.181e+18                   8.181e+18
#> [5,]             -2.762e+48                   2.762e+48
#>      TCRc_th0_m1_nld_is_amt CD80_th0_m1_nld_is_amt
#> [1,]              6.075e-08              9.144e-09
#> [2,]             -5.441e-04             -1.238e-03
#> [3,]             -1.264e+04             -5.596e+07
#> [4,]             -8.181e+18             -2.870e+29
#> [5,]             -2.762e+48             -1.347e+73
#>      CD28_CD80_th0_m1_nld_is_amt CD28_th0_m1_nld_is_amt
#> [1,]                   0.000e+00              1.417e-08
#> [2,]                   1.238e-03             -7.009e-03
#> [3,]                   5.596e+07             -6.626e+08
#> [4,]                   2.870e+29             -7.522e+30
#> [5,]                   1.347e+73             -8.026e+74
#>      CD86_th0_m1_nld_is_amt CD28_CD86_th0_m1_nld_is_amt
#> [1,]              1.834e-08                   0.000e+00
#> [2,]             -5.771e-03                   5.771e-03
#> [3,]             -6.066e+08                   6.066e+08
#> [4,]             -7.235e+30                   7.235e+30
#> [5,]             -7.892e+74                   7.892e+74
#>      OX40L_th0_m1_nld_is_amt OX40_OX40L_th0_m1_nld_is_amt
#> [1,]               2.763e-09                    0.000e+00
#> [2,]              -3.663e-05                    3.663e-05
#> [3,]              -2.466e+03                    2.466e+03
#> [4,]              -1.342e+19                    1.342e+19
#> [5,]              -3.201e+50                    3.201e+50
#>      OX40_th0_m1_nld_is_amt
#> [1,]              4.868e-09
#> [2,]             -3.662e-05
#> [3,]             -2.466e+03
#> [4,]             -1.342e+19
#> [5,]             -3.201e+50

^{Created on 2018-11-08 by the reprex package (v0.2.1)}

Hi Kyle, Thank you for fast reply.

I was able to reproduce your simulations but it is still unclear how to simulate the whole curve from 0 to 100. Playing with solver options I just could do simulations from 0 to 60 with low resolution.

Answering to your questions:

1) You are right. I created the code programmatically, not from SBML but from the another internal format. I can produce code in several formats including SBML if it helps. th1_sbml.xml.txt

2) Possibly it will be helpful if I share some example of correct simulations Th1_ln~time in two softwares for the same model. They looks very similar. In DBSolve: 2018-11-09_13-23-30 In RxODE

3) Another unusual behavior. Working with this and another models in mrgsolve I notices that the parameter delta is very critical for results of simulations. If the parameter value is not guessed the simulations will fail even for simple models, for example produces random positive and negative numbers or produces errors like in my first post. I am not sure but as I understand the parameter should regulate output time points but not solver algorithm directly. Am I wrong?

4) Working in private repository might be more convenient but I haven't got one here on github.

I don't know for sure if the model code is ok or not. But the fact that you have compartments going off the rails suggests to me that there is a coding issue. And that the simulations work in other tools doesn't mean much toward validating your implementation in mrgsolve.

My guess is that if delta affects the results (or the probability of success for the simulation) then something is wrong in the model code. Same thing if "sometimes it works, sometimes it doesn't work". That almost always happens because the model itself implements undefined behavior.

Here's an example from pbpk model showing identical results with very different delta.


library(mrgsolve)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

withr::with_dir(Sys.getenv("devmodels"), {
  mod <- mread("rifampicin/rifampicin_midazolam.cpp", rtol = 1E-10, atol = 1E-10)
})
#> Compiling rifampicin_midazolam_cpp ...
#> done.

init(mod) 
#> 
#>  Model initial conditions (N=41):
#>  name           value . name                    value .
#>  Cadipose (6)   0     | CLIV2 (25)              0     |
#>  Cent (9)       0     | CLIV3 (26)              0     |
#>  central (3)    0     | CLIV4 (27)              0     |
#>  CHC1 (15)      0     | CLIV5 (28)              0     |
#>  CHC2 (16)      0     | Cmucblood (8)           0     |
#>  CHC3 (17)      0     | Cmuscle (4)             0     |
#>  CHC4 (18)      0     | Cportal (29)            0     |
#>  CHC5 (19)      0     | Cserosa (7)             0     |
#>  CHE1 (10)      0     | Cskin (5)               0     |
#>  CHE2 (11)      0     | CYP3A4_ratio_ent (41)   1     |
#>  CHE3 (12)      0     | CYP3A4_ratio_HC1 (36)   1     |
#>  CHE4 (13)      0     | CYP3A4_ratio_HC2 (37)   1     |
#>  CHE5 (14)      0     | CYP3A4_ratio_HC3 (38)   1     |
#>  CLIV1 (24)     0     | CYP3A4_ratio_HC4 (39)   1     |
#>  name                    value
#>  CYP3A4_ratio_HC5 (40)   1    
#>  mCadipose (23)          0    
#>  mcentral (20)           0    
#>  mCmuscle (21)           0    
#>  mCskin (22)             0    
#>  Mgutlumen (2)           0    
#>  UGT_ratio_ent (35)      1    
#>  UGT_ratio_HC1 (30)      1    
#>  UGT_ratio_HC2 (31)      1    
#>  UGT_ratio_HC3 (32)      1    
#>  UGT_ratio_HC4 (33)      1    
#>  UGT_ratio_HC5 (34)      1    
#>  Xgutlumen (1)           0    
#>  . ...                   .

rif <- ev(amt = 75, ii = 24, addl = 6)
mid <- ev(amt = 3, cmt = 2)

rif_mid <- ev_seq(rif,  wait = -12,  mid)

mid <- filter(rif_mid, cmt==2)

both <- as_data_set(mid, rif_mid)

mod %>% mrgsim(both,end = 172) %>% plot(Ccentral+mCcentral~time)


add <- c(22.1, 101.3, 0.3, 160, 165.9)
out1 <- mod %>% mrgsim(both, end = 168, delta = 0.1, add = add)
out2 <- mod %>% mrgsim(both, end = 168, delta = 168, add = add)

dim(out1)
#> [1] 3369   47
dim(out2)
#> [1] 17 47

x1 <- filter(out1, time %in% add) %>% as.data.frame
x2 <- filter(out2, time %in% add) %>% as.data.frame

x1 %>% select(ID,time,Cskin, Ccentral, mCcentral, mCmuscle) 
#>    ID  time       Cskin    Ccentral    mCcentral     mCmuscle
#> 1   1   0.3 0.000000000 0.000000000 0.0000000000 0.0000000000
#> 2   1  22.1 0.000000000 0.000000000 0.0000000000 0.0000000000
#> 3   1 101.3 0.000000000 0.000000000 0.0000000000 0.0000000000
#> 4   1 160.0 0.000000000 0.000000000 0.0023572327 0.0022268926
#> 5   1 165.9 0.000000000 0.000000000 0.0008882546 0.0007681547
#> 6   2   0.3 0.204913572 0.538096723 0.0000000000 0.0000000000
#> 7   2  22.1 0.006785672 0.002279376 0.0000000000 0.0000000000
#> 8   2 101.3 0.667192305 0.229896915 0.0000000000 0.0000000000
#> 9   2 160.0 0.024110507 0.007833348 0.0005620594 0.0005702355
#> 10  2 165.9 0.003765719 0.001227367 0.0001483272 0.0001316433
x2 %>% select(ID,time,Cskin, Ccentral, mCcentral, mCmuscle) 
#>    ID  time       Cskin    Ccentral    mCcentral     mCmuscle
#> 1   1   0.3 0.000000000 0.000000000 0.0000000000 0.0000000000
#> 2   1  22.1 0.000000000 0.000000000 0.0000000000 0.0000000000
#> 3   1 101.3 0.000000000 0.000000000 0.0000000000 0.0000000000
#> 4   1 160.0 0.000000000 0.000000000 0.0023572327 0.0022268926
#> 5   1 165.9 0.000000000 0.000000000 0.0008882546 0.0007681547
#> 6   2   0.3 0.204913572 0.538096723 0.0000000000 0.0000000000
#> 7   2  22.1 0.006785672 0.002279376 0.0000000000 0.0000000000
#> 8   2 101.3 0.667192305 0.229896915 0.0000000000 0.0000000000
#> 9   2 160.0 0.024110507 0.007833348 0.0005620594 0.0005702355
#> 10  2 165.9 0.003765719 0.001227367 0.0001483272 0.0001316433

all.equal(x1,x2)
#> [1] TRUE

^{Created on 2018-11-09 by the reprex package (v0.2.1)}

I would be interested in the RxODE code, of course, something to push the solver, altough it really doesn't relate to this project.

I would also like to mention, RxODE uses a different default solver than Kyle does. In theory they should be the same, but you should compare the same solver. The lsoda errors may be the same. You can use the same solver by method="lsoda" instead of method="liblsoda" in RxODE

Also DBSolve uses LSODE instead of LSODA

You can fiddle with the solver all you want. Every time I've seen this sort of thing, there has been an implementation issue.

#>      CD28_th0_ideca_ln_is_amt CD86_th0_ideca_ln_is_amt
#> [1,]                1.417e-08                5.853e-10
#> [2,]               -8.166e-04               -1.842e-04
#> [3,]               -5.588e+06               -2.256e+06
#> [4,]               -3.710e+26               -2.269e+26
#> [5,]               -1.554e+66               -1.221e+66

I don't know if RxODE supports pre-processor directives in the model or not. But the mrgsolve implementation here uses them. It is useless comparing tools here because the implementation is different.

You are probbly correct. The solvers are different between the tools, though. RxODE does not use pre-processor directives.

I tested the another methods for solving in RxODE for the model. All available solvers show the same result: "liblsoda", "lsoda", "dop853". No errors.

I don't know for sure if the model code is ok or not. But the fact that you have compartments going off the rails suggests to me that there is a coding issue. And that the simulations work in other tools doesn't mean much toward validating your implementation in mrgsolve.

You are right I cannot be 100% confident but I still believe the model code is OK. Why? (1) The another models generated by my code show successfully the same result in RxODE and mrgsolve. (2) There is no error messages at compilation step.

I would be interested in the RxODE code, of course, something to push the solver, altough it really doesn't relate to this project.

I can also share the model in RxODE little bit later if it will help. I need some time to clear the code.

I cannot understand mrgsolve solver code deeply but intuitively I suppose possibly something happens with initial step selection. The shared model can be characterized as the model with very fast perturbations at very initial times followed by more or less slow variable changes. It seems to me untypical for regular PK/PD models previously tested.

Matthew, I have shared the model in rxode format. The compiled model code loaded as separate file model_default.txt

I hope it will be helpful.

main.R.txt model_default.txt

Thanks metelkin.

I hope you can figure out what happened.

metrumresearchgroup / mrgsolve

cannot simulate large model #414