LIF SNN

tau*dv/dt = -v V(t) = V(0)e^-t/tau input spike .. V -> V + w if V > Vthresh the vire a spike and V -> Vreset for each time period t: Update V from value at time t to value at time t+dt Process any incoming spikes Check if V crossed the threshold If so: Emit a spike Reset V V (t _ dt) = V(t) e^-dt/tau --- Code --- # Function that runs the simulation # tau: time constant (in ms) # t0, t1, t2: time of three input spikes # w: input synapse weight # threshold: threshold value to produce a spike # reset: reset value after a spike def LIF(tau=10, t0=20, t1=40, t2=60, w=0.1, threshold=1.0, reset=0.0): # Spike times, keep sorted because it's more efficient to pop the last value off the list times = [t0, t1, t2] times.sort(reverse=True) # set some default parameters duration = 100 # total time in ms dt = 0.1 # timestep in ms alpha = np.exp(-dt/tau) # this is the factor by which V decays each time step V_rec = [] # list to record membrane potentials V = 0.0 # initial membrane potential T = np.arange(np.round(duration/dt))*dt # array of times spikes = [] # list to store spike times # run the simulation for t in T: V_rec.append(V) # record V = V * alpha # integrate equations if times and t>times[-1]: # (last e)if there has been an input spike V += w times.pop() # remove that spike from list V_rec.append(V) # record V before the reset so we can see the spike if V>threshold: # if there should be an output spike V = reset spikes.append(t)

:opt no-lint

import qualified Data.ByteString as B
import Text.Printf
:l GPlot2DH
import GPlot2DH

Constants

tau    = 6e-2      -- time constant
uRest  = -70e-3    -- base voltage
r      = 1e3       -- leak resistor ???
dt     = 0.001     -- delta t
vTh    = -54e-3    -- Firing threshold
spikeV = 10e-3     -- Firing amplitude
uReset = -80e-3    -- Firing reset voltage (sets refractory period)

Main Code

{-
   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) = (curTime+dt,uReset) : solveDE (curTime+dt) endTime i uReset sts
  | (otherwise) =
    let
      (u',sts') = if not (null sts) -- check for empty spike table
                    then
                      if curTime > head sts -- are we going to spike?
                        then
                           (u + spikeV, tail sts) -- advance spike table
                        else   
                          (u, sts)  -- spike not used 
                     else
                        (u, sts)    -- spike not used
                        
      -- Do the math and recurse
      du = ((1.0/tau) * ((-1) * (u' - uRest)) + (i * r))*dt
      newU = u' + du
      newT = curTime + dt
          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.

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$

spikeTable  =  []
results = solveDE 0 0.2 1e-3 (-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.

spikeTable  =  [0.02, 0.055, 0.07]
results = solveDE 0 0.2 9e-4 (-70e-3) spikeTable
(ts,vs) = unzip results
graphResults ts vs vTh