Nonlinear 2D modeling of Neural Electromagnetic field in Matlab
I am trying to replicate the MATLAB simulation presented in this paper: http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=6399743. More specifically, I have to code the solution to this equation (http://ieeexplore.ieee.org/ielx5/6387371/6399636/6399743/html/img/6399743-eqdisp-7-small.png). Here m,n,and h are gating variables that are ordinary differential equations and are functions of time, and I have solved for them already. 394
I have also written the difference equation for this PDE. I'd like to set dE/dx(x,t) equal to 1 (denoted as g in the script) for now (and tackle that next). I'm using a finite differences method (I think)
clear all;
T=1; % Max time (s)
L=1; % Max length of neuron
dt=0.04; %time step
dx=0.01; %space step
t=[0:dt:T]; %time vector
x=[-L:dx:L]; %Space vector
M=length(x);
g=1; %The dE/dx bit
a=5; %Neuron radius
R_i=1; %Intracellular resistance
A=a./2*R_i;
I=0.1; %External Current Applied
ENa=55.17; % mv Na reversal potential
EK=-72.14; % mv K reversal potential
El=-49.42; % mv Leakage reversal potential
g_Na=1.2; % mS/cm^2 Na conductance
g_K=0.36; % mS/cm^2 K conductance
g_l=0.003; % mS/cm^2 Leakage conductance
b=4;
Cm=0.01; %Membrane capacitance
C=Cm;
%Solve for m,n,h
%Solve for m, n, h
V=-60; % Initial Membrane voltage
m=am(V)/(am(V)+bm(V)); % Initial m-value
n=an(V)/(an(V)+bn(V)); % Initial n-value
h=ah(V)/(ah(V)+bh(V)); % Initial h-value
y0=[V;n;m;h];
tspan = [0,max(t)];
%Matlab's ode45 function
[time,V] = ode45(@HH,tspan,y0);
V_hh=V(:,1);
n=V(:,2);
m=V(:,3);
h=V(:,4);
N=length(time);
plot(time,m,time,h,time,n);
V_initialgrid=zeros(M,N);
figure();
plot(time,V_hh)
V_initialgrid(1,:)=V_hh';
xlabel('time');
ylabel('length');
for j=2:N-1 %Stepping through time
B=g_Na*m(j).^3*h(j)+g_K*n(j).^4+g_l;
for i=2; %stepping thru space
V_initialgrid(i,j)=(A* V_initialgrid(i-1,j)+B* V_initialgrid(i,j)+C* V_initialgrid(i-1,j)+g);
for i=3:M-1 %stepping thru space
V_initialgrid(i+1,j)=(-dt*(g-(A.*(( V_initialgrid(i,j-1)-2.* V_initialgrid(i,j)+ V_initialgrid(i,j+1))./(dx.^2))-B.*-60))./C)+ V_initialgrid(i,j);
end
end
end
imagesc(log(abs(V_initialgrid)))
xlabel('Time');
ylabel('L');
The boundary conditions are as follows: dV/dx=dm/dx=dh/dx=dn/dx=0 when x=L or -L and dV/dt=dm/dt=dh/dt=dn/dt=0.
I know my error lies within the indexing of the solution of the difference equation. 268