Simulating multiple model parameterizations using scenarios and simfunctions

Hello,

I'm having trouble simultaneously simulating multiple parameterizations of a model using simfunctions and scenarios, in R2019b. Please refer to the example below. Is there a proper way of doing this? I think I could get around this error be inputting all of the parameterizations when calling 'createSimFunction', but this is not possible in my application because I do not know a priori how many different parameterizations I plan to simulate.

Thank you,

Abed

% The following works:
sbioloadproject('insulindemo','m1');
variants = getvariant(m1);
singleMeal = sbioselect(m1,'Name','Single Meal');
sObj = SimBiology.Scenarios;
add(sObj,'cartesian','variants',variants);
add(sObj,'cartesian','dose',singleMeal);
sObj.verify(m1);
f = createSimFunction(m1,sObj,{'[Plasma Glu Conc]','[Plasma Ins Conc]'},[]);
sd = f(sObj,24);

% But the following triggers an error:
nIndivs = 2; % number of parameterizations (let's pretend each individual has different parameterizations -- although in this example they have the same parameterization).
for iEntry = 1:sObj.NumberOfEntries
    sObj = sObj.updateEntry(iEntry,'Content',repmat(sObj.getEntry(iEntry).Content,nIndivs,1));
end
sObj.verify(m1);
sd2 = f(sObj,24);

Answer 1

[Abed Alnaif]

Hello,

Fulden emailed me an example of how to get this to work, and I'm posting it here for reference. Thank you, Florian and Fuilden, for your help with this. This is my first time using the new scenarios object, and I'm really liking it now that I experienced its flexibility and how one can construct complex combinations with very few lines of code. It's a very nice addition to SimBiology! It would be good in the future if the 'sbiofit' command allowed the user to specify different scenarios for each of the groups.

clc
clear

m = sbmlimport('lotka');
cs = getconfigset(m);
set(cs.SolverOptions, 'RelativeTolerance', 1e-6);

% Define doses
d1 = sbiodose('dose for y1');
set(d1, 'TargetName', 'y1');
set(d1, 'Amount', 50);

d2 = sbiodose('dose for y2');
set(d2, 'TargetName', 'y2');
set(d2, 'Amount', 50);

% Define variants
v1 = sbiovariant('individual1');
v1.addcontent({'species', 'y1', 'value', 900});

v2 = sbiovariant('individual2');
v2.addcontent({'species', 'y1', 'value', 900});

v3 = sbiovariant('individual3');
v3.addcontent({'species', 'y1', 'value', 900});


% Create scenarios over variants and doses:
sObj = SimBiology.Scenarios('individuals', [v1, v2]);
add(sObj, 'cartesian', 'doses', [d1, d2]);

% Print content of scenarios object
disp(sObj)

% Create SimFunction and evaluate
f = createSimFunction(m, sObj, {'y1', 'y2'}, [], 'AutoAccelerate', false);

sd = f(sObj, 10);
sbioplot(sd);

% Add additional values to scenarios
sObj.updateEntry('individuals','Content', [v1, v2, v3]);

% Print content of scenarios object
disp(sObj)

sd = f(sObj, 10);
sbioplot(sd);

Answer 2

[Florian Augustin]

Hi Abed,

The short answer why you see the error message is that the SimFunction is being created with one dose input, and it is called (after the modification of the scnearios object) with two dose inputs. In the scearios object, each element in the dose entry is treated as a separate dose input to the SimFunction. The immediate fix would be to create a new SimFunction after you updated the scenarios object.

But this may be expensive in your application and possibly unnecessary. In your example you mention different parameterizations of your doses. In this case, you may be able to use parameterized doses? You could create a scenarios object containing the parameterization parameter(s) instead of the dose objects:

sbioloadproject('insulindemo', 'm1');
variants = getvariant(m1);
singleMeal = sbioselect(m1, 'Name', 'Single Meal');
addparameter(m1, 'DoseAmount', 78, 'Units', 'gram');
set(singleMeal, 'Amount', 'DoseAmount');

sObj = SimBiology.Scenarios;
add(sObj,'cartesian', 'variants', variants);
add(sObj,'cartesian', 'DoseAmount', [75, 78, 85]);

sObj.verify(m1);
f = createSimFunction(m1, sObj, {'[Plasma Glu Conc]', '[Plasma Ins Conc]'}, singleMeal);
sd = f(sObj, 24, getTable(singleMeal));

There may be other options, but those probably depend on your application.

I hope this helps, but feel free to follow up if there's anything we can help with.

Best,

-Florian

[Abed Alnaif]

Hi Florian,

I'm actually planning to explore non-dose-related parameters. And I don't know a priori how many sets of parameters I planning to simulate. And I'd like to avoid having to recreate and reaccelerate the simfunction every time.

Is this possible? I think this would be possible if calling the simfunction with variants/doses instead. I'm not sure why it doesn't work when calling the simfunction with scenarios?

Thanks,

Abed

[Florian Augustin]

Hi Abed,

I think the answer depends on your use case. Would you be able to share an example? You should be able to do everything you can do with plain SimFunctions with scenarios as well.

Below are some general thougts that may help.

If you want to use flexible scenarios to explore different parameterizations or dosing regimen of your model (rather than different values), then you could create scenarios objects with 'fixed/template' values as placeholders for future variations. Here is an example.

m = sbmlimport('lotka');
cs = getconfigset(m);
set(cs.SolverOptions, 'RelativeTolerance', 1e-6);

d1 = sbiodose('dose for y1');
set(d1, 'TargetName', 'y1');
set(d1, 'Amount', 50);

% 'Void' dose (amount set to 9), but it will provide a hook to explore 
% dose d2 after the SimFunction was created.
d2 = sbiodose('another dose for y1');
set(d2, 'TargetName', 'y1');
set(d2, 'Amount', 0); 

sObj = SimBiology.Scenarios;
add(sObj, 'cartesian', 'y1', [850, 900, 950]);
add(sObj, 'cartesian', 'dose for y1', d1);
% Add a placeholder entry so values of y2 can be varied later. The same 
% idea works for variants. Just extend the variant content to contain
% nominal values of species/parameters you may want to vary in the future.
add(sObj, 'cartesian', 'y2', 900);
% Add a placeholder for potential future dosing of y1; here the doses are 
% put in different entries, but mixed elementwise/Cartesian combinations 
% of doses can be created as well if doses need to be applied jointly.
add(sObj, 'cartesian', 'another dose for y1', d2); 

f = createSimFunction(m, sObj, {'y1', 'y2'}, [], 'AutoAccelerate', false);
sd = f(sObj, 10);

sbioplot(sd)

% Now explore 'another dose for y1':
templateDoseEntryOnScenario = getEntry(sObj, 'another dose for y1');
templateDoseEntryOnScenario.Content.Amount = 100;

% Or explore different values of y2
sObj = updateEntry(sObj, 'y2', 'Content', [880, 920]);

sd = f(sObj, 10);

sbioplot(sd)

Unfortunately, this requires you to create a framework of dosing regimen in advance.

If the bottleneck is the acceleration of the SimFunction, then the new name-value pair 'AutoAccelerate' in createsimfunction (added in R2019b) may save some time while building the scenarios object.

-Florian

[Abed Alnaif]

Hi Florian,

That's interesting, it seems like your example works for the 'cartesian' case, but I can't get it to work when combining the entries 'elementwise'. Please refer to the example below, which is just a modification of your example to have the entries combine 'elementwise'. Can you get this one to work?

Thanks,

Abed

m = sbmlimport('lotka');
cs = getconfigset(m);
set(cs.SolverOptions, 'RelativeTolerance', 1e-6);
d1 = sbiodose('dose for y1');
set(d1, 'TargetName', 'y1');
set(d1, 'Amount', 50);
% 'Void' dose (amount set to 9), but it will provide a hook to explore 
% dose d2 after the SimFunction was created.
d2 = sbiodose('another dose for y1');
set(d2, 'TargetName', 'y1');
set(d2, 'Amount', 0); 
sObj = SimBiology.Scenarios;
add(sObj, 'elementwise', 'y1', [850, 900, 950]);
add(sObj, 'elementwise', 'dose for y1', [d1, d1, d1]);
% Add a placeholder entry so values of y2 can be varied later. The same 
% idea works for variants. Just extend the variant content to contain
% nominal values of species/parameters you may want to vary in the future.
add(sObj, 'elementwise', 'y2', [900, 901, 902]);
% Add a placeholder for potential future dosing of y1; here the doses are 
% put in different entries, but mixed elementwise/Cartesian combinations 
% of doses can be created as well if doses need to be applied jointly.
add(sObj, 'elementwise', 'another dose for y1', [d2, d2, d2]); 
f = createSimFunction(m, sObj, {'y1', 'y2'},[]);
sd = f(sObj, 10);
sbioplot(sd)
% Now we would like to explore one additional individual:
sObj.updateEntry('y1','Content',[850, 900, 950, 1000]);
sObj.updateEntry('y2','Content',[900, 901, 902, 903]);
sObj.updateEntry('dose for y1','Content',[d1, d1, d1, d1]);
sObj.updateEntry('another dose for y1','Content',[d2, d2, d2, d2]);
sd = f(sObj, 10);
sbioplot(sd)

[Florian Augustin]

Hi Abed,

You are right, unfortunately there are cases where you will have to create a new SimFunction. Your example points out one of those: if you combine different doses elementwise with value entries, then you will need to create a new SimFunction if the number of doses in the scenarios changes. Otherwise, you'd need to do some unnecessary simulations which may be more wasteful than creating a new SimFunction.

Since you mentioned that you are actually planing on exploring non-dose related parameters, here is another thought that may help. If the dosing stays the same accross all scenarios, you could do the following:

m = sbmlimport('lotka');
cs = getconfigset(m);
set(cs.SolverOptions, 'RelativeTolerance', 1e-6);

% Define doses
d1 = sbiodose('dose for y1');
set(d1, 'TargetName', 'y1');
set(d1, 'Amount', 50);

d2 = sbiodose('dose for y2');
set(d2, 'TargetName', 'y2');
set(d2, 'Amount', 50);

% Create scenarios:
sObj = SimBiology.Scenarios;
add(sObj, 'elementwise', 'y1', [850, 900, 950]);
add(sObj, 'elementwise', 'y2', [900, 901, 902]); 

% Create SimFunction and evaluate
f = createSimFunction(m, sObj, {'y1', 'y2'}, [d1, d2]);

sd = f(sObj, 10, getTable([d1, d2]));

sbioplot(sd)

% Add additional values to scenarios
sObj.updateEntry('y1','Content',[850, 900, 950, 1000]);
sObj.updateEntry('y2','Content',[900, 901, 902, 903]);

sd = f(sObj, 10, getTable([d1, d2]));

sbioplot(sd)

You will get the same behavior if you bundled the doses in a separate scenarios object:

m = sbmlimport('lotka');
cs = getconfigset(m);
set(cs.SolverOptions, 'RelativeTolerance', 1e-6);

% Define doses
d1 = sbiodose('dose for y1');
set(d1, 'TargetName', 'y1');
set(d1, 'Amount', 50);

d2 = sbiodose('dose for y2');
set(d2, 'TargetName', 'y2');
set(d2, 'Amount', 50);

% Create scenarios:
sObj = SimBiology.Scenarios('y1', [850, 900, 950]);
add(sObj, 'elementwise', 'y2', [900, 901, 902]); 

sObjForDoseBundling = SimBiology.Scenarios('dose1', d1);
sObjForDoseBundling.add('elementwise', 'dose2', d2);

sObj = sObj.add('cartesian', sObjForDoseBundling);

% Create SimFunction and evaluate
f = createSimFunction(m, sObj, {'y1', 'y2'}, []);

sd = f(sObj, 10);

sbioplot(sd)

% Add additional values to scenarios
sObj.updateEntry('y1', 'Content', [850, 900, 950, 1000]);
sObj.updateEntry('y2', 'Content', [900, 901, 902, 903]);

sd = f(sObj, 10);

sbioplot(sd)

Best,

-Florian

[Abed Alnaif]

Hi Florian,

In my application, the doses aren't the same across the scenarios. In some scenarios, a dose is not active, in which case I replaced the dose object with an empty dose object. So, I can't use 'cartesian' combinations.

Thanks,

Abed