creating a mesh figure

7 views (last 30 days)
Parul Tewatia
Parul Tewatia on 7 Feb 2019
Commented: Parul Tewatia on 7 Feb 2019
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
  2 Comments
Jan
Jan on 7 Feb 2019
Edited: Jan on 7 Feb 2019
"I want a mech figure" is not clear. You have attached two images. Which one looks similar to what you want to achieve? Which are the input values you want to draw? What have you tried so far and which problem occur? As far as I can see, the posted code create some data, but there is no relation to the actual problem, is there? These details are more important that if you work with calcium or bananas.
Parul Tewatia
Parul Tewatia on 7 Feb 2019
Edited: Parul Tewatia on 7 Feb 2019
i have edited my question a bit, i had written clearly i have already made the figure spike_freq.jpg but now want the mesh.jpg figure. the code that is commented in the end is the one i tried to create the mesh but i have not been able to do that as i am bit unclear how to define the function in case of mesh figure. hope i am clear now.

Sign in to comment.

Accepted Answer

Florian Augustin
Florian Augustin on 7 Feb 2019
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
  1 Comment
Parul Tewatia
Parul Tewatia on 7 Feb 2019
thanks a lot Florian Augustin for the answer.
i will try it out and come back to you if i face anymore problems.
thanks a lot again...big help
parul

Sign in to comment.

More Answers (0)

Categories

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

Products

Community Treasure Hunt

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

Start Hunting!