NEURON + Python Michael Hines HBP CodeJam Workshop #7 Manchester - PowerPoint PPT Presentation

NEURON + Python Michael Hines HBP CodeJam Workshop #7 Manchester 2016 NINDS

I r n e t t e e r r p p r r e e t Python t e n HOC r I Neuron specific syntax Compiled

I r n e t t e e r r p p r r e e t Python t e n HOC r I Neuron specific syntax Compiled Run: compute model

I r n e t t e e r r p p r r e e t Python t e n HOC r I Neuron specific syntax LFP Compiled RxD Run: compute model

Gather values into numpy array. Interpreter Python Compiled

Gather values into numpy array. Interpreter Python from neuron import h import numpy Compiled

Gather values into numpy array. Interpreter Python from neuron import h import numpy hv = h.Vector(size) v = hv.as_numpy() Compiled

Gather values into numpy array. Interpreter Python from neuron import h import numpy hv = h.Vector(size) v = hv.as_numpy() pv = h.PtrVector(size) pv.pset(i, _ref_hocvar) pv.gather(hv) Compiled

Gather values into numpy array. Interpreter Python from neuron import h import numpy hv = h.Vector(size) v = hv.as_numpy() pv = h.PtrVector(size) pv.pset(i, _ref_hocvar) pv.gather(hv) i=0 for sec in h.allsec(): Compiled for seg in sec: pv.pset(i, seg._ref_v) i += 1

Python API for NEURON solvers

Python API for NEURON solvers nonvint_block_supervisor dy/dt = f(y, v, t) y can affect channel conductance

Python API for NEURON solvers nonvint_block_supervisor dy/dt = f(y, v, t) y can affect channel conductance High performance evaluation of M*x = b ∂ ∂ f / y M ~ 1 − dt*

NMODL ica Pump [Ca] t n e c a j d a e n a r b m e M Bulk KINETIC pmp { ~ cabulk <−> cai (1/tau, 1/tau) ~ cai + pump <−> capump (k1, k2) ~ capump <−> cao + pump (k3, k4) ica_pmp = 2*FARADAY*(f_flux − b_flux) ~ cai << −(ica) : there is a problem here COMPARTMENT width {cai} : volume has dimensions of (um) COMPARTMENT 1 {pump capump}�: area is dimensionless COMPARTMENT 1(m) {cao cabulk} } 0.1 0.003 soma.cai( 0.5 ) mM mA/cm2 0.005 soma.ica( 0.5 ) 0.08 0.002 0.06 0.04 0 0.001 65 73 0.02 0 0 0 20 40 60 80 100 0 20 40 60 80 100 ms ms

NMODL ica Pump [Ca] t n e c a j d a e n a r b m e M Bulk KINETIC pmp { ~ cabulk <−> cai (1/tau, 1/tau) ~ cai + pump <−> capump (k1, k2) ~ capump <−> cao + pump (k3, k4) ica_pmp = 2*FARADAY*(f_flux − b_flux) ~ cai << −(ica) : there is a problem here COMPARTMENT width {cai} : volume has dimensions of (um) COMPARTMENT 1 {pump capump}�: area is dimensionless COMPARTMENT 1(m) {cao cabulk} } 0.1 0.003 soma.cai( 0.5 ) mM mA/cm2 0.005 soma.ica( 0.5 ) 0.08 0.002 0.06 0.04 0 0.001 65 73 0.02 0 0 0 20 40 60 80 100 0 20 40 60 80 100 ms ms

from neuron import nonvint_block_supervisor as nbs callables = [ setup, initialize, current, conductance, fixed_step_solve, count, reinit, fun, msolve, jacobian, abs_tolerance ] nbs.register(callables)

from neuron import nonvint_block_supervisor as nbs callables = [ setup, initialize, current, conductance, fixed_step_solve, count, reinit, fun, msolve, jacobian, abs_tolerance ] class ExtraEquations(): nbs.register(callables) def __init__(self): self.callables = [...] nbs.register(self.callables) ... def __del__(self): nbs.unregister(self.callables) ...

