Solution to an ODE (Leaky Integrate and Fire Neuron Model)

In summary, the conversation discusses the search for a general analytical solution to a common ODE in neuroscience. It is believed that such a solution does not exist, but the reasons for its difficulty are not fully understood. The equations for the model are outlined and some unsuccessful attempts at solving it are presented. Finally, a hypothesis involving the Lambert W function is proposed and tested numerically, showing promising results.
  • #1
madness
815
70
I'm looking for a general analytical solution to a particular ODE that comes up in neuroscience a lot. My feeling is that such a solution can't be obtained, otherwise someone would have presented it by now, but I don't have a good understanding of why it is so hard to solve.

The equations as they are usually given:

Let ##\dot{V} = \frac{dV}{dt}## and ##I(t)## be some function (if it helps you to find a solution, you could assume certain conditions on ##I(t)##, but I'm trying to find the solution for the general case). The model is defined by the following two equations: 1) ## \dot{V} = - V + I(t) ## while ##V(t) < 1##, and 2) if ##V(t) = 1## then ## \lim_{\epsilon\rightarrow 0} V(t + |\epsilon|) = 0 ##

Some further notation:

Let ##\{ t^{(1)}, t^{(2)}, ...\} = \{ t: V(t) = 1\} ##. The key problem is to identify this set.

Some attempts at a solution:

The equation is piecewise linear. We can solve for any interval ##\left[t_0,t\right] \subset \left(t^{(i-1)}, t^{(i)}\right)## as ## V(t) = V(t_0) e^{-t} + \int_{t_0}^t e^{-t'} I(t-t') dt' ##. So the solution in general for ##t^{(i-1)} < t < t^{(i)} ## is ## V(t) = \int_{t^{(i-1)}}^{t} e^{-t'} I(t-t') dt' ##. So the times ##t^{(i)}## could be identified by solving the equation ## 1 = \int_{t^{(i-1)}}^{t^{(i)}} e^{-t'} I(t^{(i)}-t') dt' ## - I've tried Taylor expanding ##I(t)## which allows the integral to be solved as a series but then I got stuck.

We can also compress the two usual rules of the model into a single equation ## \dot{V} = - V + I(t) - \sum_i \delta(t-t^{(i)}) ##, then we can identify the general solution ##V(t) = e^{-t}V(0) - \sum_i e^{-(t-t^{(i)})} \Theta(t-t^{(i)}) + \int_0^t e^{-t'}I(t-t')dt'## where ##\Theta(x)## is the Heaviside step function, but of course we don't know the ##t^{(i)}##s. Equivalently, we can write as ## \dot{V} = - V + I(t) - |\dot{V}|\delta(V-1) ## (I believe this is correct, I used a rule for the composition of a function with a delta function). Then because ##|\dot{V}| \ge 0 ## at the point ##V=1## we have ##|\dot{V}| = \dot{V}## and so ##(1+\delta(V-1)) \dot{V} = -V+I(t)##. This equation is nice and compact, but looks quite tricky to solve.

I'd be curious to see if anyone can provide some further insights!
 
Last edited:
  • Like
Likes Delta2
Physics news on Phys.org
  • #2
I'm pretty confused by the limit you wrote down. Is this equivalent to solving for V(t), and then deciding to replace some of the places where V(t)=0 with V(t)=1? I feel like probably not, since that wouldn't be an interesting condition to apply.
 
  • #3
Office_Shredder said:
I'm pretty confused by the limit you wrote down. Is this equivalent to solving for V(t), and then deciding to replace some of the places where V(t)=0 with V(t)=1? I feel like probably not, since that wouldn't be an interesting condition to apply.

If V(t) reaches 1 from below, it is reset to V(0) at the next instant. You solve for V(t) while V(t)<1, and if you find that V(t)=1 at any point you immediately reset V to 0.
 
  • #4
I see. I did not think correctly about what that absolute value sign was doing.
 
  • #5
I wanted to post my latest attempt at a solution. I'm not sure if it is any way helpful but I found it interesting.

Recall the model equations: ##\dot{V} = -V+I(t)## while ##V(t)<1## , and where ##V(t)=1## then ##\lim_{\epsilon\rightarrow 0} V(t+|\epsilon|) = 0##. Define ##t^{(i)}## as the times where ##V(t)=1## and is reset to 0. Define ##\Phi(t) = \int e^tI(t)dt##.

