Recalling time in the pre block "t"

1

Hello,
I was able to run the following code few days ago

mod = @model begin
        @param   begin
            TVCL      ∈ RealDomain(lower = 0.00001)
            tvKTV     ∈ RealDomain(lower = 0.00001)
            TVVc      ∈ RealDomain(lower = 0.000001, upper = 8) 
            wt_cl     ∈ RealDomain()
            wt_v      ∈ RealDomain()
            circuit   ∈ RealDomain()
            ecmo      ∈ RealDomain()
            Ω         ∈ PDiagDomain(3)
            σ_add     ∈ RealDomain(lower = 0.0001)
            σ_prop    ∈ RealDomain(lower = 0.0001)
        end
        @random begin
            η ~ MvNormal(Ω)
        end
        @covariates WT AGE_D ISMALE SCR GFR AGEYRS TYPE_MODELING IS_BLEEDING IS_CIRCUIT_CHANGE ECMO_DAYS Occassions IS_CIRCUIT_CHANGE_UPDATE IS_BLEEDING_UPDATE
        @pre begin
            KTV = tvKTV*exp(η[1])*(ECMO_DAYS/6)^ecmo
            CL     = TVCL * exp(η[2])*(1-exp(-KTV*t))*(WT/50)^wt_cl*circuit^IS_CIRCUIT_CHANGE
            Vc     = TVVc * exp(η[3])*(WT/50)^wt_v
      
        end
        @dynamics  Central1
      
        @derived begin
            CP = @. Central / Vc
            CONC = @. Normal(abs(CP), sqrt(CP^2*σ_prop+σ_add))
        end
end

However, I am receiving the following error

ERROR: ArgumentError: Using explicit time in the pre or dcp blocks with Central1 <: ExplicitModel is not support.
Stacktrace:
 [1] macro expansion
   @ /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/MxXdQ/src/dsl/model_macro.jl:1281 [inlined]
 [2] top-level scope
   @ ~/data/code/UFH_ECMO_UTAH/UFH_UTAH/ECMO_UTAH_NEW.jl:23

Could you please help ?

2

@ahmed.salem the previous use of t was a bug that has been addressed as you can see here Release Notes for Pumas 2.2.0 · Pumas

Use a differential equation model that automatically gets covered to a linear ODE. And in that you can use the model as is.

3

I see, thanks … I adjusted the model accordingly:

mod = @model begin
        @param   begin
            TVCL      ∈ RealDomain(lower = 0.0001)
            tvKTV     ∈ RealDomain(lower = 0.0001)
            TVVc      ∈ RealDomain(lower  = 0.0001) 
            wt_cl     ∈ RealDomain(upper = 2)
            wt_v      ∈ RealDomain(upper = 2)
            circuit   ∈ RealDomain(lower = 0.0001)
            ecmo      ∈ RealDomain()
            Ω         ∈ PDiagDomain(3)
            σ_add     ∈ RealDomain(lower = 0.00001)
            σ_prop    ∈ RealDomain(lower = 0.00001)
        end
        @random begin
            η ~ MvNormal(Ω)
        end
        @covariates WT AGE_D ISMALE SCR GFR AGEYRS TYPE_MODELING IS_BLEEDING IS_CIRCUIT_CHANGE ECMO_DAYS Occassions IS_CIRCUIT_CHANGE_UPDATE IS_BLEEDING_UPDATE
        @pre begin
            KTV = tvKTV*exp(η[1])*(ECMO_DAYS/6)^ecmo
            CL     = TVCL * exp(η[2])*(1-exp(-KTV*(t+ 1e-10)))*(WT/50)^wt_cl*circuit^IS_CIRCUIT_CHANGE
            Vc     = TVVc * exp(η[3])*(WT/50)^wt_v
      
        end
        @dynamics  begin
                Central'  = -(CL/Vc) * Central
        end 
      
        @derived begin
            CP = @. Central / Vc
            CONC = @. Normal(abs(CP), sqrt(CP^2*σ_prop+σ_add))
        end
end

However I am getting the following warning message and the fit function does not allow estimation

