On Certain Mathematical Problems in Two-Photon Laser Scanning Microscopy

In an earlier paper, we derived the distribution of the number of photons detected in two-photon laser scanning microscopy when the counter has a dead period. We assumed a Poisson number of emissions, exponential waiting times, and an infinite time horizon, and used an equivalent inhomogeneous Poisson process formulation. We then used that result to improve image quality as measured by the signal-to-noise ratio. Here, we extend that study in two directions. First, we treat the finite-horizon case to assess the accuracy of the simpler infinite-horizon approximation. Second, we use a direct approach by conditioning on the Poisson count for the infinite-horizon case to derive several polynomial identities.


Introduction
Two-photon laser scanning microscopy (TPLSM) is a valuable tool for imaging living tissues. It involves the periodic excitation of molecules by hundredfemtosecond laser pulses, followed by the emission of fluorescent light before the arrival of the next pulse [5]. The count of the photons is used to estimate the intensity of a pixel in the resulting image. To convert the arrival of a photon into a pulse of current, photomultiplier tubes are commonly used [4]. This conversion leads to a dead period in the photon detector, during which it is not able to record another emitted photon. The censored photons lead in turn to an undercount, resulting in images with smaller signal-to-noise ratios and poorer image quality. Figures 1 and 2 illustrate how the dead period of TPLSM affects the photon counting process when six photons are emitted [12]. If the detector did not have a dead period the count would have a Poisson distribution, and it would be easy to estimate the intensity. However, after each photon arrival, the detector of the machine has a dead period: in Figure 2 the photons arriving at t 2 , t 4 ,  The dead periods of the photon counting systems affect the measured counting distribution. There is an extensive literature on dead periods in counting processes. Here, we describe the features relevant to TPLSM, and provide references. Dead periods are classified as either paralyzable or nonparalyzable [6,9,11]. If the dead period only occurs after registration of a pulse, we have a Type I or non-paralyzable counter; and if the dead period occurs after each pulse (whether it is recorded or not), we have a Type II or paralyzable counter. The technology of TPLSM leads to at Type I counter [7]. The existence of dead periods in photon detectors makes the photon counting process challenging. The photon arrival process N (t) counts the number of arrivals during the time interval [0, t), and most of the prevoius studies assume that N (t) is a homogeneous Poisson process with rate (photon arrivals per unit time) [1,10,14,15]. Much of the literature on dead periods uses tools from renewal theory of renewal theory to deal primarily with moments (such as means and variances) rather than probabilities. The emphasis of much of this literature is on asymptotic approximations. In contrast, our work below and in [13] calculates the exact probabilities using assumptions that are well supported by the physical background; moments and other summaries can easily be derived from these probabilities.
In an earlier paper [13], we used standard model assumptions [7]: Poisson with mean α for the number of photons emitted, exponential emission waiting times (fluorescence lifetimes), and a fixed (standardized) dead period δ upon registration of a photon. We assumed an infinite time horizon for the waiting times. Using an equivalent inhomogeneous Poisson process formulation, we derived the exact distribution of all observed counts, rather than grouped counts as in [5]. We then used maximum likelihood to estimate both α and δ from the observed counts, and showed how those estimates yield improved image quality, as measured by the signal-to-noise ratio.
We now describe the details of our model. Assume that the number of excited molecules N is a Poisson random variable with mean α. When N = n ≥ 1, assume that the photon emission waiting times W 1 , . . . , W n are independent exponential random variables with time constant τ nanoseconds (ns). An equivalent formulation is based on the following fact [3]: suppose that N (t) is an inhomogeneous Poisson process (IHPP) with intensity function λ(t) for t ≥ 0; given that N (t) = n, the n ordered occurrence times have the same joint distribution as the order statistics of n independent random variables with common density Thus, in our context, for any horizon t ≤ ∞ with N (t) = n, the ordered photon emission waiting times W (n:1) < . . . < W (n:n) have the same distribution as the occurrences of an IHPP with λ(t) = αe −t . When the arrival of a photon is recorded, the detector has a dead period of length ∆ ns. For our calculations below we use a clock in units of the time constant, so the standardized dead period is the dimensionless quantity δ = ∆/τ and the waiting times have mean 1. In typical applications, δ is around 1 and α is at most 3.0, so our numerical illustrations are for those values.
In [13] we used this inhomogeneous Poisson process approach to obtain the exact distribution of the number of photons detected, D, for an infinite time horizon. Writing P T to denote probability for horizon up to time T , for k = 1, 2, . . ., and with the convention that b k=a A k is 1 whenever b < a. In this paper, we extend our earlier results in two directions. In Section 2, we address the problem of computing the probabilities for the finite-horizon case. Because the complete distribution is not analytically tractable, we use numerical methods to give guidelines on how long the horizon must be in order for the infinite-horizon distribution to be a good approximation. In Section 3, we use a second approach to the infinite-horizon case to prove certain polynomial identities.

