Reproducing Neuron Channels in matlab
Posted: Wed Nov 30, 2011 4:43 am
I'm working on a script to export QUB (http://www.qub.buffalo.edu) models to NEURON and am having trouble getting the models to match. As an intermediate step I've also been trying to reproduce them in matlab, I can get the same traces in matlab as in QUB, but when I get to NEURON (particularly when I have more than 2 states) the traces don't match.
The instance I'm looking at now is a 3 state model with an SEClamp holding at -20 mV till 100ms, than jumping to 60 mV for 400 ms. The channel states probabilities end up quite divergent, the NEURON open channel hovers around .06 for the opening 100 ms (apparently I'm not setting the STEADYSTATE properly as it hasn't settled) and drops to .004 by the end of the simulation.
In matlab it opens at .2724 and drops to only .1583, I'm not sure why the NEURON trace should be so different.
thanks for any assistance.
This is the channel I create in NEURON
And here is the matlab code
The instance I'm looking at now is a 3 state model with an SEClamp holding at -20 mV till 100ms, than jumping to 60 mV for 400 ms. The channel states probabilities end up quite divergent, the NEURON open channel hovers around .06 for the opening 100 ms (apparently I'm not setting the STEADYSTATE properly as it hasn't settled) and drops to .004 by the end of the simulation.
In matlab it opens at .2724 and drops to only .1583, I'm not sure why the NEURON trace should be so different.
thanks for any assistance.
This is the channel I create in NEURON
Code: Select all
TITLE channel model imported from QUB
NEURON {
SUFFIX QUB_K_CHANNEL
USEION k READ ek WRITE ik
RANGE gbar, gk, ik, k0_C1_O1, k1_C1_O1, k0_O1_C1, k1_O1_C1, k0_C2_C1, k1_C2_C1, k0_C1_C2, k1_C1_C2
}
UNITS {
(mA) = (milliamp)
(mV) = (millivolt)
(pS) = (picosiemens)
(um) = (micron)
}
PARAMETER {
gbar = 0.001 (pS/um2)
k0_C1_O1 = 0.0005 (/ms)
k1_C1_O1 = 0.08 (/mV)
k0_O1_C1 = 0.0006 (/ms)
k1_O1_C1 = 0.04 (/mV)
k0_C2_C1 = 0.00081 (/ms)
k1_C2_C1 = 0.57 (/mV)
k0_C1_C2 = 0.002 (/ms)
k1_C1_C2 = 3 (/mV)
}
STATE {
O1 C2 C1 }
ASSIGNED {
v (mV)
ek (mV)
ik (mA/cm2)
gk (pS/um2)
C2C1 (/ms)
C1C2 (/ms)
C1O1 (/ms)
O1C1 (/ms)
}
INITIAL {
SOLVE states STEADYSTATE sparse
}
BREAKPOINT {
SOLVE states METHOD sparse
gk = gbar*O1
ik = (1e-4)*gk*(v - ek)
}
KINETIC states {
rates(v)
CONSERVE O1 + C2 + C1 = 1
~ C2 <-> C1 (C2C1, C1C2)
~ C1 <-> O1 (C1O1, O1C1)
}
PROCEDURE rates(v(millivolt)) {
C2C1 = k0_C2_C1 * exp(k1_C2_C1*v)
C1C2 = k0_C1_C2 * exp(k1_C1_C2*v)
C1O1 = k0_C1_O1 * exp(k1_C1_O1*v)
O1C1 = k0_O1_C1 * exp(k1_O1_C1*v)
}
Code: Select all
time_end = 500; % 500 ms
v_rev = 0;
V = zeros(time_end,1);
V(1:100) = -20;
V(101:time_end) = 60;
C = [0, .001, 0]; % conductance vector in pS/um2
K = zeros(3,3,2);
K(:,:,1) = [0, 0.0006, 0.00081; 0.0005, 0, 0; .002, 0, 0]; % the k0
K(:,:,2) = [0, 0.04, 0.57; 0.08, 0, 0; 3, 0, 0]; % k1
Q = zeros(size(K,2),size(K,2));
p = zeros(size(K,2),time_end+1);
I = zeros(time_end);
% calculate the inital Q matrix
for x=1:size(K,2)
for y=1:size(K,2)
Q(x,y) = K(x,y,1) * exp(K(x,y,2)*V(1)); % q_ij = k0*exp(k1*v)
end
end
for x=1:size(K,2) % set the diagonals
Q(x,x) = 0;
Q(x,x) = - sum(Q(:,x));
end
% solve for steady state
p(:,1) = [Q; ones(1, length(Q))] \ [zeros(length(Q), 1); 1];
for t=1:time_end
% generate the new transition matrix
for x=1:size(K,2)
for y=1:size(K,2)
Q(x,y) = K(x,y,1) * exp(K(x,y,2)*V(t)); % q_ij = k0*exp(k1*v)
end
end
for x=1:size(K,2) % set the diagonals
Q(x,x) = 0;
Q(x,x) = - sum(Q(:,x));
end
A = expm(Q);
p(:,t+1) = A*p(:,t);
p(:,t+1) = p(:,t+1)/sum(p(:,t+1)); % re-normalize
I(t) = 1e-4 * C * p(:,t+1) * (V(t)-v_rev); % current in mA/cm2
end
figure(1)
plot(V)
figure(2)
plot(I)
figure(3)
plot(p(2,:))