Simulate a model neuron

Question

An Izhikevich neuron is a simple yet quite effective model of a biological neuron, designed for use in a discrete time-stepping simulation. In this golfing challenge, you will be implementing this model.

Parameters

This model involves only 7 variables organized into 2 differential equations, compared to the dozens of parameters of a physiologically accurate model.

• v and u are the two state variables of the neuron. Here, v is the "fast" variable representing the cell potential over time, and u is the "slow" variable representing certain membrane properties. The v variable is the most important one, as this is the output of the simulation.
• a, b, c, and d are fixed constants that describe the properties of the neuron. Different types of neurons have different constants, depending on the desired behavior. Notably, c is the reset potential, which is the membrane potential the cell returns to after spiking.
• I represents the input current to the neuron. In network simulations, this will change over time, but for our purposes we will treat I as a fixed constant.

The Model

This model has very simple pseudocode. First, we take the constant values of abcd and use them to initialize v and u:

v = c
u = b * c

Next, we loop through the simulation code as many times as desired. Each iteration represents 1 millisecond of time.

for 1..t:
  if v >= 30:    # reset after a spike
    v = c
    u = u + d
  v += 0.04*v^2 + 5*v + 140 - u + I
  u += a * (b*v - u)
  print v

Certain real-world implementations include additional steps for numerical accuracy, but we aren't including those here.

Input

As input, your program/function should take the values of a, b, c, d, I, and t (the number of time steps to simulate). Once set, none of these parameters will be changing during our simple simulation. The order of input does not matter: you can specify the order in which your program takes these parameters.

Output

Output will be a list of numbers representing the cell's membrane potential (given by variable v) over the course of the simulation. The list can be in any appropriate format.

You have the choice of whether to include the 0th value of the simulation (the initial configuration before any time has passed) in your output. For example, for an input of 0.02 0.2 -50 2 10 6 (for a b c d I t), an output of either

-50
-40
-16.04
73.876224
-42.667044096
-25.8262335380956
29.0355029192068

or

-40
-16.04
73.876224
-42.667044096
-25.8262335380956
29.0355029192068

is acceptable.

Your values do not have to be exactly the same as those above, depending on how your language handles floats.

Reference Implementation

Here is a TIO implementation I wrote in Perl to demonstrate the model. The parameters are that of a "chattering" neuron from the paper linked above, and this serves as a demonstration of how this model is able to recreate some of the more complex properties of neurons, such as alternating between states of high and low activity. If you look at the output, you can see where the neuron immediately spikes several times, but then waits a while before spiking several more times (despite the cell input voltage I being constant the whole time).

Answer 1

R, 110 99 bytes

Anonymous function that takes 6 arguments. Nothing fancy, just a straightforward port of the reference implementation. The updating of u, v, and the printing of v have all been combined into a single line, thanks to the fact that R's print returns the value that is being printed, so you can use it in assignment. Many thanks to Giuseppe for saving 11 bytes!

pryr::f({v=c;u=b*c;for(i in 1:t){if(v>=30){v=c;u=u+d}
u=a*b*(v=print((.04*v+6)*v+140+I-u))-a*u+u}})

Try it online!

Answer 2

Clean, 150 145 140 138 bytes

import StdEnv
$a b c d i t=map snd(iterate(\(u,v)#(w,n)=if(30.0<v)(c,u+d)(v,u)
#y=0.04*w*w+6.0*w+140.0-n+i
=(a*b*y-a*n+n,y))(b*c,c))%(0,t)

Try it online!

Defines the function $ :: Real Real Real Real Real Int -> [Real], implementing the algorithm as described in the OP, starting from the 0th term.

Answer 3

Python 2, 100 bytes

a,b,c,d,I,t=input();v=c;u=b*c
exec"if v>=30:v=c;u+=d\nv=v*v/25+6*v+140-u+I;u+=a*(b*v-u);print v\n"*t

Try it online!

Saved 2 bytes thanks to user71546.

Answer 4

JavaScript (Node.js), 107 ... 103 101 bytes

Contributed by @apsillers

(a,b,c,d,I,t)=>[...Array(t)].map(_=>(v<30||(v=c,u+=d),v=v*(v/25+6)+140-u+I,u+=a*(b*v-u),v),u=b*(v=c))

Try it online!

Original Approach: 105 103 bytes. -1 byte Thanks Arnauld, and -2 bytes Thanks @Kamoroso94.

