creating a mesh figure

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

[Jan]

"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]

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.

Answer 1

[Florian Augustin]

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

[Parul Tewatia]

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