% This file  looks at two coupled MS cells 
% as a model of border zone arrhythmia.


clear
time=cputime;

% parameter values
gna = 3.33;
gk = 0.2;
%gk = 3.16;
vgate = 0.13;
kappa = 20;
topen = 120;
tclose = 150;

vna = 1;
vk = .12; %gives a "dead" cell
%vk = .06; % gives a modestly ischemic cell
%vk = 0.02;
%oscillations occur with coup =0.03, but not with coup = 0;
coup = .15;
coup = 0.05  %this is the coupling strength, with this = 0 there is no coupling
coup = 0
M = 500;
kstep = 20;
dt = 0.2;
N=100;

%phase plane
Vt=(1:N-1)/N;
V1n = gk./(gna*Vt.*(vna-Vt));
V2n = gk.*(Vt-vk)./(gna*Vt.^2.*(vna-Vt));
Wn = 1/2-tanh(kappa*(Vt-vgate))/2;

V1=0; 
W1 = .9945;
v1(1,1)=V1;
w1(1,1)=W1;
V2=0.177; 
W2 = .133;
v2(1,1)=V2;
w2(1,1)=W2;
Tt(1,1)=0;



t=0;

% The algorithm  (here the integration uses a simple Euler step method)

for n=2:M
for k = 1:kstep
fv1= 1/2+tanh(kappa*(V1-vgate))/2;
fv2 = 1/2+tanh(kappa*(V2-vgate))/2;
th1 = topen +(tclose-topen)*fv1;
th2 = topen +(tclose-topen)*fv2;
hinf1 = 1-fv1;
hinf2 = 1-fv2;

G1 = (hinf1-W1)/th1;
G2 = (hinf2-W2)/th2;


% update W
W1 = W1 + dt*G1;
W2 = W2 + dt*G2;
      
m1 = min(max(V1,0),1);
m2 = min(max(V2,0),1);
	  
	  
F1  = gna*W1.*m1.^2.*(vna-V1)-gk*V1 +coup*(V2-V1); %This is the "normal cell
      
F2  = gna*W2.*m2.^2.*(vna-V2)+gk*(vk-V2) +coup*(V1-V2)/2;  %this is the ischemic cell

%update V
V1 = V1+dt*F1;
V2=V2+dt*F2;

t=t + dt;
v1(n,1)=V1;
w1(n,1)=W1;
v2(n,1)=V2;
w2(n,1)=W2;
Tt(n,1)=t;
end
end
Tl = t;
%for n=1:M
figure(1);
title(['coupled cells']);

  subplot(2,2,1)
  plot(Tt,v1,'b',Tt,w1,'g')
  title(['normal cell']);
   axis([0 Tl 0 1])
  subplot(2,2,2)
  plot(Tt,v2,'b',Tt,w2,'g')
  title(['"ischemic" cell']);
  axis([0 Tl 0 1])
  subplot(2,2,3)
  plot(Vt,V1n,'g--',Vt,Wn,'r--',v1,w1,'b')
  axis([0 1 0 1])
  
   subplot(2,2,4)
  plot(Vt,V2n,'g--',Vt,Wn,'r--',v2,w2,'b')
  axis([0 1 0 1])

  pause(0.05)

%end

TL = t;
time=cputime-time        % outputs the cpu time taken to solve the equations

time_steps = M;

