runModel.m

function [firing_rate] = runModel(g)
% this function will simulate the signaling-electrophysiology model
% the ode's for the model are in the file odemodel.m
%%%% input: g, a vector of 6 values weighting the reference levels of AT1R
%%%%%%% and five channels: gNa, gKdr, gKa, gKahp, gCaL
%%%%%%% code to run reference simulation:
%%%%%%% [firing_rate] = runModel([1,1,1,1,1,1]); 
%%%% output: firing rate (firing_rate)

options = odeset( 'RelTol', 1e-9, 'AbsTol', 1e-9 );

% simulation time parameters
ti      = -1;               % initial simulation time (s)
dt      = 0.5*10^(-3);      % time step (0.0005 s)
tf      = 300;              % final simulation time (s)
tspan   = ti:dt:tf;         % simulation time span (s)

% load initial conditons
y0      = LoadInitialConditions;

% modify AT1R level according to input g(1)
% these parameters are part of the AT1R complex
y0(94)  = y0(94)*g(1);
y0(116) = y0(116)*g(1);
y0(117) = y0(117)*g(1);
y0(118) = y0(118)*g(1);
y0      = abs(y0);

% load model parameters
% input g is supplied to LoadParameters so that conductances are modified
params  = LoadParameters(g);

% AngII stimulus, applied continuously starting at t = 0
AngII100 = [ti:0.01:(tf+1); ti:0.01:(tf+1)]';
AngII100(AngII100(:,1)<0,2) = 0;
AngII100(AngII100(:,1)>0,2) = 0.1;

% run simuations
[t, y] = ode15s(@odemodel, tspan, y0, options, params, AngII100);

% determine firing rate
[tp_fr fr index]    = freq(t,y(:,179));
firing_rate         = interp1(tp_fr, fr,tspan);

%%%%% error checking %%%%%
% replace firing rate output with zeros simulation errors occur

%%% check the length of the firing rate vector
%%% if too short, replace with zeros
if (length(firing_rate) < length(tspan))
    firing_rate = zeros(length(tspan),1);
end

%%% check for nan's and replace with zeros
indNAN = find(isnan(firing_rate));
if (isempty(indNAN)>0)
    firing_rate(indNAN) = 0;
end

figure;
plot(tspan,firing_rate)
xlabel('time (s)');
ylabel('firing rate (Hz)')

end

function [tp, f, index] = freq(t,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to calculate firing frequency
%
%  inputs: t - time vector
%          x - oscillating function (indices should match t)
%
%  outputs: tp - time points 
%           f  - corresponding frequency (at that time point)
%           index - indices of the time vector
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=1;
index=1;
tp=t(1);

for i=1:length(t)-1
    if (x(i)<-0&x(i+1)>-0), 
        tp(k)=t(i); 
        k=k+1; 
        index = i;    
    end
end

try
    f=[1./diff(tp)];
    f=[f, f(end)];
    tp=[tp];
catch
     tp=[t(1) t(end)];
     f=[0, 0];
end

end