Distribution of D for a finite horizon
Using the same approach as in [13], we study P T (D = d) for T < ∞. The main reason for studying the finite horizon case is that, of course, actual experiments have finite length. In [5], the horizon was T = 13.2 ns, and the infinite-horizon approximation (1.1) was quite accurate there. Our calculations below will provide guidelines about how large T should be in order for that approximation to be accurate. We will see below that the finite horizon calculations are considerably more complicated than the ones leading to (1.1); in fact, they appear to be analytically intractable for d ≥ 2. One reason for the difficulty is the possibility of a dead period starting very close to the boundary T . We present numerical calculations for d ≤ 2 to support our recommendations.
The event [D = 0] happens if and only if no photon arrives at the detector before T . The probability of this event is Next, there are two possible configurations for a single photon recorded at time t: its dead period ends either before (t ≤ T − δ) or after (T − δ < t < T ) the recording period [0, T ]. The two cases are depicted in Figure 3; if it ends before T , then there must be no emissions between the end of the dead period and T . In the first case, we compute the probabilities of the events in the intervals thus: the probability of no emissions in [0, t) is an emission in [t, t + dt) is αe −t dt, and no emissions after the dead period up to T is Combining these terms, the probability of the first case, w 1 , is In the second case, the photon is recorded between T − δ and T , so that its dead period goes beyond T . As for the first case, the probability of no Combining these terms, the probability of the second case, w 2 , is Summing (2.1) and (2.2), we get the probability of recording one photon: Once again, as T → ∞, P T (D = 1) → P ∞ (D = 1). We expect that the contribution of (2.1) to be much larger than that of (2.2). To assess their relative magnitudes, we set δ = 1 and present our numerical results in Table 1. This table shows that for small rates α, (2.2) is negligible for T ≥ 4; and for larger α, it is negligible for T ≥ 2. We now assess how large T must be in order to use the infinite observation time probabilities as approximations for the finite horizon case. We compare the probabilities P T (D = 0), P T (D = 1), and an approximation to P (D = 2) with  Tables 2, 3, and 4, respectively; once again with δ = 1, and also with τ = 1.
In Table 2, we see that as T increases, the infinite horizon approximation always underestimates the finite horizon probability of [D = 0] because the longer we wait the more likely we are to see a photon. Moreover, as α gets larger P T (D = 0) converges faster to P ∞ (D = 0) as T → ∞. Next, from Table  3, we conclude that when α ≤ 1 and observation period is small, the infinite horizon approximation overestimates the probability of [D = 1]. For α > 1, it underestimates P T (D = 1), especially when T is small. The finite-horizon calculations for [D = 2] follow the same steps as the ones for the cases above. They are quite a bit more involved, leading to several double integrals some of which cannot be evaluated in elementary terms. Ignoring edge effects near T , the approximation of P T (D = 2) is Thus, we have Note that as T → ∞, P T (D = 2) → P ∞ (D = 2).
In Table 4, we compare this approximation of P T (D = 2) with the corresponding infinite-horizon values. In this table, we only consider T > 2 ns because it is impossible to see 2 photons when T ≤ 2 ns. From Table 4, we see that the infinite-horizon approximation is an overestimate (underestimate) for small (large) α; it is good approximation when T ≥ 6 ns. In summary, these tables show that it is reasonable to assume an infinite observation period when T ≥ 7 or 8 ns, depending on the accuracy desired. As noted before, Driscoll, et al. [5], used the simpler approximation because they used T = 13.2 ns.

Polynomial identities arising in photon counting problems
In this section, we use a second approach to derive the distribution of D for an infinite horizon, which leads to three sets of polynomial identities. In the course of using the first approach described above in [13] we needed the first set: Similarly, when using the approach described above, we need the second set: Unlike the first two sets of identities, the third set of identities do not appear in the steps of the approach describe in Section 2. They appear in the derivation of the distribution of D by using our second approach, which is described in the following section. For integer k ≥ 1, let B k (x) = 1 + x + · · · + x k , and let B k = B k (e δ ). For d = 3, 4, . . ., the third set of polynomial identities is Bi , for j = 0, 0, for j = 1, . . . , d − 2 and we have the convention that b k=a A k is 1. It turns out that these identities I 3 reduce to I 1 for j = 0, 1 and to I 2 for j = 2, . . . , d − 2. Now, we present our second approach to derive the distribution of D. For few cases, we provide the details of the derivation that arises these polynomials identities given in I 1 , I 2 , and I 3 .

Assumptions
We start with a Poisson(α) number of photons actually emitted and their arrivals following exponential waiting times with parameter τ . As before, we rescale the dead period δ = ∆/τ since it is convenient to use a clock with time constant units, and assume τ = 1 without loss of generality. Suppose that the photon detector has a dead period of ∆ ns; the photons that are emitted during a dead period do not affect that dead period, so this model is a Type I counter. The rescaled dead period is δ = ∆/τ . We also assume that the observation period is infinite in Section 3.2. We will later change this assumption with finite horizon in Section 3.2. Note that, below we use the notation at the end of Section 1.

