Re: [Help-mcsim] Simple PBPK Model Help- Integration/Convergence Issue

[fredomatic]

Hi!

Yes, there can be a subtle interaction between the times at which the
integrator is required to give outputs and the time step. The
integrator does not overshoot an output time or discontinuity (to
precisely avoid stepping over discontinuities). So the output timing
can limit the last integration step taken.

However, that interaction does not matter at all if the system is
stable. The fact is that your system is not stable, because in dynamics,
you are using some output variables in differentials before updating
them with the state variables (e.g. you compute
Cnv_rob AFTER using it to get dt(Arob_n). This is dangerous because
output variables are not dealt with by the integrator (only state
variables are managed by it) and their value when entering CalcDeriv may
be off trajectory (because of a too long attempted time step).
It is your responsibility to update them properly before use. I will
make it clear in the user's manual.
See the revised model below.

Note that may want to check mass balance when developing the model (your
differential equations are not symmetrical and I added a
Q_total output for check). The dynamics of exchange between plasma and
rest of the body are very fast, and that does not help stability.

On my machine, the revised model works fine a outputs intervals ranging
from 0.0005 to 8.

I hope this helps.

Best regards.

F. Bois



#------------------------------------------------------------------------------
# RA Vreeland - 2017 - US FDA /NCTR
# PBPK model of nitrate and nitrite in humans
# Compartments for nitrate : plasma, thyroid in diffusion limited, rest
of body, saliva, GI tract(not including large intestines)
# Compartment for nitrite: plasma, saliva, GI, rest of body
#------------------------------------------------------------------------------
#
#
# n = nitrate, rf = nitrite ('reduced form'); g = (stomach, liver, small
intestines)
# s = saliva; t = thyroid; a = arterial; v = venus; Hb = hemoglobin;
metHb = methemoglobin
# rob = rest of body;
#
# As far as I can tell, the model is working, and its a software
limiting thing

States = {Aplasma_n, Arob_n, AIVD_n};

Outputs = {Cn_a, Cnv_total, Cn_rob, Cnv_rob, TestPrint, time, Q_total};

Inputs = {BW, IV_n_g, IVDose_n};


      MW_n;                   # g/mol or mg/mmol
      t_len;
      IV_n_mol;               # mol
      DosePerHour;
      P_nrob;

      volF_plasma;
      volF_rob_n;

      vol_plasma;
      vol_rob_n;

      Q_cardiac;
      Q_rob_n;

      QF_cardiac;
      QF_rob_n;

Initialize {

      MW_n = 62.004;                      # g/mol or mg/mmol
      t_len = 1.0;
      IV_n_mol = IV_n_g/MW_n;             # g * (mol/g)
      DosePerHour = IV_n_mol/t_len;

      P_nrob = 1.0;

      volF_plasma = 0.044;
      volF_rob_n = 0.956;

      vol_plasma = volF_plasma * BW;
      vol_rob_n = volF_rob_n * BW;

      QF_cardiac = 15.0;
      QF_rob_n = 1.0;   # 1 - (sum of other non plasma compartments) -->
only rob compartment currently

      Q_cardiac = QF_cardiac * pow(BW, 0.75);
      Q_rob_n = QF_rob_n * Q_cardiac;
      }; # End Initialize

Dynamics {

      # compute those before using them:
      Cn_a = Aplasma_n / vol_plasma;
      Cn_rob = Arob_n / vol_rob_n;
    # Cnv_rob = Arob_n / (vol_rob_n * P_nrob);
      Cnv_rob = Cn_rob / P_nrob;
    # Cnv_total = (Q_rob_n * Cnv_rob)/Q_cardiac ;

# Dose #
      dt(AIVD_n) = IVDose_n;
      TestPrint = dt(AIVD_n);

# Plasma #

      dt(Aplasma_n) = Q_cardiac * (Cnv_total - Cn_a ) + IVDose_n;

      dt(Arob_n) = Q_rob_n * (Cn_a - Cnv_rob);

      Q_total = Aplasma_n + Arob_n; # checking mass balance, should be =
AIVD_n

      time = t;

}; # End Dynamics

End
=======================================================================================



