function A2ONCB
%%% AII-ONCB network for a single AII and a single ONCB
%%% Cell 1 is AII; Cell
% Units: mF, mV,ms,S,ohm,mA, cm
clc
clear all
global numcell armsection R capacitance % Number of cells and compartments, R and capacitance
global gbar_na gbar_M gbar_A G_gap Gpas epas ena ek % conductances
global vhalfh_na vhalfm_na kh_na km_na % Na current activation/inactivation
global vhalfm_M km_M % M-type K (slow) current activation
global vhalfh_A vhalfm_A vhalfw_A kh_A km_A kw_A f % A-type K (fast) current activation/inactivation
global mtau_na htau_na mtau_M mtau_A % Time constants
global IinjectAII IinjectONCB % Injected current
global epas2
% set the values of injected currents with InjectVector (a list of values indicates switching to other values during the run)
% Note that the plots only include data after a transient given by cuttime.
% Set the conductance of the gap junction between AII and ONCB using the switch ONCB:
% ONCB=1 for Fig.14 and for Fig.5B top. ONCB=2,3,4 for remaining panels of Fig.5B
%%%%%%%%%% -------------------------- Switches -------------------------------------%%%%%%%%
plotswitch = 14; % 14 for Fig 14, 5 for Fig 5
scan = 0; % 0 if looking at voltage traces; 1 if scanning Vm (Fig2A); 2 if scanning VhlafM (Fig3B)
ONCB =1; % ONCB type (coupling strength)
RD = 1; % 1 if RD; 0 if wildtype
ONCBnum = 1;% 5; %
duration=1000; % duration of a run for a given parameter value
IinjectAII = 0;
IinjectONCB = 0;
% Current injection into ONCB
if plotswitch == 5
InjectVector = [0];
elseif plotswitch == 14
InjectVector = [0:2e-9:20*1e-9];
end
Vector = InjectVector;
tol_steady = 0.1; %(ms)
lookcell = 2; % 1 for cell 1 (AII), 2 for cell 2 (ONCB)
numcell = 2;
plotsection =1; %1 is Soma; 2 is cable; 3 is IS
v_init = -62;%(mV)
v_init2 = -50;%(mV)
% set coupling between AII and ONCB
if ONCB == 1
gjCondR= 750*1e-12 *ONCBnum; %(S)
elseif ONCB == 2
gjCondR= 300*1e-12 *ONCBnum; %(S)
elseif ONCB == 3;
gjCondR= 150*1e-12 *ONCBnum;
else
gjCondR= 100*1e-12 *ONCBnum; %(S)
end
gjCondIS= 1e-5*1e-12; %(S)
gjCondDC= 1e-5*1e-12; %(S)
%%%%%%%%%% -------------------------- Parameters Unchanged --------------------------------%%%%%%%%
armsection = 1; % # of sections in the cable, 11 in Mark's code
somadiam=25*1e-4; %(cm)
somalength=25*1e-4; %(cm)
armdiam=0.3*1e-4; %(cm)
armlength=32*1e-4; %(cm)
handdiam=2*1e-4; %(cm)
handlength=2*1e-4; %(cm)
SA_soma= 1963.5*1e-8; %(cm^2)
SA_arm=30.1593*1e-8; %(cm^2)
SA_hand=12.5664*1e-8; %(cm^2)
SA_ONCB= 439.8*1e-8; %(cm^2)
global_ra = 150; %(ohm-cm)
ena = 50; %(mV)
ek = -77; %(mV)
Cm = 1e-3; %(mf/cm^2)
%%
%%% --------- Na and K currents activation/inactivation -------- %%%
vhalfh_na = -49.5; %(mV)
vhalfm_na = -48; %(mV)
kh_na = 2; %(mV)
km_na = 5; %(mV)
mtau_na = 0.01; %(ms)
htau_na = 0.5; %(ms)
vhalfm_M = -40; %(mV)
km_M = 4; %(mV)
mtau_M = 50;%(ms)
vhalfh_A = -40.5; %(mV)
vhalfm_A = -10; %(mV)
vhalfw_A= -45; %(mV)
kh_A = 2; %(mV)
km_A = 7; %(mV)
kw_A = 15; %(mV)
f=0.83;
mtau_A = 1; %(ms)
gbar_na_soma = 0;
gbar_na_arm = 0;
gbar_na_IS = 0.2; %(mho/cm^2)
gbar_M_soma = 0;
gbar_M_arm = 0;
gbar_M_IS = 0.03; %(mho/cm2)
gbar_A_soma = 0.004; %(mho/cm2)
gbar_A_arm = 0;
gbar_A_IS = 0.08; %(mho/cm2)
if numcell == 2
gbar_na_soma2 = 0;
gbar_na_arm2 = 0;
gbar_na_IS2 = 0; %(mho/cm^2)
gbar_M_soma2 = 0;
gbar_M_arm2 = 0;
gbar_M_IS2 = 0; %(mho/cm2)
gbar_A_soma2 = 0; %(mho/cm2)
gbar_A_arm2 = 0;
gbar_A_IS2 = 0; %(mho/cm2)
end
%%%%%%%% ----------------------------- Leak current for AII ------------------------------------ %%%%%%%
if RD == 0
Rm = 40000; %(ohm*cm^2)
epas = -40; %(mV)
gpas=1/Rm; %(S/cm^2)
Rm2 = 12000;
epas2 = -30;%-50;
gpas2=1/Rm2;
elseif RD == 1
Rm = 40000; %(ohm*cm^2)
epas = -65;%-30; %(mV)
gpas=1/Rm; %(S/cm^2)
Rm2 = 12000;
epas2 = -35;%-50;
gpas2=1/Rm2;
end
%%%%%%%% -------------------------- Resistance between sections -------------------------------- %%%%%%%
R=zeros(armsection+1,1);
for num=1:armsection+1
if num==1
R(num)=(global_ra/(pi*(somadiam/2)^2))*somalength/2 + (global_ra/(pi*(armdiam/2)^2))*armlength/armsection/2;
elseif num==armsection+1
R(num)=(global_ra/(pi*(handdiam/2)^2))*handlength/2 + (global_ra/(pi*(armdiam/2)^2))*armlength/armsection/2;
else
R(num)=(global_ra/(pi*(armdiam/2)^2))*armlength/armsection;
end
end
%%%%%%%% ------------------------------ Conductances Set Up -------------------------%%%%%%%%
gbar_na=zeros(numcell,armsection+2);
gbar_M=zeros(numcell,armsection+2);
gbar_A=zeros(numcell,armsection+2);
capacitance=zeros(numcell,armsection+2);
G_gap=zeros(numcell,armsection+2);
Gpas=zeros(numcell,armsection+2);
for ii=1:numcell
if ii == 2
gbar_na_soma = gbar_na_soma2;
gbar_na_arm = gbar_na_arm2;
gbar_na_IS = gbar_na_IS2; %(mho/cm^2)
gbar_M_soma = gbar_M_soma2;
gbar_M_arm = gbar_M_arm2;
gbar_M_IS = gbar_M_IS2; %(mho/cm2)
gbar_A_soma = gbar_A_soma2; %(mho/cm2)
gbar_A_arm = gbar_A_arm2;
gbar_A_IS = gbar_A_IS2; %(mho/cm2)
gpas = gpas2;
SA_soma = SA_ONCB;
end
for jj=1:armsection+2
if jj==1
gbar_na(ii,jj) = gbar_na_soma * SA_soma;
gbar_M(ii,jj) = gbar_M_soma * SA_soma;
gbar_A(ii,jj)= gbar_A_soma * SA_soma;
capacitance(ii,jj) = Cm*SA_soma;
G_gap(ii,jj) = gjCondR;
Gpas(ii,jj) = gpas * SA_soma;
elseif jj==armsection+2
gbar_na(ii,jj) = gbar_na_IS * SA_hand;
gbar_M(ii,jj) = gbar_M_IS * SA_hand;
gbar_A(ii,jj) = gbar_A_IS * SA_hand;
capacitance(ii,jj) = Cm*SA_hand;
G_gap(ii,jj) = 0;%gjCondIS;
if ii == 1
Gpas(ii,jj) = gpas * SA_hand;
elseif ii == 2
Gpas(ii,jj) = 0;
end
else
gbar_na(ii,jj) = gbar_na_arm * SA_arm/armsection;
gbar_M(ii,jj) = gbar_M_arm * SA_arm/armsection;
gbar_A(ii,jj) = gbar_A_arm * SA_arm/armsection;
capacitance(ii,jj) = Cm*SA_arm/armsection;
G_gap(ii,jj) = 0;%gjCondDC;
if ii == 1
Gpas(ii,jj) = gpas * SA_arm/armsection;
elseif ii == 2
Gpas(ii,jj) = 0;
end
end
end
end
%%%%%% ------------------- Initial conditions --------------- %%%%%%%%%%
%%--------------- Initial conditions for Na
minit_na = 1/(1+exp(-(v_init-vhalfm_na)/km_na));
hinit_na = 1/(1+exp((v_init-vhalfh_na)/kh_na));
%%--------------- Initial conditions for M-type K
minit_M = 1/(1 + exp(-(v_init-vhalfm_M)/km_M));
%%--------------- Initial conditions for A-type K
minit_A = 1/(1+exp(-(v_init-vhalfm_A)/km_A));
hinit_A = f*(1/(1 + exp((v_init-vhalfh_A)/kh_A)))+(1-f);
%%%%%% ------------------- Initial conditions for cell 2 --------------- %%%%%%%%%%
if numcell == 2
minit_na2 = 1/(1+exp(-(v_init2-vhalfm_na)/km_na));
hinit_na2 = 1/(1+exp((v_init2-vhalfh_na)/kh_na));
minit_M2 = 1/(1 + exp(-(v_init2-vhalfm_M)/km_M));
minit_A2 = 1/(1+exp(-(v_init2-vhalfm_A)/km_A));
hinit_A2 = f*(1/(1 + exp((v_init2-vhalfh_A)/kh_A)))+(1-f);
end
%%%%%% ------------------- Initial conditions Set Up ------------------ %%%%%%%%%%
initial = zeros(numcell, armsection+2, 7);
for ii = 1:numcell
if ii==1
for jj = 1:armsection+2
initial(ii,jj,1) = minit_na;
initial(ii,jj,2) = hinit_na;
initial(ii,jj,3) = minit_M;
initial(ii,jj,4) = minit_A;
initial(ii,jj,5) = hinit_A;
initial(ii,jj,6) = hinit_A;
initial(ii,jj,7) = v_init;
end
elseif ii==2
for jj = 1:armsection+2
initial(ii,jj,1) = minit_na2;
initial(ii,jj,2) = hinit_na2;
initial(ii,jj,3) = minit_M2;
initial(ii,jj,4) = minit_A2;
initial(ii,jj,5) = hinit_A2;
initial(ii,jj,6) = hinit_A2;
initial(ii,jj,7) = v_init2;
end
end
end
initial = reshape(initial,numcell*(armsection+2)*7,1);
vmaxmat = zeros(1,length(Vector ));%[];
vminmat = zeros(1,length(Vector ));%[];
frequencymat = zeros(1,length(Vector ));%[];
ampmat = zeros(1,length(Vector));
Tend = 0;
for ii = 1:length(InjectVector)
IinjectONCB = InjectVector(ii);
time_start = 0;
time_end = duration;%000;
limit_time = 1100;
cuttime = 500;
freq = 0;
amp = 0;
while time_end < limit_time
%%%%%%%%%%% ------------------- Solving Diff Eqs ------------------------ %%%%%%%%%%
options = odeset('RelTol', 1e-8);
[T,Y1] = ode45(@solve_A2ONCB, [time_start time_end], initial,options);
Y = zeros(numcell,armsection+2,7,length(T));
for i = 1:length(T)
Y(1:numcell,1:armsection+2,1:7,i) = reshape(Y1(i,:),numcell, armsection+2,7);
end
initial_new = reshape(Y(:,:,:,end),numcell*(armsection+2)*7,1);
initial = initial_new;
plotcell =1;
v= Y(plotcell,plotsection,7,:);
v_fix=reshape(v,size(T));
vIS= Y(plotcell,3,7,:);
vIS_fix=reshape(vIS,size(T));
%%%%% -------------- If there are more than one cell -------------- %%%%%%%
if numcell == 2
plotcell2 = 2;
v2= Y(plotcell2,plotsection,7,:);
v2_fix=reshape(v2,size(T));
end
v_AII = v_fix;%[v_AII;v_fix];
v_ONCB = v2_fix;%[v_ONCB;v2_fix];
%%% ----------------- Cutting of the transient ------------------%%%
savenow = [];
for i = 1:length(T)
if T(i) >= cuttime
savenow = [savenow i];
else
savenow = savenow;
end
start = min(savenow);
end
if lookcell == 1
v = v_fix(start:end);
elseif lookcell == 2
v = v2_fix(start:end);
end
t = T(start:end);
%%% ----------------- Finding local max/minima ----------------------%%%
if lookcell == 1
[vmax,imax,vmin,imin] = extrema(v_fix(start:end));
elseif lookcell == 2
[vmax,imax,vmin,imin] = extrema(v2_fix(start:end));
end
imax = sort(imax);
imin = sort(imin);
%% No spiking nor bursting
if length(vmax)<=2 || length(vmax)<=2
time_end = time_end+500;
time_start = time_start+500;
cuttime = time_start;
if time_end>=limit_time
tmax = T(imax);
Vmax = vmax;
tmin = T(imin);
Vmin = vmin;
if length(vmax)<1 || length(vmax)<1
freq = 0; %(Hz)
amp = 0;
t_cross = [];
else
freq = 0; %(Hz)
amp = max(Vmax)-min(Vmin);
t_cross = [];
end
break
end
else
%% Quadratic interpolation
tmax = zeros(1,length(imax));
Vmax = zeros(1,length(imax));
iMax = zeros(1,length(imax));
for index = 1:length(imax)
iMax = imax(index);
x = [t(iMax-1) t(iMax) t(iMax+1)];
y = [v(iMax-1) v(iMax) v(iMax+1)];
[p,S,mu] = polyfit(x,y,2);
A_new = p(1); B_new = p(2); C_new = p(3); mu1 = mu(1); mu2 = mu(2);
c=A_new*(mu1)^2/(mu2)^2-(B_new*(mu1/mu2))+C_new;
b = (B_new*mu2 - 2*A_new*mu1)/(mu2)^2;
a = A_new/(mu2)^2;
tmax(index) = -b/(2*a);
Vmax(index) = c-b^2/(4*a);
end
tmin = zeros(1,length(imin));
Vmin = zeros(1,length(imin));
iMin = zeros(1,length(imin));
for index = 1:length(imin)
iMin = imin(index);
x = [t(iMin-1) t(iMin) t(iMin+1)];
y = [v(iMin-1) v(iMin) v(iMin+1)];
[p,S,mu] = polyfit(x,y,2);
A_new = p(1); B_new = p(2); C_new = p(3); mu1 = mu(1); mu2 = mu(2);
c=A_new*(mu1)^2/(mu2)^2-(B_new*(mu1/mu2))+C_new;
b = (B_new*mu2 - 2*A_new*mu1)/(mu2)^2;
a = A_new/(mu2)^2;
tmin(index) = -b/(2*a);
Vmin(index) = c-b^2/(4*a);
end
%%
crossline = (max(Vmax)+min(Vmin))/2;
count1 = 1;
count2 = 1;
Vmin_re = [];
while count1 <= length(tmin)
if Vmin(count1)<crossline
Vmin_re(count2) = Vmin(count1);
count2=count2+1;
end
count1=count1+1;
end
count1 = 1;
count2 = 1;
Vmax_re = [];
while count1 <= length(tmax)
if Vmax(count1)>crossline
Vmax_re(count2) = Vmax(count1);
count2=count2+1;
end
count1=count1+1;
end
Vmax = Vmax_re;
Vmin = Vmin_re;
%%
t_cross = [];
v_cross = (max(vmax)+min(vmin))/2;
for iter = 1:length(v)-1
if (v(iter)-v_cross)<=0 && (v(iter+1)-v_cross)>=0
t_cross = [t_cross ((v_cross-v(iter))*(t(iter+1)-t(iter))/(v(iter+1)-v(iter)))+t(iter)];
else
t_cross = t_cross;
end
end
V_cross = ones(1,length(t_cross))*v_cross;
if length(t_cross)<=2
time_end = time_end+500;
time_start = time_start+500;
cuttime = time_start;
% cuttime = cuttime+500;
if time_end>=limit_time
if length(t_cross)<=1
freq = 0; %(Hz)
amp = Vmax(end)-Vmin(end);
else
period = (t_cross(2)-t_cross(1))*1e-3;
freq=1/period;
amp = Vmax(end)-Vmin(end);
end
break
end
else
dtcross = zeros(1,length(t_cross)-1);
for step = 1:length(t_cross)-1
dtcross(step) = t_cross(step+1)-t_cross(step);
end
vari = zeros(1,length(dtcross)-1);
for step = 1:length(dtcross)-1
vari(step) = abs(dtcross(step+1)-dtcross(step));
if vari(step) < tol_steady
% cuttime = tmax(step);
period = dtcross(step+1)*1e-3; %(s)
freq = 1/period; %(Hz)
amp = Vmax(step+1)-Vmin(step+1);
break
end
end
%% Increase the simulataion time if :
% 1) it does not go to a steady state
% 2) it does not have more than one period in it
% Need a limit
if min(vari)>=tol_steady
time_end = time_end+500;
time_start = time_start+500;
cuttime = time_start;
if time_end>=limit_time
period = dtcross(end)*1e-3; %(s)
freq = 1/period; %(Hz)
amp = Vmax(end)-Vmin(end);
break
end
else
break
end
end
end
end
frequencymat(ii) = freq;
ampmat(ii) = amp;
%%%%%%%%%%------------------------------------ Decide on Plots -------------------------------------%%%%%%%%%%
if plotswitch == 14
figure(1)
hold on
subplot(3,1,1), hold on, plot(T+Tend,IinjectONCB*ones(1,length(T)))
xlabel('Time (ms)','FontWeight','bold'), ylabel('I_{inject} (mA)','FontWeight','bold')
subplot(3,1,2), hold on, plot(T+Tend,v_AII)
xlabel('Time (ms)','FontWeight','bold'), ylabel('V_{AII} (mV)','FontWeight','bold')
subplot(3,1,3), hold on, plot(T+Tend,v_ONCB)
xlabel('Time (ms)','FontWeight','bold'), ylabel('V_{ONCB} (mV)','FontWeight','bold')
elseif plotswitch == 5
figure(2)
hold on, plot(T+Tend,v_AII,'b', T+Tend,v_ONCB,'r'), ylim([-68 -40]), xlim([0 900])
xlabel('Time (ms)'); ylabel(' Voltage (mV)');
legend('AII','ONCB')
end
Tend=Tend+T(end);
end
end
%%%%%%%%%%%%----------- MATLAB reconstruction of Mark's AII model on NEURON --------------%%%%%%%%%%
% Units: mF, mV,ms,S,ohm,mA, cm
function dy = solve_A2ONCB(t,y1)
global numcell armsection R capacitance % Number of cells and compartments, R and capacitance
global gbar_na gbar_M gbar_A G_gap Gpas epas ena ek % conductances
global vhalfh_na vhalfm_na kh_na km_na % Na current activation/inactivation
global vhalfm_M km_M % M-type K (slow) current activation
global vhalfh_A vhalfm_A vhalfw_A kh_A km_A kw_A f % A-type K (fast) current activation/inactivation
global mtau_na htau_na mtau_M mtau_A % Time constants
global IinjectAII IinjectONCB
global epas2
% global inject_start inject_end
y = reshape(y1,numcell,armsection+2,7);
dy = zeros(numcell,armsection+2,7);
%%%%%%%%%%%% ---------------------------- Differential Equations ------------------------- %%%%%%%%%%%
for i=1:numcell
for j=1:armsection+2
if i ==1 && j == 1% && t>=inject_start && t<inject_end% i == 17
ie = IinjectAII;
elseif i == 2 && j == 1
ie = IinjectONCB;
else
ie = 0;
end
dy(i,j,1) = (1/mtau_na)*(1/(1+exp(-(y(i,j,7)-vhalfm_na)/km_na))-y(i,j,1)); % m_na
dy(i,j,2) = (1/htau_na)*(1/(1+exp((y(i,j,7)-vhalfh_na)/kh_na))-y(i,j,2)); %h_na
dy(i,j,3) = (1/mtau_M)*(1/(1 + exp(-(y(i,j,7)-vhalfm_M)/km_M))-y(i,j,3)); %m_M
dy(i,j,4) = (1/mtau_A)*(1/(1 + exp(-(y(i,j,7)-vhalfm_A)/km_A))-y(i,j,4)); %m_A
dy(i,j,5) = (1/(-20/(1+exp(-(y(i,j,7)+35)/6))+25))*(f*(1/(1 + exp((y(i,j,7)-vhalfh_A)/kh_A)))+(1-f)-y(i,j,5)); %h0_A
dy(i,j,6) = (1/max(100,((y(i,j,7)+17)^2/4+26)))*(f*(1/(1 + exp((y(i,j,7)-vhalfh_A)/kh_A)))+(1-f)-y(i,j,6)); %h1_A
if i == 1
if j==1
isection = 1/R(j)*(y(i,j+1,7)-y(i,j,7));
elseif j==armsection+2
isection = 1/R(j-1)*(y(i,j-1,7)-y(i,j,7));
else
isection = 1/R(j)*(y(i,j+1,7)-y(i,j,7))+1/R(j-1)*(y(i,j-1,7)-y(i,j,7));
end
elseif i == 2
isection = 0;
end
if numcell == 1
igap = 0;
else
if i==1
igap = G_gap(i,j)*(y(i+1,j,7)-y(i,j,7));
elseif i==numcell
igap = G_gap(i,j)*(y(i-1,j,7)-y(i,j,7));
else
igap = G_gap(i,j)*(y(i+1,j,7)-y(i,j,7))+G_gap(i,j)*(y(i-1,j,7)-y(i,j,7));
end
end
if i==2
im = gbar_na(i,j)*y(i,j,1).^3.*y(i,j,2).*(y(i,j,7) - ena)+gbar_M(i,j)*y(i,j,3)*(y(i,j,7) - ek)+1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A))*gbar_A(i,j)*y(i,j,4)*y(i,j,5)*(y(i,j,7) - ek) + (1-1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A)))*gbar_A(i,j)*y(i,j,4)*y(i,j,6)*(y(i,j,7) - ek)+Gpas(i,j)*(y(i,j,7)-epas2);
elseif i == 1
im = gbar_na(i,j)*y(i,j,1).^3.*y(i,j,2).*(y(i,j,7) - ena)+gbar_M(i,j)*y(i,j,3)*(y(i,j,7) - ek)+1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A))*gbar_A(i,j)*y(i,j,4)*y(i,j,5)*(y(i,j,7) - ek) + (1-1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A)))*gbar_A(i,j)*y(i,j,4)*y(i,j,6)*(y(i,j,7) - ek)+Gpas(i,j)*(y(i,j,7)-epas);
end
dy(i,j,7) = (1/capacitance(i,j))*(-im + isection + igap +ie);
end
end
dy = reshape(dy,numcell*(armsection+2)*7,1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xmax,imax,xmin,imin] = extrema(x)
%EXTREMA Gets the global extrema points from a time series.
% [XMAX,IMAX,XMIN,IMIN] = EXTREMA(X) returns the global minima and maxima
% points of the vector X ignoring NaN's, where
% XMAX - maxima points in descending order
% IMAX - indexes of the XMAX
% XMIN - minima points in descending order
% IMIN - indexes of the XMIN
%
% DEFINITION (from http://en.wikipedia.org/wiki/Maxima_and_minima):
% In mathematics, maxima and minima, also known as extrema, are points in
% the domain of a function at which the function takes a largest value
% (maximum) or smallest value (minimum), either within a given
% neighbourhood (local extrema) or on the function domain in its entirety
% (global extrema).
%
% Example:
% x = 2*pi*linspace(-1,1);
% y = cos(x) - 0.5 + 0.5*rand(size(x)); y(40:45) = 1.85; y(50:53)=NaN;
% [ymax,imax,ymin,imin] = extrema(y);
% plot(x,y,x(imax),ymax,'g.',x(imin),ymin,'r.')
%
% See also EXTREMA2, MAX, MIN
% Written by
% Lic. on Physics Carlos Adrián Vargas Aguilera
% Physical Oceanography MS candidate
% UNIVERSIDAD DE GUADALAJARA
% Mexico, 2004
%
% nubeobscura@hotmail.com
% From : http://www.mathworks.com/matlabcentral/fileexchange
% File ID : 12275
% Submited at: 2006-09-14
% 2006-11-11 : English translation from spanish.
% 2006-11-17 : Accept NaN's.
% 2007-04-09 : Change name to MAXIMA, and definition added.
xmax = [];
imax = [];
xmin = [];
imin = [];
% Vector input?
Nt = numel(x);
if Nt ~= length(x)
error('Entry must be a vector.')
end
% NaN's:
inan = find(isnan(x));
indx = 1:Nt;
if ~isempty(inan)
indx(inan) = [];
x(inan) = [];
Nt = length(x);
end
% Difference between subsequent elements:
dx = diff(x);
% Is an horizontal line?
if ~any(dx)
return
end
% Flat peaks? Put the middle element:
a = find(dx~=0); % Indexes where x changes
lm = find(diff(a)~=1) + 1; % Indexes where a do not changes
d = a(lm) - a(lm-1); % Number of elements in the flat peak
a(lm) = a(lm) - floor(d/2); % Save middle elements
a(end+1) = Nt;
% Peaks?
xa = x(a); % Serie without flat peaks
b = (diff(xa) > 0); % 1 => positive slopes (minima begin)
% 0 => negative slopes (maxima begin)
xb = diff(b); % -1 => maxima indexes (but one)
% +1 => minima indexes (but one)
imax = find(xb == -1) + 1; % maxima indexes
imin = find(xb == +1) + 1; % minima indexes
imax = a(imax);
imin = a(imin);
nmaxi = length(imax);
nmini = length(imin);
%%%% ----- Commented Out on 052012 by Hannah ------------------%%%%%
% Maximum or minumim on a flat peak at the ends?
if (nmaxi==0) && (nmini==0)
if x(1) > x(Nt)
xmax = x(1);
imax = indx(1);
xmin = x(Nt);
imin = indx(Nt);
elseif x(1) < x(Nt)
xmax = x(Nt);
imax = indx(Nt);
xmin = x(1);
imin = indx(1);
end
return
end
%% Maximum or minumim at the ends?
% if (nmaxi==0)
% imax(1:2) = [1 Nt];
% elseif (nmini==0)
% imin(1:2) = [1 Nt];
% else
% if imax(1) < imin(1)
% imin(2:nmini+1) = imin;
% imin(1) = 1;
% else
% imax(2:nmaxi+1) = imax;
% imax(1) = 1;
% end
% if imax(end) > imin(end)
% imin(end+1) = Nt;
% else
% imax(end+1) = Nt;
% end
% end
%%%%% -----------------------------------------------------------%%%%%
xmax = x(imax);
xmin = x(imin);
% NaN's:
if ~isempty(inan)
imax = indx(imax);
imin = indx(imin);
end
% Same size as x:
imax = reshape(imax,size(xmax));
imin = reshape(imin,size(xmin));
% Descending order:
[temp,inmax] = sort(-xmax); clear temp
xmax = xmax(inmax);
imax = imax(inmax);
[xmin,inmin] = sort(xmin);
imin = imin(inmin);
% Carlos Adrián Vargas Aguilera. nubeobscura@hotmail.com
end