modfit = fit(mod, estimation, param, Pumas.FOCEI())
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
Iter     Function value   Gradient norm 
     0              Inf              NaN
 * time: 4.100799560546875e-5
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
FittedPumasModel

Successful minimization:                     false

Likelihood approximation:               Pumas.FOCE
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
Log-likelihood value:                         -Inf
Number of subjects:                            159
Number of parameters:         Fixed      Optimized
                                  0             12
Observation records:         Active        Missing

I tried even to initialize the ODE at very small time 1e-10 but I am not able to get it work in Pumas, although I did the same run in another software and it ran successfully with the same results as when I used to use the closed form solution. So, I do not think there is a problem with model specification. I am wondering if this problem is due to the settings of the ODE solver? please let me know

4

Did you try findinfluential? Maybe there’s an issue with one of the covariates you use.

5

Thanks for the reply. I tried to use this function. It seems that ID = 20 is the problem

modfit = findinfluential(mod, estimation, param, Pumas.FOCEI(), k = 100)
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /builds/PumasAI/PumasSystemImages-jl/.julia/packages/SciMLBase/mndcy/src/integrator_interface.jl:345
100-element Vector{NamedTuple{(:id, :nll), Tuple{String, Float64}}}:
 (id = "20", nll = Inf)
 (id = "112", nll = 304.7092273740628)
 (id = "125", nll = 262.8198478415608)
 (id = "248", nll = 261.1547795464976)
 (id = "185", nll = 218.91950414263349)
 (id = "25", nll = 212.6464723493648)
 (id = "181", nll = 207.2863831956399)
 (id = "232", nll = 187.74712293115962)
 (id = "123", nll = 184.22096771056218)
 (id = "74", nll = 182.15866723353736)
 (id = "139", nll = 181.86417917383557)
 (id = "19", nll = 178.08662883120795)
 (id = "192", nll = 175.88409508526817)
 (id = "131", nll = 175.62432969437899)
 (id = "200", nll = 165.37169632282868)

I had a closer look at this subject’s data

julia> @chain df begin
               @rsubset :ID ∈ [20]
               select(_,:ID,:TIME,:evid,:mdv, :AMT,:RATE,:cmt,:WT,:ECMO_DAYS,:IS_CIRCUIT_CHANGE)
       end
15×10 DataFrame
 Row │ ID     TIME      evid   mdv    AMT         RATE        cmt    WT       ECMO_DAYS  IS_CIRCUIT_CHANGE 
     │ Int64  Float64   Int64  Int64  Float64?    Float64?    Int64  Float64  Int64      Int64             
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │    20   0.0          1      1      121.38       71.4       1     3.57          4                  0
   2 │    20   1.36667      0      0  missing     missing         1     3.57          4                  0
   3 │    20   1.7          1      1      803.25       89.25      1     3.57          4                  0
   4 │    20  10.7          1      1      535.5       107.1       1     3.57          4                  0
   5 │    20  15.7          1      1     3727.08      128.52      1     3.57          4                  0
   6 │    20  17.6167       0      0  missing     missing         1     3.57          4                  0
   7 │    20  17.7333       0      0  missing     missing         1     3.57          4                  0
   8 │    20  23.95         0      0  missing     missing         1     3.57          4                  0
   9 │    20  25.9          0      0  missing     missing         1     3.57          4                  0
  10 │    20  41.8667       0      0  missing     missing         1     3.57          4                  0
  11 │    20  43.95         0      0  missing     missing         1     3.57          4                  0
  12 │    20  44.7          1      1      835.38      139.23      1     3.57          4                  0
  13 │    20  49.9833       0      0  missing     missing         1     3.57          4                  0
  14 │    20  50.7          1      1     3070.2       153.51      1     3.57          4                  0
  15 │    20  65.7833       0      0  missing     missing         1     3.57          4                  0

I am not able to identify any red flag. Also, as I mentioned previously, I did not experience this issue when using the closed form solution or when I fitted the data using another software. Also, I get the same warning even when running simulations.
Could this issue happen since there are close times like 1.36 and 1.7 when using ODE? thanks

6

It’s a bit special that the infusions are back-to-back. We’ll try to see if we can reproduce the issue with a similar regimen.

