7 views (last 30 days)

Show older comments

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

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

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

Start Hunting!