creating a mesh figure
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Parul Tewatia
el 7 de Feb. de 2019
Comentada: Parul Tewatia
el 7 de Feb. de 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 comentarios
Jan
el 7 de Feb. de 2019
Editada: Jan
el 7 de Feb. de 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.
Respuesta aceptada
Florian Augustin
el 7 de Feb. de 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
Más respuestas (0)
Ver también
Categorías
Más información sobre Biotech and Pharmaceutical en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!