For ##t^{(i-1)} < t < t^{(i)} ##, ##V(t) = e^{-t}\int_{t^{(i-1)}}^t e^{t'} I(t') dt'##, so we have ##1 = e^{-t^{(i)}}\int_{t^{(i-1)}}^{t^{(i)}} e^{t'} I(t') dt' = e^{−t^{(i)}} (\Phi(t^{(i)})−\Phi(t^{(i−1)}))\implies e^{t^{(i)}}−\Phi(t^{(i)})=\Phi(t^{(i−1)})##. So if we define ##\Psi(t) = (e^t-\Phi(t))^{-1}\circ\Phi(t)##, then ##t^{(i)}=\Psi(t^{(i−1)})=\Psi(\Psi...\Psi(t^0)...)## .
 
Last edited:
  • #6
madness said:
I wanted to post my latest attempt at a solution. I'm not sure if it is any way helpful but I found it interesting.

Recall the model equations: ##\dot{V} = -V+I(t)## while ##V(t)<1## , and where ##V(t)=1## then ##\lim_{\epsilon\rightarrow 0} V(t+|\epsilon|) = 0##. Define ##t^{(i)}## as the times where ##V(t)=1## and is reset to 0. Define ##\Phi(t) = \int e^tI(t)dt##.

For ##t^{(i-1)} < t < t^{(i)} ##, ##V(t) = e^{-t}\int_{t^{(i-1)}}^t e^{t'} I(t') dt'##, so we have ##1 = e^{-t^{(i)}}\int_{t^{(i-1)}}^{t^{(i)}} e^{t'} I(t') dt' = e^{−t^{(i)}} (\Phi(t^{(i)})−\Phi(t^{(i−1)}))\implies e^{t^{(i)}}−\Phi(t^{(i)})=\Phi(t^{(i−1)})##. So if we define ##\Psi(t) = (e^t-\Phi(t))^{-1}\circ\Phi(t)##, then ##t^{(i)}=\Psi(t^{(i−1)})=\Psi(\Psi...\Psi(t^0)...)## .

I suggest you test your hypothesis with a real example. For example, let ##I(t)=t##. Then if :
##\Psi(t) = (e^t-\Phi(t))^{-1}\circ\Phi(t)##
I get:
##\Phi(t)=\int te^t=te^t-e^t##
Now let ##e^t-(te^t-e^t)=u##. Surprisingly, we can solve for t in terms of u via the lambert W function, albeit multi-valued. We obtain:
##
t=2+W\left(\frac{-u}{e^2}\right)
##
so that
##\Psi(t) = (e^t-\Phi(t))^{-1}\circ\Phi(t)=2+W\left(-\frac{te^t-e^t}{e^2}\right)##

Now, solve the DE numerically and compute the spikes. Do they then agree with ##\Psi(\Psi...\Psi(t^0)...)##?
 
  • #7
Cool, nice work! See attached code and plot :). I took the first brach of the W function and plotted the real part.

Code:
clear all

%% solve by simulation

t(1) = 0;
V(1) = 0;

dt = 0.001;

Nt = 100000;

for i=1:Nt
    
    V(i+1) = V(i) + (-V(i) + (i*dt))*dt;
    
    if V(i+1) >= 1 
        V(i+1) = 0;
        
        t = [t,i*dt];
        
    end
end
  

%% solve by Lambert W

clear tl

tl(1) = 0;

for i=1:(length(t)-1)

tl(i+1) = 2+lambertw(0,-(tl(i) * exp(tl(i)) - exp(tl(i)))/exp(1)^2);

end

%% Plot

hold on
scatter(1:length(t), t)
scatter(1:length(real(tl)), real(tl))
legend('Numerical (dt = 0.001)', 'Iterated Lambda W (real part)')
xlabel('Spike Number')
ylabel('Spike Time')
[\CODE]
[ATTACH type="full"]270734[/ATTACH]
 

Attachments

  • LIF_Analytical.png
    LIF_Analytical.png
    5.9 KB · Views: 235
  • #8
Here's my check of your formula in Mathematica. I had to remove the negative sign in the W-function in order for the predicted spike values to agree with the DE. Not sure why:

