MATLAB: How to optimize parallel simulations of a SimBiology model

MATLABParallel Computing ToolboxparforsbioaccelerateSimBiologysimfunction

In March of 2020, there was a question about parfor loops and simbiology. (https://www.mathworks.com/matlabcentral/answers/510565-issues-with-doses-variants-in-parfor-loop-in-simbiology?s_tid=srchtitle)
I currently use a parfor loop, but my model is fairly large and uses a lot of molecules — around 40 species and 40 reactions. As a result, it takes a very long time to obtain my results. I was wondering if you think SimFunction (along with 'UseParallel' and 'AutoAccelerate') would be appropriate for me to use. — Below, you'll find that I use a parfor appoach and make sure the files are accessible through the latter two functions. After seeing the aforementioned post, I'm wondering if I should be using SimFunction, sbiosimulate, sbioaccelerate, sbioensemblerun, or another method to optimize my approach.
Here is the general structure for what I am trying to accomplish:
I am using SimBiology 2019b. The solver I need to use is SSA, monte carlo style. I have a chemical reaction network set up. I add an event at 3000 seconds to set the concentration of an 'input' species to a given value.
I have a species in the chemical reaction network that I view as the 'output' species. My goal is to run 100 simulations for each given input amount. Then, I would like to increase the input amount over 100 steps in a uniformly increasing manner at each step and observe the behavior of the 'output' species.
Here is the structure of my code:
for i = 1:100 %100 steps with uniform increasing 'input' species amounts
    sbioloadproject('model.sbproj');
    cs = getconfigset(m1, 'default');
    cs.SolverType = 'ssa';
    solver = cs.SolverOptions;
    solver.LogDecimation = 1000;
    inputSpecies = sbioselect(m1,'Type','species','Name','inputSpecies'); %The 'input' species
    inputSpecies.InitialAmount = 0;
    parameterTimeunits = addparameter(m1,'Timeunits',1.0,'ValueUnits','minute','ConstantValue',true);
    parameterMolunits = addparameter(m1,'Molunits',1.0,'ValueUnits','molecule','ConstantValue',true);
    
    inputSpecies = addevent(m1, '(time/Timeunits) >= 3000', ['inputSpecies = (1500000*' num2str(i) ')*Molunits']);
    sbiosaveproject modelParallel.sbproj m1 inputSpecies
    clear m1;
    
    %simulate the model
    parfor j=1:100 %Number of simulation runs per input step
        [m1,inputSpecies] = loadSimBiologyModel('modelParallel.sbproj');
        [t,x,names] = sbiosimulate(m1);
        parsave(['StochasticTest_Input' num2str(i) '_Instance' num2str(j)], inputSpecies, t, x, names);
    end
end
function [m2,inputSpecies2] = loadSimBiologyModel(filename)
    sbioloadproject(filename);
    m2 = m1;
    inputSpecies2 = inputSpecies;
end
function parsave(fname,inputSpecies,t,x, names)
    save(fname,'inputSpecies','t','x','names');
end
%I'll then take all these datafiles, interpolate them, and combine it into one dateset.
%Then I'll observe the amount of output over time for the varying range of input.
With my current approach, I haven't become very familiar with parameters or doses. Would you be able to provide me some insights on A.) if SimFunciton (along with 'UseParallel' and 'AutoAccelerate') is appropriate and if so, B.) since the parameters and doses don't seem applicable to my setup, how can I adapt the function inputs appropriately? Any tips, suggestions, or critiques are happily welcomed. If you need any more info, please let me know.

Best Answer

Hi Colton,
Yes, using SimFunction would be better as it offers easy and efficient parallelization and avoids repeated compilation of the model. Also, instead of saving 100 variations of the model, you can update the event to set inputSpecies to some new parameter value that is an input to the SimFunction. For example, your event could look like:
addevent(m1,'(time/Timeunits) >= 3000', 'inputSpecies = newParam');
you could create simfunction as:
F = createSimFunction(m1,{'newParam'}, {'output'}, 'UseParallel', true);
and you could then run the simfunction as:
simdata = F(phi, t_stop)
where phi is an array of input amount values.
If you could avoid using events as you mentioned and just change the initial conditions to steady state values instead, this should make the simulations a lot faster. Once you eliminate events, you can also look into using exlptau or impltau instead of ssa to speed up simulations. These are faster than SSA but introduce approximations which you can read about here: https://www.mathworks.com/help/simbio/ug/stochastic-solvers.html
Additional things to consider:
  • Since you’re only interested in one output, it would be most efficient to only log this one output species when creating the SimFunction. This is the observables input to createSimFunction.
  • Lastly, SSA is probably reporting more simulation output times than you may be interested in. This could also slow things down due to the extra memory usage. So you may want to investigate reducing the number of times using OutputTimes (can be an input to SimFunction) or LogDecimation.
Fulden
Related Question