SimBiology - sundials - ODE solver - CVODE

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

Answer 1

[Arthur Goldsipe]

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);

[mary]

Thanks Arthur for your answer. In my case, I have 500 reactions (which is actually the reduced model. The detailed one has more than thousands of reactions) which make it very difficult to write them down one by one. However, the function source_terms returns the variations of species concentrations by multiplying matrixes.

Could I conclude form your reply that there is no mean to use such a source term function within SimBiology ?

Best regards,

Mary

Answer 2

[Arthur Goldsipe]

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.

[mary]

Hi Arthur,

Unfortunately, for confiential reasons, I am not allowed to share these files. However, the structure of the latter is not very complicated. I have a matrix of reaction rate coefficients which depends on the temperature. A matrix filled by 0, 1 and -1 as stoichiometric coefficients of reactifs and products and species' concentration. The multiplication of these three gives the variation of species' concentration.

Do you have any example where SimBiology calls a source term function to calculate the appropriate terms ?

Regards,

Mary

[Arthur Goldsipe]

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.