Mathematica:
myF[t_] := 2 + ProductLog[(t Exp[t] - Exp[t])/E^2];
tPrev = 0;
tNext = myF[0];
theSpikeData = Table[
   tNext = myF[tPrev] // N;
   mySol =
    First@NDSolve[{v'[t] + v[t] == t, v[tPrev] == 0},
      v, {t, tPrev, tNext}];
   solTrace[t_] = Evaluate[Flatten[v[t] /. mySol]];
   thePlot = Plot[solTrace[t], {t, tPrev, tNext}];
   result = {solTrace[t], tPrev <= t < tNext};
   tPrev = tNext // N;
   result, {i, 0, 10}];

endT = (Last@theSpikeData)[[2, 5]];
theSpikeTrain[t_] = Piecewise@theSpikeData;
Plot[theSpikeTrain[t], {t, 0, endT}, ExclusionsStyle -> {Red, Blue}]
leaky integrate and fire.jpg
 
  • #9
aheight said:
Here's my check of your formula in Mathematica. I had to remove the negative sign in the W-function in order for the predicted spike values to agree with the DE. Not sure why:

Mathematica:
myF[t_] := 2 + ProductLog[(t Exp[t] - Exp[t])/E^2];
tPrev = 0;
tNext = myF[0];
theSpikeData = Table[
   tNext = myF[tPrev] // N;
   mySol =
    First@NDSolve[{v'[t] + v[t] == t, v[tPrev] == 0},
      v, {t, tPrev, tNext}];
   solTrace[t_] = Evaluate[Flatten[v[t] /. mySol]];
   thePlot = Plot[solTrace[t], {t, tPrev, tNext}];
   result = {solTrace[t], tPrev <= t < tNext};
   tPrev = tNext // N;
   result, {i, 0, 10}];

endT = (Last@theSpikeData)[[2, 5]];
theSpikeTrain[t_] = Piecewise@theSpikeData;
Plot[theSpikeTrain[t], {t, 0, endT}, ExclusionsStyle -> {Red, Blue}]
View attachment 270742

You're right, I made a sign error - it should be ##\Psi(t) = (\Phi(t) - e^t)^{-1}\circ \Phi(t)##. Nice plot. You can see in my plot that the the discrepancy between the two solutiuons accumulates over time (in your code you start the DE from the analytical solution after each spike, in mine I do a spike-reset when the DE voltage exceeds 1). Presumably the Lambert W solution is more accurate, but in Matlab it actually seems to take longer to run than the numerical solution of the DE! (I suppose the time the DE takes depends on the chosen step size, so it could be faster or slower depending on chosen accuracy).
 

Related to Solution to an ODE (Leaky Integrate and Fire Neuron Model)

1. What is the Leaky Integrate and Fire Neuron Model?

The Leaky Integrate and Fire (LIF) Neuron Model is a mathematical model used to describe the behavior of a single neuron. It takes into account the leaky nature of the neuron's membrane, which allows ions to flow in and out, and the integration of incoming signals from other neurons. When the membrane potential reaches a certain threshold, the neuron fires an action potential, or electrical impulse, down its axon.

2. How is the LIF Neuron Model used to solve ODEs?

The LIF Neuron Model is a system of Ordinary Differential Equations (ODEs) that can be solved using numerical methods. The equations describe the change in membrane potential over time and can be solved to predict when a neuron will fire an action potential. By solving these equations, we can gain insight into the behavior of neurons and how they respond to different stimuli.

3. What are the key parameters in the LIF Neuron Model?

The key parameters in the LIF Neuron Model include the membrane capacitance, the membrane resistance, the resting potential, the threshold potential, and the reset potential. These parameters determine the behavior of the neuron and can be adjusted to simulate different types of neurons or to study the effects of different stimuli on the neuron's firing behavior.

4. How accurate is the LIF Neuron Model in predicting neuron behavior?

The LIF Neuron Model is a simplified model that does not take into account all the complexities of a real neuron. Therefore, it is not always accurate in predicting neuron behavior. However, it can provide valuable insights and is often used as a starting point for more complex models. The accuracy of the model also depends on the accuracy of the parameters used and the assumptions made.

5. What are the limitations of the LIF Neuron Model?

One of the main limitations of the LIF Neuron Model is that it does not take into account the spatial distribution of ion channels on the neuron's membrane. It also does not consider the effects of synaptic inputs from other neurons. Additionally, the model assumes that the neuron is a point source and does not take into account the dendritic structure of real neurons. These limitations make the model less accurate in predicting the behavior of real neurons.

Similar threads

  • Differential Equations
Replies
1
Views
804
  • Differential Equations
Replies
1
Views
830
  • Differential Equations
Replies
5
Views
702
  • Differential Equations
Replies
1
Views
716
Replies
3
Views
847
Replies
3
Views
1K
Replies
1
Views
1K
  • Differential Equations
Replies
16
Views
946
Replies
7
Views
2K
  • Differential Equations
Replies
2
Views
2K
Back
Top