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

[madness]

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!

 

[Office_Shredder]

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.

 

[madness]

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.

 

[Office_Shredder]

I see. I did not think correctly about what that absolute value sign was doing.

 

[madness]

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)...)$ .

 

[aheight]

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)...)$?

 

[madness]

Cool, nice work! See attached code and plot :). I took the first brach of the W function and plotted the real part.

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]

 

[aheight]

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:

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}]

 

[madness]

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:

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).