In the other software, do you happen to know how the ODE was solved? Numerical integration or matrix exponential?

7

t is in an exp function call in the pre block, I don’t think any software would be able to use matrix exponentials, would it?

8

In NONMEM, you can refer to time in the PK block, but it will be considered piece-wise constant between events which allows for using the matrix exponential. Other software packages might use a similar interpretation.

2 Likes
9

I am not sure about the ODE solver type. But one thing to note is that the other software has only one option to use SAEM algorithm but I used FOCEI in Pumas but still this issue happen even when I do simulations so I think it is due to the ODE solver when including this subject in the data.

10

I was able to recreate locally. Let me try to figure out what is required for this to run. The issue is related to the fact that infusions stop just as another one starts. I’ll return when I have further information :slight_smile:

1 Like
11

Okay, so what you’re hitting is the age old problem of computers having to represent numbers. Long story short, you never have full precession in calculations on a computer, so for example 0.2+0.1 ends up not being 0.3

julia> 0.2+0.1
0.30000000000000004

Something similar is happening here. You can mitigate it by adding an ever so slightly small number to your rates, such that the infusion ends a very short time earlier. For example you can add 1e-8 to your subject’s rates, and it works.

df=DataFrame(:ID=>fill(20, 15), :dv=>fill(missing,15), :time=>[0.0,1.36667,1.7,10.7,15.7,17.6167,17.7333,23.95,25.9,41.8667,43.95,44.7,49.9833,50.7,65.7833], :evid=>[1,0,1,1,1,0,0,0,0,0,0,1,0,1,0],:AMT=>[121.38,missing, 803.25,535.5,3727.08,missing,missing,missing,missing,missing,missing,835.38,missing,3070.2,missing], :RATE=>[71.4+1e-8,missing,89.25+1e-8,107.1+1e-8,128.52+1e-8,missing,missing,missing,missing,missing,missing,139.23+1e-8,missing,153.51+1e-8,missing], :cmt=>fill(1,15),:WT=>fill(3.57,15), :ECMO_DAYS=>fill(4,15),:IS_CIRCUIT_CHANGE=>fill(0,15))

sub =  read_pumas(ans, id=:ID, time=:time, evid=:evid, amt=:AMT, rate=:RATE, cmt=:cmt, covariates=[:WT, :ECMO_DAYS, :IS_CIRCUIT_CHANGE])


mod = @model begin
        @param   begin
            TVCL      ∈ RealDomain(lower = 0.0001)
            tvKTV     ∈ RealDomain(lower = 0.0001)
            TVVc      ∈ RealDomain(lower  = 0.0001) 
            wt_cl     ∈ RealDomain(upper = 2)
            wt_v      ∈ RealDomain(upper = 2)
            circuit   ∈ RealDomain(lower = 0.0001)
            ecmo      ∈ RealDomain()
            Ω         ∈ PDiagDomain(3)
            σ_add     ∈ RealDomain(lower = 0.00001)
            σ_prop    ∈ RealDomain(lower = 0.00001)
        end
        @random begin
            η ~ MvNormal(Ω)
        end
        @covariates WT IS_CIRCUIT_CHANGE ECMO_DAYS
        @pre begin
            KTV = tvKTV*exp(η[1])*(ECMO_DAYS/6)^ecmo
            CL     = TVCL * exp(η[2])*(1-exp(-KTV*(t)))*(WT/50)^wt_cl*circuit^IS_CIRCUIT_CHANGE
            Vc     = TVVc * exp(η[3])*(WT/50)^wt_v
      
        end
        @dynamics  begin
                Central'  = -(CL/Vc) * Central
        end 
      
        @derived begin
            CP = @. Central / Vc
            CONC = @. Normal(abs(CP), sqrt(CP^2*σ_prop+σ_add))
        end
end

julia> simobs(mod, sub, init_params(mod))

Hope this helps. Please ask if you need a clarification.

2 Likes
12

I see, I tried yesterday to round the times also and I was able to get the estimation process to work including the whole dataset. So, I believe that what you are proposing makes sense. Also, I will try to add these little values to the infusions and will let you know. Thanks so much!

1 Like