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

[Vreeland, Ryan *]

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