from neuron import nonvint_block_supervisor as nbs callables = [ setup, initialize, Contribution to current balance eqns current, conductance, fixed_step_solve, count, reinit, fun, msolve, jacobian, abs_tolerance ] nbs.register(callables)

from neuron import nonvint_block_supervisor as nbs callables = [ setup, initialize, current, conductance, y(t) −> y(t + dt) fixed_step_solve, count, reinit, fun, msolve, jacobian, abs_tolerance ] nbs.register(callables)

from neuron import nonvint_block_supervisor as nbs callables = [ setup, initialize, current, conductance, fixed_step_solve, Additional equations for CVODE count, reinit, fun, msolve, jacobian, abs_tolerance ] nbs.register(callables)

Contribution to current balance eqns. Subtract current from rhs numpy array. def current(self, rhs): self.pv.gather(self.hv) n = self.n self.g = self.gkbar*n*n*n*n i = self.g*(self.v + 77.) rhs[self.nodeindices] −= i def conductance(self, d): d[self.nodeindices] += self.g Add conductance to matrix diagonal numpy array.

Contribution to current balance eqns. def current(self, rhs): self.pv.gather(self.hv) n = self.n self.g = self.gkbar*n*n*n*n i = self.g*(self.v + 77.) rhs[self.nodeindices] −= i def conductance(self, d): Copy seg.node_index() during setup. d[self.nodeindices] += self.g

Contribution to current balance eqns. def current(self, rhs): Voltage needed ... self.pv.gather(self.hv) n = self.n self.g = self.gkbar*n*n*n*n ... to compute current i = self.g*(self.v + 77.) rhs[self.nodeindices] −= i def conductance(self, d): d[self.nodeindices] += self.g

Contribution to current balance eqns. def current(self, rhs): self.pv.gather(self.hv) Instance gating states... n = self.n ... needed for current as well self.g = self.gkbar*n*n*n*n i = self.g*(self.v + 77.) rhs[self.nodeindices] −= i def conductance(self, d): d[self.nodeindices] += self.g

CVODE interface Where our portion of the CVODE def count(self, offset): state vector begins. self.offset = offset return len(self.v) def fun(self, t, y, ydot): last = self.offset + len(self.n) self.pv.gather(self.hv) self.n = y[self.offset:last] if ydot == None: return ninf, nrate = self.ninftau(self.v) ydot[self.offset:last] = (ninf − self.n)*nrate def msolve(self, dt, t, b, y): last = self.offset + len(self.n) self.pv.gather(self.hv) x,nrate = self.ninftau(self.v) b[self.offset:last] /= 1. + dt*nrate

CVODE interface def count(self, offset): self.offset = offset return len(self.v) y’, y are numpy arrays def fun(self, t, y, ydot): last = self.offset + len(self.n) self.pv.gather(self.hv) Our portion of y self.n = y[self.offset:last] if ydot == None: return ninf, nrate = self.ninftau(self.v) Our portion of y’ = f(y, t) ydot[self.offset:last] = (ninf − self.n)*nrate def msolve(self, dt, t, b, y): last = self.offset + len(self.n) self.pv.gather(self.hv) x,nrate = self.ninftau(self.v) b[self.offset:last] /= 1. + dt*nrate

CVODE interface def count(self, offset): self.offset = offset return len(self.v) def fun(self, t, y, ydot): last = self.offset + len(self.n) self.pv.gather(self.hv) self.n = y[self.offset:last] if ydot == None: return ninf, nrate = self.ninftau(self.v) ydot[self.offset:last] = (ninf − self.n)*nrate Solve M*x = b x−>b def msolve(self, dt, t, b, y): y is rarely used last = self.offset + len(self.n) self.pv.gather(self.hv) x,nrate = self.ninftau(self.v) M is implicitly 1 + dt/ntau(v) b[self.offset:last] /= 1. + dt*nrate