Deriving distribution of D using alternate approach
Let W (n:1) < . . . < W (n:n) be the ordered waiting times when N = n ≥ 1. These order statistics can be represented thus [2]: where the X k are independent exponential random variables with time constant 1. The first few probabilities are easy. As before, because the observation period is infinite, P (D = 0) = P (N = 0) = e −α .
Next, for N ≥ 1, write Using the lack of memory property of the exponential distribution, we have P (D = 1|N = n) = P (W (n:n) ≤ W (1:n) + δ) = (1 − e −δ ) n−1 ; and since the generating function of N is E(u N ) = e −α(1−u) , we have For larger values of d we decompose the event (D = d). For the d-tuple κ where κ = κ(d) = (k 1 , . . . , k d ) with 1 = k 1 < · · · < k d ≤ n, we observe W (k 1) < · · · < W (k d ) if and only if all occur. By representation (3.4), for i = 1, . . . , d − 1, Q κ i depends only on (X ki+1 , · · · , X ki+1 ), and Q κ d depends on X k d +1 , . . . , X n ; hence, the events {Q κ i : 1 ≤ i ≤ d} are independent. By the lack of memory property, Summing over all the d-tuples, we have With this expression, we use to derive the marginal distribution of D.
Derivation of P (D = 2). First, sum over k d in equation (3.5) to get P (D = 2): To get (3.6), we set m = k d − k d−1 − 1 and then use the binomial theorem. Next, let d = 2, and substitute k 1 = 1. Then the conditional probability for d = 2 is Notice that this expression is zero for n = 1 (see (3.2)). Thus, we have The equivalence of the coefficients of ζ 0 in this expression, and for all larger d, is a special case of I 1 : see the representation given in (3.1).
Derivation of P (D = 3). We now sum over k d , k d−1 in equation (3.5). Now we define m = k d−1 − k d−2 − 1 and use the binomial theorem again to get Thus, In this case, the coefficients c i are the polynomial identities given in I 3 : see (3.3). Thus, c 1 = 0 and c 0 = e 2δ(n+1)
Notice that I 1 representation in (3.1) appears when we calculate h 0 which is the coefficient of the ζ 0 = e −α term. Also, h i terms, the coefficients of α j e −α j! for j = 1, 2, · · · , d − 1, are the polynomial identities given in I 2 : see (3.2).
Derivation of P (D = s). First, sum over k d , k d−1 , . . . , k s−2 in equation (3.5) and use binomial theorem: (3.7) Note that we have the convention that b k=a z k is 0 whenever b < a. When obtaining equation (3.7), we apply equation (3.3) for the coefficients of n−k d−2 j for j = 0, 1, . . . , s − 1. Also, for d = s, k d−s+1 = k 1 = 1, so substitute k 1 = 1 in (3.7). This yields Next we use (3.8) and apply the polynomial identities given in (3.1) and (3.2) in each step to get the marginal distribution of D = s: see (1.1).

Conclusions
In this paper, we consider extension of our model assumption and propose alternate approach to derive the distribution of D under the old and new assumption. In our earlier work [13], the distribution of D was derived by using inhomogeneous Poisson approach under the infinite time period assumption. In this work, we first address the finite horizon assumption for this model. Next, we obtain the distribution of D using an equivalent approaches by using lengthy analysis of the order statistics formula. This model is also studied under finite and infinite time horizon. When finite time assumption made in both approaches, we only calculate explicitly P (D = 0) and P (D = 1), and sketch the higher count cases because they are too involved. This shows that we can use the infinite time approximation when the observation period is long enough. Thus, we give an approximation for P (D = 1) for the finite case, and assess its accuracy. We then compute the probabilities for various values of T : our numerical work shows that when the observation period is at least 8 ns, we can safely assume an infinite horizon instead of a finite one, leading to simpler expressions for the probabilities. Finally, we used two different approaches to determining the distribution of D in order to derive a new classes of polynomial identities.
In order for this equation to hold, the coefficient of the polynomials, α i /i! must be zero, completing the proof.
Proof. I 3 identities. The proof is just a rearrangement of indexes. First, ignore e (d−1)(n+1)δ term to get that (3.11) is proportional to Then, write I 3 identities in terms of A terms by using B n = e nδ (1 + e −δ + · · · + e −nδ ) = e nδ A n+1 /A 1 : (3.11) The right side of (3.11) with j = 0 is I 1 , with d − 1 instead of d when ignoring −e −jδ A d−1−j 1 term. Next, let k = l + 1 and ignore −e −jδ A d−1−j 1 term to get that (3.11) is proportional to (3.12) The right side of (3.12) with j = 1 is I 1 , with d − 2 instead of d. Also, the left side of (3.12) with j = 2, . . . , d − 2 is I 2 , with d − 1 instead of d. This completes the proof.