How can I optimize parallel simulations of a SimBiology model?

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.

[Colton Harper]

Here is an additional note. If there are limitations with triggering an event in one of the possible functions to speed this process up, I think I can get around using an event altogether.

The reason I use events is to trigger the input amount when the system reaches steady state. However, I was just thinking I can adapt the model species' initial amounts to match the average steady state values, essentially starting the simulation at the average steady state. Then, instead of triggering an event at time 3000, I can just set the initial concentration of the 'input' species accordingly.

[Fulden Buyukozturk]

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, then 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

[Fulden Buyukozturk]

I should have also pointed out that if you could go around using events and just change the initial conditions to steady state values instead, this should make the simulations a lot faster.

One way to change the initial conditions to steadystate conditions would be to run sbiosteadystate and commit the steady state variant to the model. Note that this would require temporarily changing to a deterministic solver.

PS. I will also copy and paste my initial comment above to the answer field.

Answer 1

[Fulden Buyukozturk]

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

[Colton Harper]

Thank you, Fulden, for such a detailed, informative, and quick response!

I ran an ensemble of just 3 runs using my previous method and compared it to your suggested simfunction implementation and reducing the observable species. The parfor method took about 330 seconds. The simfunction method only took 115 seconds!

I'll try some of your other suggestions too -- starting from steady state species' amounts, changing LogDecimation and OutputTimes values, etc.

When trying to run an ensemble using SimFunction is it efficient to just specify the duplicate input amounts for the ensemble in phi?

For example 3 different input values, 2 simulation runs per value:

phi = [1.5e6; 1.5e6; 3e6; 3e6; 4.5e6; 4.5e6]

Thanks again!

[Arthur Goldsipe]

I can chime in to answer your followup question. Yes, it should be efficient to duplicate your parameter values to get ensemble results.

Glad to hear that you're seeing a significant speedup already. I would say that next biggest improvement would come from eliminating the event and switching and changing your initial conditions to steady state values. Switching the solver type to one of the tau leap solvers would be the next thing I'd suggest trying. This will also reduce the number of logged points, which would means it might not make that big of a difference to use LogDecimation. Also, I shouldn't have told Fulden to suggest using OutputTimes. I forgot that OutputTimes is only available for deterministic solvers.

-Arthur

[Colton Harper]

Excellent, thank you for following up, Arthur. This is very helpful!