Scanning Simbiology variants in Matlab for virtual population simulations

Hello,

I have a Simbioogy variant that I want to scan it's parameter values for 2-fold changes over 10 runs. First, how can I scan all parameters included in the variant for 2-fold changes? Can this be done using Simbiology variants or do I need to select each parameter from my model using: p1= sbioselect(m2, 'Name', 'kNeu2NA0', 'Type', 'parameter'); %m2 is my model

kNeu2NA0 is a parameter value I'm selecting.

I know this can be simulated using the Simbiology model analyzer app but I would like to use Matlab code for now. Also, is using a simbiology variant object a good way of creating 10 different patient parameters in a for loop?

Thanks in advance,

some code:

%simulating the virtual population for 5FU

sbioloadproject human_model_5FU;

v=getvariant(m2);

set(v(1), 'Active'); %set variant that contains top 5 param.

p1= sbioselect(m2, 'Name', 'kNeu2NA0', 'Type', 'parameter');

p2= sbioselect(m2, 'Name', 'kE2D0', 'Type', 'parameter');

p3= sbioselect(m2, 'Name', 'kDgCXC0', 'Type', 'parameter');

p4= sbioselect (m2, 'Name', 'kNeu_LP2L0', 'Type', 'parameter');

%two fold changes for first variable

for i= 1:10

a= [p1.value/2 p1.value*2];

b= [p2.value/2 p2.value*2];

c= [p3.value/2 p3.value*2];

d= [p4.value/2 p4.value*2];

r(i)=(a(2)-a(1).*(rand+a(1)));

t(i)=(b(2)-b(1).*(rand+b(1)));

u(i)=(c(2)-c(1).*(rand+c(1)));

p(i)=(d(2)-d(1).*(rand+d(1)));

vObj = addvariant(m2, 'vsim');

addcontent(vObj, {'parameter', 'kNeu2NA', 'Value', r(i)});

addcontent(vObj{'parameter', 'kE2D0', 'Value', t(i)});

addcontent(vObj{'parameter', 'kDgCXC0', 'Value', u(i)});

addcontent(vObj{'parameter', 'kNeu_LP2L0', 'Value', p(i)});

vObj.Active = true;

end

Answer 1

[Arthur Goldsipe]

There are lots of ways you could do this. And you certainly could use variants. But I personally would do this by creating a SimFunction from the model and using Scenarios to represent the parameter samples. This is similar to the code that is used in the Model Analyzer. In fact, I created this example by viewing and modifying code created by the Model Analyzer.

sbioloadproject lotka m1

paramNames = ["x", "y1", "y2"];
outputNames = ["y1", "y2"];
simfunc = createSimFunction(m1, paramNames, outputNames, [])

simfunc = 
SimFunction

Parameters:

     Name     Value       Type    
    ______    _____    ___________

    {'x' }       1     {'species'}
    {'y1'}     900     {'species'}
    {'y2'}     900     {'species'}

Observables: 

     Name        Type    
    ______    ___________

    {'y1'}    {'species'}
    {'y2'}    {'species'}

Dosed: None

paramValues = simfunc.Parameters.Value;

samples = SimBiology.Scenarios();
for i = 1:numel(paramNames)
    probdist = makedist("Uniform", "lower", 0.5*paramValues(i), "upper", 2*paramValues(i));
    add(samples, "elementwise", paramNames(i), probdist, "Number", 10);
end

simdata = simfunc(samples, 2);
sbioplot(simdata);

I hope this helps.

-Arthur

[Fearass Zieneddin]

Thank you so much that really helps! Is there a better way to see the responses seperately for example I'd like to see each species response on a seperate plot. Also, how do I specify the dose in the simfunction? I know the ''[ ]'' represents the dosed species in:

simfunc = createSimFunction(m1, paramNames, outputNames, []);

corect?

Thanks

[Arthur Goldsipe]

One way you could put each plot into a separate response would be to extract and plot each response from the SimData array. Something like this (continuing from the code in my answer):

for i = 1:numel(outputNames)
    sbioplot(selectbyname(simdata, outputNames(i)))
end

The approach I would suggest for dosing depends on exactly what differs (if anything) about the doses from one simulation to the next. If you want a separate dose for each, you could just include the doses in the Scenarios and use the Scenarios to create and execute the SimFunction. Again, continuing from my sample code, here's how you could create a separate dose for each simulation:

for i = 1:10
    doses(i) = sbiodose("d" + i);
    doses(i).TargetName = "z";
    doses(i).Amount = 10 + i;
end
samples2 = copy(samples);
add(samples2, "elementwise", "doses", doses);

simfunc = createSimFunction(m1, samples2, outputNames, []);
simdata = simfunc(samples2, 2);

Another way to represent the above example would be to use a parameterized dose. This involves having a parameter for the dose property of interest, say "doseAmount" for the dose amount. And then you can add this parameter to your Scenarios:

samples3 = copy(samples);
m1.addparameter('doseAmount');
paramDose = sbiodose("paramDose");
paramDose.TargetName = "z";
paramDose.Amount = "doseAmount";
add(samples3, "elementwise", "doseAmount", 11:20);

doseTable = getTable(paramDose);
simfunc = createSimFunction(m1, samples3, outputNames, paramDose);
simdata3 = simfunc(samples3, 2, doseTable);

[Fearass Zieneddin]

I actuly realized that I don't need to simulate 10 different doses. I just had to specify the dose programmaticly without the for loop. Something else I am trying is to plot the mean populaiton behavior. I am looking at the addobservable function but can't get it to work becuase I am using the wrong syntax. Any example using the same code as before?

sdout = addobservable(simdata,outputNames(i), mean(simdata))

Thanks!

[Arthur Goldsipe]

So you want to average the results of 10 simulations? You can't use an observable for that. An observable can only operate on a single simulation at a time. So you could use an observable to calculate statistics separately for each simulation, but not statistics of multiple simulations.

I would probably average simulation results with code like the following:

% Resample all simulations to a common time vector
stopTime = simdata(1).Time(end);
timeVector = linspace(0,stopTime,1001);
simdataResampled = resample(simdata, timeVector);

% "Stack" the matrices of simulation results and average them
stackedData = cat(3, simdataResampled.Data);
meanData = mean(stackedData, 3);

% Plot the results
plot(timeVector, meanData);
legend(simdata(1).DataNames);

[Fearass Zieneddin]

Thank you so much! I was able to get the mean simulation behavior to compare with min and max behavior with similar code like this:

stopTime = simdata(1).Time(end);
timeVector = linspace(0,stopTime,843);
simdataResampled = resample(simdata, timeVector);
% "Stack" the matrices of simulation results and average them
stackedData = cat(3, simdataResampled.Data);
meanData = mean(stackedData, 3);
maxData = max(stackedData, [],3);
minData = min(stackedData, [],3);
for i= 1:numel(outputNames)
    figure(i);
    plot(timeVector, meanData(:,i),'color', 'red');
    hold on;
    plot(timeVector, maxData(:, i), 'color', 'blue');
    hold on;
    plot(timeVector, minData(:, i), 'color', 'green');
    hold on;
    xlabel('time (day)');
    set(gcf,'color','w');
    set(gca,'fontsize',11)
    set(gca,'box','off');
end
hold off;