hhred.m

利用Matlab模拟细胞HH 模型中的所有参数变化和动态Action Potential 生成过程

% 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 TNDHH

IEXHH = 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');
	end
end