On 23/06/17 23:29, Vreeland, Ryan * wrote:
> Hello,
> 
>  
> 
> I have been having trouble the last few days with my PBPK model. I’ve
> been running a .sim file. The error that keeps popping up is
> 
>  
> 
> Doing analysis - 1 normal experiment
> 
>  1lsodes-- at current t (=r1), mxstep (=i1) steps  
> 
>          taken on this call before reaching tout    
> 
>  in above message, i1=       500
> 
> in above message, r1=      0.1039293261229
> 
> Warning: Integration failed - No output generated
> 
>  
> 
> I have simplified my model down to two compartments to try and figure
> this error out. A colleague and I kept looking over the code to see if
> there is anything wrong, but we do not see anything wrong.  Some
> observations I have made when attempting to figure this out is, using
> smaller doses, or different body weight seems to help and allow the
> variables to print to 24 hours.  I have also tried changing the
> tolerances on the integration function, seemed to help in some cases.
> One thing that bugs me, I have been using PrintStep to print out the
> variables, and for the same tolerances, bodyweight, and dose, using a
> <time-step> of 1 throws the error, while using 0.1 or 0.01 works just
> fine. I had thought that <time-step> is only for telling what time the
> system prints out the variables.  
> 
> I feel like it is being really finicky about converging with certain
> numbers, maybe it can’t use very big numbers.
> 
> I am hoping someone has seen this error before and or can help me figure
> out what is wrong with the file. Thank you all for your time.
> 
>  
> 
> Here is my model file.
> 
>  
> 
> #------------------------------------------------------------------------------
> 
> # RA _Vreeland_ - 2017 - US FDA /NCTR
> 
> # PBPK model of _nitrate_ and _nitrite_ in humans
> 
> # Compartments for _nitrate_ : plasma, _thyroid_ in diffusion limited,
> rest of body, saliva, GI tract(not including large intestines)
> 
> # Compartment for _nitrite_: plasma, saliva, GI, rest of body
> 
> #------------------------------------------------------------------------------
> 
> #
> 
> #
> 
> # n = _nitrate_, _rf_ = _nitrite_ ('reduced form'); g = (stomach, liver,
> small intestines)
> 
> # s = saliva; t = _thyroid_; a = arterial; v = _venus_; _Hb_ =
> _hemoglobin_; metHb = _methemoglobin_
> 
> # rob = rest of body;
> 
> #
> 
> # As far as I can tell, the model is working, and its a software
> limiting thing
> 
>  
> 
> States = {Aplasma_n, Arob_n, AIVD_n};
> 
>  
> 
> Outputs = {Cn_a, Cnv_total, Cn_rob, Cnv_rob, TestPrint, time};
> 
>  
> 
> Inputs = {BW, IV_n_g, IVDose_n};
> 
>  
> 
>  
> 
>       MW_n;                   # g/_mol_ or _mg_/_mmol_
> 
>       t_len;
> 
>       IV_n_mol;               # _mol_
> 
>       DosePerHour;
> 
>       P_nrob;
> 
>      
> 
>       volF_plasma;           
> 
>       volF_rob_n;
> 
>      
> 
>       vol_plasma;
> 
>       vol_rob_n;
> 
>      
> 
>       Q_cardiac;
> 
>       Q_rob_n;
> 
>            
> 
>       QF_cardiac;            
> 
>       QF_rob_n;              
> 
>  
> 
> Initialize {
> 
>  
> 
>       MW_n = 62.004;                          # g/_mol_ or _mg_/_mmol_
> 
>       t_len = 1.0;     
> 
>       IV_n_mol = IV_n_g/MW_n;             # g * (_mol_/g)
> 
>       DosePerHour = IV_n_mol/t_len;
> 
>      
> 
>       P_nrob = 1.0;
> 
>      
> 
>       volF_plasma = 0.044;         
> 
>       volF_rob_n = 0.956;          
> 
>      
> 
>       vol_plasma = volF_plasma * BW;        
> 
>       vol_rob_n = volF_rob_n * BW;             
> 
>  
> 
>      
> 
>       QF_cardiac = 15.0;                 
> 
>       QF_rob_n = 1.0;                     # 1 - (sum of other non plasma
> compartments) --> only rob compartment currently
> 
>      
> 
>       Q_cardiac = QF_cardiac * _pow_(BW, 0.75);
> 
>       Q_rob_n = QF_rob_n * Q_cardiac;    
> 
>                  
> 
>       }; # End Initialize
> 
> Dynamics {
> 
>  
> 
>       Cn_a = Aplasma_n/vol_plasma;
> 
>       Cn_rob = Arob_n/vol_rob_n;
> 
>      
> 
> # Dose #
> 
>       _dt_(AIVD_n) = IVDose_n;
> 
>       TestPrint = _dt_(AIVD_n);
> 
>      
> 
> # Plasma # 
> 
>      
> 
>       _dt_(Aplasma_n) = Q_cardiac * (Cnv_total -Cn_a ) + (_dt_(AIVD_n));   
> 
>     _dt_(Arob_n) = Q_rob_n * (Cn_a - Cnv_rob);
> 
>    
> 
>     Cnv_total = (Q_rob_n * Cnv_rob)/Q_cardiac ;
> 
>       Cnv_rob = Arob_n/(vol_rob_n*P_nrob);
> 
>      
> 
>      
> 
>       time = t;
> 
> }; # End Dynamics
> 
>  
> 
> End
> 
>  
> 
> This is the output file for reference
> 
>  
> 
> OutputFile("NuclearOption.out");
> 
>   Integrate(_Lsodes_, 1e-6, 1e-12, 1);
> 
>   #Integrate(_Lsodes_, 1e-5, 1e-7, 1); #default settings of integration
> 
>   
> 
>  # 3 _mg_/_kg_ dose with _Wagners_ study. Not even close to that, right now
> 
>  
> 
> Simulation {
> 
>       IV_n_g = 0.06196;  # So doing 61.96 dose with BW 75 _kg_ and 1e-6,
> 1e-12, 1 ; print step of 0.1
> 
>       IVDose_n = PerDose(DosePerHour, 24.0, 0, 1);   
> 
>       BW = 75.0;
> 
>       # Seems like the time on the PrintStep step is when each of the
> calculations occur, not just when it prints out
> 
>       # Well, using print, and 60 _mg_ dose with 1e-6, 1e-12 causes it
> to fail, but using print step with printing every 0.1 steps it works, go
> bloody figure
> 
>      
> 
>       PrintStep(IVDose_n, AIVD_n, TestPrint, Cn_a, Cnv_total, Cn_rob,
> Cnv_rob, Arob_n, Aplasma_n, 0.0, 24.0, 0.1 );
> 
>       #PrintStep( time, 0.0, 24, 1);
> 
>      
> 
>     }
> 
>  
> 
> END
> 
>  
> 
>  
> 
> 
> 
> _______________________________________________
> Help-mcsim mailing list
> address@hidden
> https://lists.gnu.org/mailman/listinfo/help-mcsim
>