(a,b,c,d,I,t)=>{for(u=b*(v=c);t--;){v<30||(v=c,u+=d);v=v*(v/25+6)+140-u+I;u+=a*(b*v-u);console.log(v)}}

Try it online!

Or if popping alerts is OK, then 101 ... 99 97 bytes (-1 byte Thanks Arnauld, -2 bytes Thanks @Kamoroso94):

(a,b,c,d,I,t)=>{for(u=b*(v=c);t--;){v<30||(v=c,u+=d);v=v*(v/25+6)+140-u+I;u+=a*(b*v-u);alert(v)}}

var u, v;
var f = 
(a,b,c,d,I,t)=>{for(u=b*(v=c);t--;){v<30||(v=c,u+=d);v=v*(v/25+6)+140-u+I;u+=a*(b*v-u);alert(v)}}

function run() {
 f(...["a", "b", "c", "d", "I", "t"].map(x => document.getElementById(x).value * 1));
}

a = <input id="a" value="0.02"><br>
b = <input id="b" value="0.2"><br>
c = <input id="c" value="-50"><br>
d = <input id="d" value="2"><br>
I = <input id="I" value="10"><br>
t = <input id="t" value="6"><br>
<input type="button" value="Run" onclick="run()">

Answer 5

Ruby, 94 bytes

->a,b,c,d,i,t{v=c
u=b*c
t.times{v>=30?(v=c;u+=d):0
v+=0.04*v**2+5*v+140-u+i
u+=a*(b*v-u)
p v}}

Try it online!

Another straightforward port of the reference implementation, a lambda accepting 6 arguments.

Answer 6

Haskell, 112 111 bytes

(a#b)c d i t|let r(v,u)|v>=30=r(c,u+d)|p<-0.04*v^2+6*v+140-u+i=(p,u+a*(b*p-u))=fst<$>take t(iterate r$r(c,b*c))

Try it online!

Does not output the zero case. Assumes that c is never >=30 since that wouldn't make sense.

Never thought I'd have to use a where clause in a code golf but there are just too many variables.

EDIT: Thanks @Lynn for taking off a byte! I forgot that you can put let statements in guards. Sure kills the readability though

Answer 7

Element, 81 bytes

_a;_b;_3:b~*u;_d;_I;_'[3:\.04*5+*140+u~-+I~++4:\
.`30<!b~*u~-+a~*u~+[d~+]u;[#2:]]

Try it online!, Esolangs page

Explanation:

_a;_b;_3:b~*u;_d;_I;_'[ ... ]

This part of the program takes input. It stores constants a, b, d, and I into variables. The input for c is never stored in a variable, but rather remains on the main stack throughout execution. Three copies are made: one on top for initializing u, one in the middle to serve as the initial v, and one on the bottom to serve as the constant c. The input for t is thrown immediately onto the control stack to serve as the basis of the FOR loop (the [...]) surrounding the rest of the program.

3:\.04*5+*140+u~-+I~++4:

This part of the program takes the current value of v and calculates the new value, and then four copies of the new v value are made.

\
.`

The first copy of v has a newline appended and is printed.

30<!

The second copy of v is used to test if the neuron has spiked. The result of this test is put on the control stack for later use.

b~*u~-+a~*u~+

This part computes the "delta u", meaning the amount to add to u.

[d~+]

This IF block adds d to the above sum if the neuron is spiking. This combines what would normally be two assignments into a single assignment.

u;

This stores the updated value of u.

[#2:]

This IF block is a continuation of the above IF block. If the neuron is spiking, delete the current value of v (which is now on the top of the main stack), and replace it with a duplicate of c (which has been on the bottom of the main stack this entire time).

And that's basically all there is to it. One minor note is that this thing leaks memory: it takes an extra "# in there to delete the top of the control stack (the evaluated IF condition) after each loop iteration.

Although I wouldn't call Element the most graceful golfing language, this challenge does allow me to showcase an interesting feature: due to the division between the main stack and control stack, I can take an IF statement and split the condition and body into multiple parts, interlaced with unconditional code.

Answer 8

MATLAB, 111 bytes

function z(a,b,c,d,I,t)
v=c;u=b*c;for i=1:t if v>=30 v=c;u=u+d;end
v=.04*v^2+6*v+140-u+I
u=u+a*(b*v-u);
end
end

Pretty simple implementation, can probably be further improved.