text: 
  from numpy import *
  from pylab import plot, show, find, rand, randn
  
  # Created by Eugene M. Izhikevich, February 25, 2003
  # Excitatory neurons   Inhibitory neurons
  Ne = 800
  Ni = 200
  
  re = rand(Ne)
  ri = rand(Ni)
  
  a = r_[0.02 * ones(Ne), 0.02 + 0.08 * ri]
  b = r_[0.2 * ones(Ne), 0.25 - 0.05 * ri]
  c = r_[-65 + 15 * re**2, -65 * ones(Ni)]
  d = r_[8 - 6 * re**2, 2 * ones(Ni)]
  S = c_[0.5 * rand(Ne + Ni, Ne), -rand(Ne + Ni, Ni)]
  
  v = -65 * ones(Ne + Ni)# Initial values of v
  u = b*v                   # Initial values of u
  firings = zeros((0,2))
  
  for t in xrange(1000):# simulation of 1000 ms
      I = r_[5 * randn(Ne), 2 * randn(Ni)] # thalamic input
      fired = find(v >= 30)# indices of spikes
      if any(fired):
          firings = vstack((firings, c_[t + 0 * fired, fired]))
          v[fired] = c[fired]
          u[fired] = u[fired] + d[fired]
          I = I + S[:,fired].sum(1)
      v = v + 0.5 * (0.04 * v**2 + 5 * v + 140 - u + I)
      v = v + 0.5 * (0.04 * v**2 + 5 * v + 140 - u + I)
      u = u + a*(b*v - u)
  
  plot(firings[:,0], firings[:,1], '.')
  show()
  