Scanning Simbiology variants in Matlab for virtual population simulations

1 view (last 30 days)
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

Accepted Answer

Arthur Goldsipe
Arthur Goldsipe on 3 Sep 2021
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
  5 Comments
Arthur Goldsipe
Arthur Goldsipe on 7 Sep 2021
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
Fearass Zieneddin on 9 Sep 2021
Edited: Fearass Zieneddin on 9 Sep 2021
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;

Sign in to comment.

More Answers (0)

Categories

Find more on Perform Sensitivity Analysis in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!