Solution (source code)

= Solution

For a constant-rate <Poisson process>, let $T$ be the waiting time from a spike to the next spike. By <independent increments>, there is no spike in an interval of length $u$ with probability $e^{-ru}$. Hence, for $r>0$,
$$
\boxed{F_T(u)=\mathbb P(T\leq u)=
\begin{cases}0,&u<0,\\1-e^{-ru},&u\geq0.\end{cases}}
$$
Its density is $re^{-ru}$, so the interspike intervals are independent <exponential distributions> with mean $1/r$. If $U$ is uniform on $(0,1)$, <inverse transform sampling> gives $T=-\log U/r$. Repeatedly add independent draws to the current spike time to simulate a constant-rate <spike train>.

For a prescribed time-dependent nonnegative rate $r(t)$, use an <Inhomogeneous Poisson process>. Conditional on the present event at $t_i$, its next waiting time has survival probability
$$
\mathbb P(T_i>u\mid t_i)=\exp\left[-\int_{t_i}^{t_i+u}r(s)\,ds\right].
$$
Thus its conditional CDF is $1$ minus that expression for $u\geq0$. Define the cumulative intensity $\Lambda(t)=\int_0^t r(s)\,ds$. For each new independent uniform draw, set $E_i=-\log U_i$ and solve
$$
\boxed{\Lambda(t_{i+1})-\Lambda(t_i)=E_i.}
$$
This transforms independent unit-rate exponential intervals in integrated-intensity time into physical spike times. Use a generalized inverse when the rate vanishes on an interval: no spikes are placed in a flat part of $\Lambda$. On a finite simulation window, stop when the next intensity target exceeds the available $\Lambda$. If the total future intensity is finite, a target beyond it means that no further spike occurs.

Numerically, integrate the known rate and invert the accumulated integral, or integrate exactly across piecewise-constant rate bins. \b[Drawing an exponential interval using only $r(t_i)$ is not exact when the rate changes during the interval.] The Poisson model also omits the <neuronal refractory period> and other spike-history effects.