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,:))
```