23 views (last 30 days)

Hi,

Could someone explain me how I can use the sundials solver (CVODE), through my use of simBiology toolbox ?

I know ODE solvers of Matlab and I used to use ODE15s. However, the latter fails in integration of stiff problems. So, I tried to use CVODE solver. The same issue occurs (if you are interested to know more about that please see my post on sundials forum : http://sundials.2283335.n4.nabble.com/force-the-CVODE-solver-to-continue-despite-the-error-test-failed-repeatedly-or-with-h-hmin-td4655576.html ). So, I ended up with the following troubleshooting solution being suggested for the very same solver included in SimBiology package ( https://nl.mathworks.com/help/simbio/ug/troubleshooting-simulation-errors.html ). Here you can see that they talk about two options of MaximumNumberOfLogs and MaximumWallClock, the use of which can help to avoid this problem. So, as I can not find these two options in CVODE directly downloaded from https://computation.llnl.gov/projects/sundials/sundials-software , I thought maybe these two options are only available in the SimBiology. Now, the problem is the following : there is no enough information/tutorial on this tool. How I can chage the following line to make use of SimBiology ?

[T,Y]=ode15s(@source_terms,[0 tend],yinitial,options);

I appreciate a lot your help,

Best regards,

Mary

Arthur Goldsipe
on 15 Apr 2019

Hi Mary,

You are correct that MaximumNumberOfLogs and MaximumWallClock are SimBiology-specific options.

In order to simulate your model in SimBiology, you will have to translate the differential equations inside your function source_terms into a SimBiology model. We have designed SimBiology primarily for simulating systems of reactions, so the most common way to build a SimBiology model is by constructing a series of reactions. You can find a simple example of that here. If you prefer to translate a series of differential equations into a SimBiology model, you may prefer to define your model in terms of rate rules. You can find an example of that here.

Once you have defined your model (let's say, in variable modelObj), here's sample code that sets the solver to sundials, set MaximumNumberOfLogs and MaximumWallClock, and simulate the model:

configsetObj = getconfigset(modelObj);

configsetObj.SolverType = 'sundials';

configsetObj.MaximumNumberOfLogs = 100;

configsetObj.MaximumWallClock = 30;

[time, states] = sbiosimulate(modelObj);

Arthur Goldsipe
on 15 Apr 2019

Hi Mary,

It is technically possible to use arbitrary MATLAB functions with SimBiology. So you could in theory create a SimBiology model that calls source_terms to calculate the appropriate terms. But I would not recommend this approach. I think it would end up being very hard to read and extremely inefficient (that is, slow to simulate).

So, for all practical purposes, I would say you should not directly use source_terms with SimBiology. Instead, I would focus on finding an efficient way to convert source_terms into a SimBiology model.

Can you share your full code for source_terms? Depending on how that's written, I might be able to offer suggestions for writing a program to automatically create a SimBiology model from the file.

Arthur Goldsipe
on 15 Apr 2019

Let me offer an alternative suggestion. If you already have reaction rates and stoichiometries stored in matrices, then it should be relatively straightforward to use these matrices to create a SimBiology model programmatically. Here's an example that recreates the Lotka-Volterra predator-prey model:

%% Define variables

speciesNames = {'y1', 'y2', 'z'};

speciesInitialAmounts = [900 900 0];

reactantStoichiometry = [1 1 0; 0 1 1; 0 0 0];

productStoichiometry = [2 0 0; 0 2 0; 0 0 1];

parameterNames = {'c1', 'c2', 'c3'};

parameterValues = [10 0.01 10];

%% Create model

% Create blank model

modelObj = sbiomodel('Lotka-Volterra');

% Add species

for i = 1:numel(speciesNames)

speciesObj(i) = addspecies(modelObj, speciesNames{i}, speciesInitialAmounts(i));

end

% Add parameters

for i = 1:numel(parameterNames)

addparameter(modelObj, parameterNames{i}, parameterValues(i));

end

% Create reactions

for i = 1:numel(parameterNames)

reactionObj = addreaction(modelObj, 'null -> null');

kineticLawObj = addkineticlaw(reactionObj, 'MassAction');

kineticLawObj.ParameterVariableNames = parameterNames{i};

for j = 1:numel(speciesNames)

if reactantStoichiometry(j,i) > 0

addreactant(reactionObj, speciesObj(j), reactantStoichiometry(j,i));

end

if productStoichiometry(j,i) > 0

addproduct(reactionObj, speciesObj(j), productStoichiometry(j,i));

end

end

end

%% Configure simulation options and simulate the model

configsetObj = getconfigset(modelObj);

configsetObj.SolverType = 'sundials';

configsetObj.MaximumWallClock = 30;

configsetObj.MaximumNumberOfLogs = 100;

configsetObj.StopTime = 10;

[time, states] = sbiosimulate(modelObj);

plot(time,states)

Depending on the details of the temperature-sensitivity of your parameters, you might be able to create an initial assignment rule that initializes them property.

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply TodayMore Answers in the SimBiology Community

## 0 Comments

Sign in to comment.