MATLAB: Creating a mesh figure

mesh figureSimBiology

hello all,
I have written a code where i give in calcium input to simbiology model and see the efects of how different number of spikes and frequency effect the output. spike train is my own function where i give in calcium as a rule. i have already made a graph like spike_freq.jpg attached. where i compare different spike numbers and frequencies in output species. this is the one that follows the code.
now I want an additional figure as shown in mesh.jpg where i can compare the amplitude and integral at each spike with all frequencies.
the code i tried for this is the one commented at end.
any help will be highly appreciated.
thanks
parul
Below is the code I used
function  Bito_like_prediction
proj=sbioloadproject('CaMK');
modelobj=proj.m1;
%startTime is time for relaxation
startTime=30;
stopTime=40;
times=0:0.01:stopTime;
frequencies=[2 3.3 5 6.7 10 20];
nFreq=length(frequencies);
maxSpikeNum= [10 20 30];
nSpi=length(maxSpikeNum);
max_step=0.005;
 y=zeros(1,nFreq,nSpi);
%y1=zeros(times,nfreq,nspikes);
output_species=[6 12 33]; 
configsetObj = getconfigset(modelobj);
set(configsetObj, 'StopTime', times(end));
set(configsetObj, 'SolverType','ode15s');
set(configsetObj.SolverOptions, 'OutputTimes', times);
set(configsetObj.SolverOptions, 'MaxStep', max_step);
set(configsetObj.SolverOptions, 'AbsoluteTolerance', 1.0e-8);
set(configsetObj.SolverOptions, 'RelativeTolerance', 1.0e-6);
set(configsetObj.RuntimeOptions, 'StatesToLog', 'all');
%  %set all initial amounts to 0
 for ii=1:size(modelobj.species,1),
      modelobj.species(ii).InitialAmount=0;
 end
%  
% %activate all rules, except the last one(?)
for ii=1:size(modelobj.Rules,1),
     modelobj.Rules(ii).Active=1;
end
modelobj.species(1).InitialAmount=10000; 
modelobj.species(6).InitialAmount=50; 
modelobj.species(7).InitialAmount=4000; 
modelobj.species(17).InitialAmount=20000; 
modelobj.species(32).InitialAmount=5000; 
modelobj.Rules(49).Active=1;
modelobj.Parameters(104).Value=startTime;
plot_counter=1;
  for  j=1:nSpi
      for i=1:nFreq
    %if frequencies(i)>0
      modelobj.Rules(49).Active=1;
      modelobj.Parameters(100).Value=frequencies(i);
      modelobj.Parameters(107).Value=maxSpikeNum(j);
      [t,y,names]=sbiosimulate(modelobj);
        
    
   figure (1) 
subplot(nSpi,nFreq,plot_counter)
plot(t(3000:end),y(3000:end,output_species))
legend(names{output_species(1)},names{output_species(2)}, names{output_species(3)})
title(['Frequency: ' num2str(frequencies(i)) 'Nr spikes: ' num2str(maxSpikeNum(j))]);
axis([30 40 0 15000])
plot_counter=plot_counter+1;
% [X,Y] = meshgrid(nSpi,nFreq);
% %amplitude is Z
%  z(1,:,i) = trapz(t,y(:,:,i));
% C = del2(Z(1,:,i));
% 
% figure(2)
% mesh(X,Y,Z,C,'FaceLighting','gouraud','LineWidth',0.3)
% %surf(X,Y,F)
end
  end
end

Best Answer

Hi Parul,
If I understand correctly what you are trying to do, then you are very much on the right track. So it is probably only a matter of reorganizing the code for when data is computed (to avoid overwriting variables) and when it is plotted.
Here is a snippet of code that may help. I am using a SimFunction to run the model simulations, but you can equally use sbiosimulate as you do in your code.
% Load SimBiology model:
mStruct = sbioloadproject('lotka'); 
lotkaModel = mStruct.m1;
 
% Configure solver:
configSet = getconfigset(lotkaModel);
configSet.SolverType = 'sundials';
configSet.SolverOptions.AbsoluteTolerance = 1e-8;
configSet.SolverOptions.RelativeTolerance = 1e-6;
 
% Define values for species x and z to scan over:
xValues = linspace(0.8,   1, 20); 
zValues = linspace(  0, 0.2, 3); 
 
% Create cross product of all combinations of values for x an z:
[X, Z] = meshgrid(xValues, zValues);
xzValueCombinations = [X(:), Z(:)];
 
% Create SimFunction for simulation:
paramNames  = { 'x', 'z'}; 
observables = {'y1', 'y2'};
dosingInfo  = [];
simFun = createSimFunction(lotkaModel, paramNames, observables, dosingInfo);
 
% Run simulations:
stopTime = 10;
simData  = simFun(xzValueCombinations, stopTime);
 
% Prepare variable for storing integral values:
integralY1 = nan(size(X));
integralY2 = nan(size(X));
 
% Compute integral of time courses for y1 and y2:
numSimulations = numel(simData);
for i = 1:numSimulations
    time = simData(i).Time;
    data = simData(i).Data;
    integralY1(i) = trapz(time, data(:,1));
    integralY2(i) = trapz(time, data(:,2));
end
 
% Plot integral values
figure(1); clf; 
 
subplot(1,2,1);
mesh(X, Z, integralY1);
xlabel(paramNames{1});
ylabel(paramNames{2});
zlabel(['integral of ', observables{1}]);
 
subplot(1,2,2);
mesh(X, Z, integralY2);
xlabel(paramNames{1});
ylabel(paramNames{2});
zlabel(['integral of ', observables{2}]);
I hope this helps. Let me know if this does not answer your question.
Best,
-Florian
Related Question