Schr¨odinger’s renewal Ansatz, Eq. (2.36), reduces the computation of the first-passage-time density of the Ornstein-Uhlenbeck process to the numerical so-lution of the Abel integral equation
g(t) = Z t
0
K(t, s)
√t−sF(s) ds . (A.1)
The inhomogeneity g(t) and the kernel K(t, s) are as defined by Eqs. (2.34) and (2.35). The existence of a unique solutionF(t) is assured, see Chapter2.2.
Below, a block-by-block method for the computation of F(s) is presented, which is a corrected form of the method given in Chapter 10.3 of Linz (1985).
It is validated in Section A.2.
A.1 Algorithm
Starting from a known initial pointF0 =F(t= 0), the block-by-block method iteratively computes the solution Fj =F(tj) at equidistant points tj =jh in blocks of two points. LetF0, F1, . . . , F2m be known, so that F2m+1 and F2m+2
are to be computed as next block. Then, Eq. (A.1) can be written as g2m+` =
Z t2m+`
0
K(t2m+`, s)F(s)
√t2m+`−s ds=I2m,`+
Z t2m+`
t2m
K(t2m+`, s)F(s)
√t2m+`−s ds , (A.2) with `= 1, 2. The first term on the right-hand side contains only the known portion of F(t) for s < t2m. Dividing the range of integration into m intervals of length 2h and interpolating K(t, s)F(s) quadratically within each interval,
122 Numerical Evaluation of First-Passage-Time Densities
one obtains after some algebra I2m,`= This leads to the integration weights1
wm,`,k =
The integral over the unknown portion of the function in Eq. (A.2), i.e.
over [t2m, t2m+`], is evaluated in the same manner. For ` = 1, the value of F(s) at the interval midpoint t2m+1
2 is obtained by quadratic interpolation F2m+1
2 = 38F2m+ 34F2m+1− 18F2m+2 .
Inserting into Eq. (A.2), one obtains two linear equations for the desired
func-1In the corresponding equation (8.35) inLinz (1985), p. 137, the different cases are not sufficiently separated. This is corrected here.
A.1 Algorithm 123
tion values F2m+1 and F2m+2
g2m+1−I2m,1 = α1(h2)K2m+1,2mF2m
Solving for the unknowns yields the iteration rule (m≥1) F2m+1 = 1
The precision of the algorithm is limited by that of the coefficients given in Eq. (A.8). While each of the two summands in ˜αn, ˜βn, and ˜γn are ∼O(n5/2),
124 Numerical Evaluation of First-Passage-Time Densities
h 1.0·100 2.0·10−1 1.0·10−1 1.0·10−2 1.0·10−3 4.0·10−4
∆r 3.8·10−2 5.0·10−4 7.3·10−5 8.2·10−8 8.2·10−11 5.2·10−12
tCPU[s] 4.3·10−1 4.2·101 2.5·102
Table A.1: Error ∆rof the block-by-block method for constant input vs. analytical solution of Eq. (A.10) for different step sizes h. The bottom line gives the required CPU time in seconds, with empty fields indicating tCPU <0.1s. The noise amplitude was σ= 0.1 and 0≤t≤20. All computations were performed on a COMPAQ XP1000 workstation.
their differences are only∼O(1/√
n). Using standard floating-point arithmetic with 16-digit precision, the coefficients will have only five significant digits.
Therefore, the coefficients have been tabulated to 20 significant digits using the MAPLE computer algebra package (Plesser 1999).
The block-by-block algorithm requires 4m+8 kernel evaluations to compute points F2m+1 and F2m+2 of the solution. The number of kernel evaluations are required to compute the solution at M points thus grows ∼ M2. Due to the complicated structure of the weights and the recursive nature of the algorithm, no detailed error analysis is known for this method (Linz 1985, Ch. 10.4).
This algorithm has been implemented for general inhomogeneity g(t) and kernel K(t, s) as programabel. The code has been designed to be as clear as possible. For the Ornstein-Uhlenbeck process, optimized versions have been developed for both periodic (p fptd) and aperiodic (ap fptd) stimuli. These implementations enforce non-negative solutions F to improve numerical sta-bility. Source code for both programs is available in electronic form (Plesser 1999).
A.2 Validation
Constant input
The first-passage-time density of the Ornstein-Uhlenbeck process for input of the form
I(t) = 1 + 2cet (A.9)
is given by
ρ(t) = 2 e2t
√πσ2(e2t−1)32 exp
"
−(1 +c−ce2t)2 σ2(e2t−1)
#
. (A.10)
Here, σ is the input noise level, the reset voltage is v(t = 0) = vR = 0 and the absorbing threshold is at v = 1; see Chapter 3.2. For c = 0, this FPTD
A.2 Validation 125
h 2.0·10−1 1.0·10−1 4.0·10−2 2.0·10−2 1.0·10−2 4.0·10−3 (a) 2.3·10−4 4.4·10−5 3.3·10−6 4.2·10−7 5.2·10−8 3.3·10−9 (b) 5.3·10−5 7.0·10−6 4.6·10−7 5.7·10−8 7.0·10−9 4.4·10−10 (c) 8.3·10−4 1.7·10−4 1.4·10−5 1.9·10−6 2.4·10−7 1.5·10−8 (d) 5.5·10−4 9.6·10−5 7.3·10−6 9.4·10−7 1.2·10−7 7.4·10−9 (e) 5.7·10−2 3.4·10−2 1.3·10−2 2.7·10−3 4.5·10−4 3.2·10−5 (f) 1.3·10−2 6.4·10−3 7.0·10−4 8.5·10−5 1.4·10−5 9.7·10−7
Table A.2: Error ∆r of the FPT density ρj from the block-by-block method for periodic input at different step sizes h. Input was slow (a), (b), fast (c), (d), and suprathreshold (e), (f) with weak noise for the first, strong noise for the second of each pair. The error given is relative to the numerical solution obtained with h= 2·10−3. The FPT density was computed for 0≤t≤100 in all cases. Parameters: (a)µ= 0.9,q= 0.1, Ω = 0.05/π, σ= 0.001, (b)µ= 0.9,q= 0.1, Ω = 0.05/π,σ= 0.01, (c)µ= 0.9,q= 0.1, Ω =π,σ= 0.05, (d)µ= 0.9,q= 0.1, Ω =π,σ= 0.1, (e)µ= 1.05,q= 0.5, Ω = 0.1/π,σ= 0.1, (f)µ= 1.05, q= 0.5, Ω = 0.1/π,σ= 1.
can be evaluated directly by p train, whence this case is used as reference to validate the program.
Table A.1 gives the root mean square error
∆r = v u u t
N
X
j=1
(ρj −ρ(tj))2 (A.11)
of the numerical solution ρj obtained from p train against the analytical so-lution ρ(t) of Eq. (A.10) for a wide range of step sizesh. The error decreases nearly with the third power of the step size (∆r ∼ h2.92, correlation coef-ficient r > 0.999), while the required CPU time grows inversely quadratic (tCPU ∼h−2). The implementation is thus validated for constant input.2 Periodic input
No first-passage time densities have been published for the Ornstein-Uhlenbeck process with time-dependent input. Thus, one may only test for the conver-gence of the algorithm as the step size is refined, and check for consistency with simulated spikes trains. The former is done here, while the latter is postponed to Appendix B.
The test stimulus used is I(t) = µ+qcos Ωt. Results are given in Ta-ble A.2 for six different input conditions, comprising slow and fast, sub- and suprathreshold currents, both for weak and strong noise. Across all conditions
2Linz (1985) mentions an error scaling∼h32 in Chapter 10.4, which is obviously wrong, compare Table 10.3 of that book.
126 Numerical Evaluation of First-Passage-Time Densities
the error scales roughly as ∆r ∼ h2.9 with correlation coefficient r > 0.999, indicating stability. As was to be expected, the error is larger for weak noise and/or fast stimuli than for strong noise and/or slow stimuli. But even in the
“worst case” (e) in Table A.2, i.e. a suprathreshold stimulus combined with weak noise, does the error decay nicely for h < 0.2. The algorithm is thus well suited to the numerical computation of first-passage-time densities of the modulated Ornstein-Uhlenbeck process.
When choosing the stepsize, two points should be kept in mind. First of all, the FPT density will be modulated on the timescale of the input periodT, requiring h T. This timescale is set by the cut-off frequency for aperiodic input. Secondly, the FPT density has a sharp peak for suprathreshold input in combination with a small noise amplitude or for very strong noise. The stepsize may have to be very small to resolve this peak properly and to avoid numerical instabilities.