creating a mesh figure

3 次查看（过去 30 天）
Parul Tewatia
Parul Tewatia 2019-2-7
评论： Parul Tewatia 2019-2-7
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 个评论
Jan
Jan 2019-2-7
编辑：Jan 2019-2-7
"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 2019-2-7
编辑：Parul Tewatia 2019-2-7
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Florian Augustin
Florian Augustin 2019-2-7
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 个评论
Parul Tewatia
Parul Tewatia 2019-2-7
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Biotech and Pharmaceutical 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by