% HHRED is a script that does simulation of reduced, 2 state-variable, HH equations from% resting initial conditions with current pulse injection. Uses hhodered.m.% The state variables are [V, n]. m = minf(V) and h = nvtohv(n).% Globals are parameters of current (see iexthh())global IEXHH TSTHH TNDHHIEXHH = 0; TSTHH = 0; TNDHH = 50;% HH channel parameters (mS/cm^2, mV, uFd/cm^2). NOTE these parameters% should not be set by constants anywhere! SETHHPARMS and GETHHPARMS should% always be used.% gnamx=120 ; Ena=55 ; gkmx=36 ; Ek=-72 ; gleak=0.3 ; El=-49.4 ; C=1 ;gnamx=120 ; Ena=55 ; gkmx=36 ; Ek=-72 ; gleak=0.3 ; El=-49.4 ; C=1 ;sethhparms([gnamx; Ena; gkmx; Ek; gleak; El; C]);gunk = 'X';while gunk~='Q' & gunk~='q'	% Get current	IEXHH = input('Current (uA/cm^2): ');	if IEXHH~=0		TSTHH = input('Starting at ?? ms: ');		TNDHH = input(' . . and ending at ?? ms: ');	end	% and duration	tmax = input('Max simulation time (ms): ');	% Get default I.C.s, shall we change them?	[tspan, yv0, opt1] = hhodered(0,0,'init'); tspan(2) = tmax;	opt = odeset(opt1, 'AbsTol', [1e-4; 1e-6], 'RelTol', 1e-5);	yvinit = input(sprintf( ...'Current I.C.s: V=%g;  n=%g.\nEnter new ICs, in order (1-2 as ''-60.1 0.4''): ', ...        yv0),'s');    if ~isempty(yvinit)        [newics, count] = sscanf(yvinit,'%g %g',2);        yv0(1:count)=newics(1:count);    end	fprintf('ICs to be used: V=%g, n=%g.\n', yv0)	% Do simulation and print final stateyv	clf;	[tim, yv] = ode15s('hhodered', tspan, yv0, opt);	fprintf('Final state: V=%g, n=%g.\n',yv(size(yv,1),:));	gunk = 'V';	% Display loop	while gunk~='Q' & gunk~='q' & gunk~='A' & gunk~='a'		if gunk~='q' & gunk ~='Q' & gunk~='a' & gunk~='A'			subplot(2,1,1); cla; legend off;			subplot(2,1,2); cla; legend off;% Plot state variables			subplot(2,1,1);plot(tim, yv(:,1)); legend('V');			thax=axis; axis([thax(1), thax(2), -80, 40]);			title(sprintf('Hodgkin Huxley reduced model, Iext=%g.',IEXHH))			ylabel('Membrane pot., mV'); hold on			subplot(2,1,2);plot(tim, yv(:,2));            hold on; plot(tim, minfhh(yv(:,1)),'g'); plot(tim, nvtohv(yv(:,2)),'r'); hold off            legend('n','m','h')			ylabel('HH variables')			xlabel('Milliseconds')		end				if gunk=='V' | gunk=='v'			subplot(2,1,2); plot(tim, yv(:,2));            hold on; plot(tim, minfhh(yv(:,1)),'g'); plot(tim, nvtohv(yv(:,2)),'r'); hold off            legend('n','m','h')			ylabel('HH variables')		elseif gunk=='G' | gunk=='g'			subplot(2,1,2); cla; plot(tim, gkmx*yv(:,2).^4,'b'); hold on			plot(tim, gnamx*minfhh(yv(:,1)).^3.*nvtohv(yv(:,2)),'c'); legend('Gk', 'Gna')			ylabel('Conductances, mS/cm^2')		elseif gunk=='I' | gunk=='i'			subplot(2,1,2); cla;			IK = gkmx*yv(:,2).^4.*(yv(:,1)-Ek); plot(tim, IK, 'b'); hold on			INa = gnamx*minfhh(yv(:,1)).^3.*nvtohv(yv(:,2)).*(yv(:,1)-Ena); plot(tim, INa,'c')			npt = size(yv,1); Icap = C*((yv(2:npt,1)-yv(1:npt-1,1))./ ...			    (tim(2:npt)-tim(1:npt-1))); plot(tim(1:npt-1), Icap, 'k')			Ileak = gleak*(yv(:,1)-El); plot(tim, Ileak,'g')			legend('Ik', 'Ina', 'Icap', 'Ileak')			ylabel('Current, uA/cm^2')		end		xlabel('Milliseconds')		gunk = input('Again, plot Vars, G-conductances, I-currents, Quit? ','s');	endend