LIF SNN¶

Rev 0.3

In [11]:
import qualified Data.ByteString as B
import Text.Printf
:l GPlot2DH
import GPlot2DH

Constants

In [12]:
r         = 100e6     -- leakage resistor
c         = 300e-12   -- storage cap
tau       = r*c       -- time constant
uRest     = -70e-3    -- base voltage
dt        = 0.001     -- delta t
vTh       = -54e-3    -- Firing threshold
outSpikeV = 20e-3     -- Firing amplitude
inSpikeV  = 10e-3     -- Input firing amplitude
uReset    = -80e-3    -- Firing reset voltage (sets refractory period)

Main Code

In [13]:
-- Do the simple difference equation math and
-- return a new time and a new voltage

update :: Float -> Float -> Float -> (Float,Float)
update cTime u i =
      let du = ((1.0/tau) * ((-1) * (u - uRest)) + (i * r))*dt
      in (cTime + dt, u + du)

{-
   solveDE parameters,
       current time - advances recursively
       stop time    - time to stop sim
       current      - zero for one shot, non-zero creates pulses
       voltage      - present voltage
       spikeTable   - advances recursively
-}     

solveDE :: Float -> Float -> Float -> Float -> [Float] -> [(Float, Float)]
solveDE curTime endTime i u sts
  | curTime >= endTime = []
  | u >= vTh = 
           let (tA,uA) = update curTime outSpikeV i -- spike output
               (tB,uB) = update tA uReset i         -- reset neuron
           in  (tA,uA) : (tB,uB) : solveDE tB endTime i uB sts
                                           
  | (otherwise) =
    let
       (u',sts') = case sts of
           [] -> (u, sts) -- no more input spikes
           otherwise -> if curTime > head sts -- ready to spike?
                        then (u + inSpikeV, tail sts)
                        else (u,sts)
                        
      -- Do the math and recurse
       (newT, newU) = update curTime u' i  -- generate new time and voltage
    in (newT,newU) : solveDE newT endTime i newU sts'
  
-- Show the output
showResults :: [Float] -> [Float] -> IO ()
showResults ts vs = do
  putStrLn "Action potential simulation using the LIF model with Inputs:"
  putStrLn "\n    Time(ms)    Membrane Potential (mV)"
  mapM_ (\(t, v) -> (printf "%10.2f       %10.2f\n" (t*1000) (v*1000))) $ zip ts vs

graphResults :: [Float] -> [Float] -> Float -> IO (B.ByteString)
graphResults ts vs thresh = do
    writeFile "lifwi_solution01.dat" 
      (concatMap (\(a,b) -> show (a*1000) ++ "  " ++ show (b*1000) ++ "\n") $ zip ts vs)
    runGNUPlot "lifwi_solution01.dat" 
                  "Leaky Integrate and Fire Neuronal Model with Inputs"
                   "Time (ms)" "Activation Voltage (mV)" 
                       "lifwi_solution01.png" (thresh * 1000)

Simple Spiking
$I(t) = 0$
Note that the reset value after a spike is lower than the resting value. With no spike inputs the voltage recovers exponentially to the resting value.

In [14]:
spikeTable  =  [0.02, 0.055, 0.07]
results = solveDE 0 0.2 0.0 (-70e-3) spikeTable
(ts,vs) = unzip results
graphResults ts vs vTh

Generate a Frequency
This gets interesting!
Nature created an Ion channel opening to a frequency .. pretty cool.
Setting the current to non-zero creates a spike train proportional to the current. The example below sets,
$I(t) = 1\ ma$

In [15]:
spikeTable  =  []
results = solveDE 0 0.2 1e-8 (-70e-3) spikeTable
(ts,vs) = unzip results
graphResults ts vs vTh

Current and Spikes
Don't know if nature does this .. It gets a bit complicated.

In [16]:
spikeTable  =  [0.02, 0.055, 0.07]
results = solveDE 0 0.2 1e-8 (-70e-3) spikeTable
(ts,vs) = unzip results
graphResults ts vs vTh