What does this message mean, "CVODES returned -4..." ?

Hi, all

When I use SimBiology Analyzer, I get this message:

This would happen when I use a customized function like:

function [F]=calculateF(Drug_1,CD_1,Kc,DrugMW,CDMW)
%Drug,CD in unit of mg/mL, Kc in unit of L/mmol
Drug=Drug_1*1000/DrugMW;
CD=CD_1*1000/CDMW;
F=(-(Kc*(CD-Drug)+1)+((Kc*(CD-Drug)+1)^2+4*Kc*Drug)^(1/2))/(2*Kc*Drug);
end

This function is to calculate the free fraction of dissolved drug. This function works when I use MATLAB (not in SimBiology).

What does this message mean? I don't know what "h" is.

Answer 1

[Rick Paxson]

Hello,

First let me comment on the message you are seeing. 'h' is the time step the ODE solver, in this case CVODES, is attempting to take, and failing. It fails because it has reached a lower bound (5.99618e-17) and has not been able to satisfy error tolerances. This typically happens when a solution if blowing up and could indicate a problem with the model formulation or its initial conditions.

Things you could try to debug the problem:

1. Change the solver and / or loosen the solver tolerances and try to run the model again.
2. If you are using release 19b, download the SimBiology Debugger from here. It is easy to use and will help you figure out what is happening.

Hope this helps.

[Wei Wang]

Hello,

Thanks for your suggestions.

I have tried suggesiton 1. The result of F vs. time is tend to be a abnormal shape so I give up.

Then I am trying suggestion 2. I am not similar with using code to manipulate SimBiology model so I come up against with some problems.

After I load my model, I type this code

verify(m1)

and I get warnings and errors like this (Chinese characters in orange means "warning", in red means " error"):

Function which is related to 'rule_18' is this:

function [IntestinalAGSolubility] = calculateSolubilityCD(...
    CDConcentration,CDMW,AGMW)
% This is a function to calculate solubility of a drug,AG, according to concentration
% of cyclodextrin (CD). Equation is transformed from phase solubility test.
% AGMW: AG molecular weight; CDMW: CD molecular weight;
% CDConcentration (mg/mL); IntestinalAGSolubility (mg/mL)
% 0.4653: slope; 0.4409: intercept (mmol/L)
if (0.4653*CDConcentration*1000/CDMW+0.4409)*AGMW>2990.1
    IntestinalAGSolubility=2.9901;
elseif (0.4653*CDConcentration*1000/CDMW+0.4409)*AGMW<=2990.1
    IntestinalAGSolubility=(0.4653*CDConcentration*1000/CDMW+0.4409)*AGMW/1000;
end

I call this function in SimBiology like this:

Con0 and MW0 is used to nondimensionlize parameters so I think variables in function are dimensionless. Does error happens because I didn't use command code "sbiodose" to add dose before I verify model? My dose schedule is some complicated so I don't know how to translate them to code. Can I generate code for my simulation program in SimBiology? And what confuse me is that the errors are only mention 'rule_12' to 'rule_18',but in fact I have similar rules form 'rule_11' to 'rule_19', why does this happens?

[Wei Wang]

To complement the warnings and errors I received below, there is some orhter statment when I use verify command:

which is: "

Test of RHS MATLAB code failed

Cause:

when call function 'calculateSolubilityCD', output parameter 'IntestinalAGSolubility' (may include other parameters) is not be assigned value. "

However, if I only use this function in my simulation program, simulation can run normally even if there is warnings in Equation view in SimBiology

If I combine another function in my model (because I try to bypass the tolerance problem):

function DrugBound=calculateDrugBound(Drug_1,CD_1,Kc,DrugMW,CDMW)
%Drug,CD in unit of mg/mL, Kc in unit of L/mmol
if CD_1==0
    DrugBound=0;
elseif CD_1>0
    Drug=Drug_1*1000/DrugMW;
    CD=CD_1*1000/CDMW;
    DrugBound=((Kc*(CD+Drug)+1)-((Kc*(CD-Drug)+1)^2+4*Kc*Drug)^(1/2))/(2*Kc);
end

I will get this message when run simulation:

The problem that parameter not be assigned value happens again.

This problem only takes place after I updated MATLAB to 2020a. I think these functions can work in 2019b. Is there anyway to fix it?

[Wei Wang]

Sorry for my long bothering and sorry for my question on MATLAB version.

I input the same function to project in 2019b and result in the same message. After I delete this part of function,

if CD_1==0
    DrugBound=0;
elseif CD_1>0

,I fix this problem.

But in terms of using Debugger tool, there are always warnings.

[Florian Augustin]

Hi,

It looks like the custom functions contain if statements that are not exhaustive. For example, calculateSolubilityCD(0, 0, 0) does not return a value because (0.4653*CDConcentration*1000/CDMW+0.4409)*AGMW evaluates to NaN. Could you replace the elseif statement with just else? Similarly in other custom functions. If CD_1 is negative (mind that small negative values, e.g. -1e-16, are possible due to numerical errors), then the return value of calculateDrugBound is undefined.

Once the simulations make it past t=0 without errors, the next step would be to find out why those function are called with unexpected values.

-Florian

[Wei Wang]

Hi,

Oh, yes. When I change elseif to if, I can run the simulation and result is good.

Now I am trying to use SimBiology Debugger, but I am confused with it:

It says my drug solubility in duodenum is NAN, due to the CDMolecularWeight is 0 as you warned me. But I don't know why it is possibly equal to 0. Its value is set to be 1 in model and another certain huge number in Variants. And I have cleared and reload my project before I start debug my model. Is there anything else I missed?

[Florian Augustin]

Hi,

To find out how the intial value is set within the model itself (e.g. by rules) you can inspect the component CDMolecularWeight in the debugger by selecting it as a custom component in the top bar of the app:

Then click visualize and click on "more details". This should give you a list of rules/reactions that determine the value of CDMolecularWeight. Note, that in the current version of the debugger, doses and variants are not visualized in this table/graph (although the reported values of components take doses and variants into account). Maybe this gives a hint why the value of CDMolecularWeight is 0.

You can uncheck response and expression plots for more real estate for the model component dependency graph.

-Florian

[Wei Wang]

Hi, Florian,

I find where the problem is. Yesterday I did not save the model after I changed the value of CDMolecularWeight from 0 to 1, so I received the result in Debugger despite I had cleared and reloaded the model. Today I retry it and find the problem disappears, then I realize I made it wrong.

Thanks for your help.

Wang Wei