Creating a 3-D bar graph with parameter sensitivities with respect to all model species

I would like to create a 3D bar graph that displays parameter sensitivities for different species in a simbiology model. I tried exporting the simbiology task code to MATLAB but I cannot run the script even after specifying the correct model file. I would like to do this with MATLAB code to create high quality publication figures for the sensitivity analysis. I cannot see any data generated using my code below:

sbioloadproject("GI_model.sbproj");
%v=getvariant(m2);
%set(v(1), 'Active'); %set variant that contains top 5 param.
cs = getconfigset(m4, 'active');
set(cs, 'StopTime', 56);
cs.TimeUnits= ('day');
%mathworks command line
%[t,r,outputFactors,inputFactors] = getsensmatrix(simdata,outputFactorNames,inputFactorNames);
outputNames= ["A", "B", "C", "D", "E", "F", "G", "H", "I"];
inoutNames= ["k1", "k2", "k3", "k4", "k5", "k6", "k7"];
varObj= m4.getvariant;
doseObj=m4.getdose;
simdata= sbiosimulate (m4,[],varObj(3), doseObj(2));
[t,r,outputFactors,inputFactors] = getsensmatrix(simdata,outputNames,inputNames);

Also attached is a figure that I would like to get using my code.

Thanks in advance!

Answer 1

[Arthur Goldsipe]

I'm guessing the configset is not configured for sensitivity analysis. You probably need to do something like the following before calling sbiosimulate:

cs.SolverOptions.SensitivityAnalysis = true;
sensOpts = cs.SensitivityAnalysisOptions;
sensOpts.Outputs = sbioselect(m1, "Name", outputNames);
sensOpts.Inputs = sbioselect(m1, "Name", inputNames);

[Fearass Zieneddin]

It worked but now how can I plot using bar3 commnad to show the result in a similar way to the figure above? I am using this function:

[t,r,outputFactors,inputFactors] = getsensmatrix(simdata,outputNames,inputNames);

[Arthur Goldsipe]

That depends... The sensitivities you've extracted are a function of time, but it looks like you want to replace each of these vectors with some sort of scalar "summary." So how do you want to summarize them? One common approach is to replace them with something like AUC (area under the curve). You could do that with something like the following:

for i = size(r,2):-1:1
    for j = size(r,3):-1:1
        auc(i,j) = trapz(t,r(:,i,j));
    end
end

[Fearass Zieneddin]

Yes, this is what I was trying to do I now have a figure similar to the one presented above. I still see my code creates a different result than the simbiology sensitivity task in the Model Analyzer App. Any idea what I might be doing wrong? You can see above I've mentioned the correct variant and dose objects before running the simulation. I am running the code for 3 different sesntivity options: 'full' 'half' and 'none'

But my resuts still look different from the model analyzer task.

Thanks again!

[Arthur Goldsipe]

My suggestion is to go back to the exported task code to see what SimBiology is doing. Ideally, I'd like to figure out why you couldn't run that code, since it would likely save you time in the future if you can just use exported code as the starting point for publication-quality plots. But for now you can just inspect that code to see what's different from the code you've written yourself, and copy and paste the bits you want to try to reuse.

[Fearass Zieneddin]

Here's a copy of my full code. The only difference is that I am converting the senstivity matrix output to AUC so I can get a 3D plot. Does the simbiology task compute the sensitivity of each output with respect to each input at all the time points that it steps through?

sbioloadproject("model.sbproj");
cs = getconfigset(m4, 'active');
set(cs, 'StopTime', 56);
cs.TimeUnits= ('day');
cs.SolverOptions.SensitivityAnalysis = true;

%mathworks command line
%[t,r,outputFactors,inputFactors] = getsensmatrix(simdata,outputFactorNames,inputFactorNames);
outputNames= ["A","B1", "B2","B3","B4", "B5", "C1", "C2", "C3", "C4"];
inputNames= ["fu_ratio", "k1", "k2", "k3", "k4", "k5", "k6"];

varObj= m4.getvariant;
doseObj=m4.getdose;
sensOpts = cs.SensitivityAnalysisOptions;

sensOpts.Inputs = sbioselect(m4, "Name", inputNames);
sensOpts.Outputs = sbioselect(m4, "Name", outputNames);
sensOpts.Normalization="None";
simdata= sbiosimulate (m4,[],varObj(3), doseObj(2));
[t,r,outputFactors,inputFactors] = getsensmatrix(simdata,outputNames,inputNames);
for i = size(r,2):-1:1
    for j = size(r,3):-1:1
        auc(i,j) = abs(trapz(t,r(:,i,j)));
    end
end
figure(1); bar3(auc);
xlabel('parameters'); ylabel('State variables');
set(gca, 'YTick', 1:length(outputNames), 'YTickLabel', outputNames);
set(gca, 'XTick', 1:length(inputNames),'XTickLabel',inputNames);

[Arthur Goldsipe]

Yes, the simbiology task computes the sensitivity of each output with respect to each input at all the time points that it steps through.

[Arthur Goldsipe]

Would you be willing to send me a message over MATLAB Answers with your contact information? I've asked some of my colleagues for help answering your questions. My main job is as a developer on the SimBiology team, and I think one of our customer-facing engineers would be better able to assist you.

[Fearass Zieneddin]

Another question: lets say we have species x1,x2, x3.. etc and to better represent all those we have a rule such: Σx= x1+x2+x3. Is there a way I can specify this rule as sensitivity output in my code? I have many different species that can reprsent a group of cells or protiens in my model and we are intrested in the outcome of varying model parameters on this group of cells or protiens.

I really apperciate your feedback!

[Arthur Goldsipe]

No, it's not currently possible to (directly) do sensitivity analysis on someting determined by a repeated assignment rule like this. The best workaround I've come up with so far is to do sensitivity analysis on x1, x2, and x3. Then, compute the sensitivity of Σx after the simulation as the sum of the relevant sensitivities. Sorry I can't offer anything better.

[Fearass Zieneddin]

I think I found the issue in the code: it looks like when using the 'full' or 'half' normalization options the auc matrix returns NaN values. The r(:,i,j) matrix shows some numbers any idea what I am doing wrong here?