For a given spike train s = (s₁, s₂, ⋯ , s_k) in the time domain [0, T], denote s₀ = 0 and s_k+1 = T. Using the notion of the inter-spike intervals (ISI), this point process sequence can be equivalently represented as a vector u=(u1,u2,⋯,uk+1)T=(s1-s0,s2-s1,⋯,sk+1-sk)T, where ^T indicates the transpose operation. This ISI vector is in a simplex Sk+1 in ℝ^k+1, where Sk+1={u=(u1,u2,⋯,uk+1)T|∑i=1k+1ui=T,ui>0,i=1,⋯,k+1}.

By applying the Isometric LogRatio (ILR) transformation [10], the ISI vector is mapped from the simplex Sk+1 to an unconstrained Euclidean space ℝ^k. Specifically, for any u=(u1,u2,⋯,uk+1)T∈Sk+1,

where g(u) is the geometric mean of u, and Ψ ∈ ℝ^k×(k+1) represents the ILR transform matrix. This matrix satisfies two conditions: 1) ΨΨT=Ik, and 2) ΨTΨ=Ik+1-1k+11k+11k+1T, where Ik∈ℝk×k is the identity matrix and 1_k+1 is a column vector of ones in ℝ^k+1. We point out that this transformation is a bijection between Sk+1 and ℝ^k, and its inverse takes the following form:

In this section, we will focus on the distribution of u^* ∈ ℝ^k from spike trains following a homogeneous Poisson process (HPP). As pointed out in Qi et al. [6], given the cardinality k, the ISI of an HPP is uniformly distributed on the simplex space Sk+1 and its density function has the form fu|k(u|k)=k!Tk. By applying the change-of-variable technique, the density function of u^* can be written as follows:

where c_k is the normalizing constant for the density.

It was shown in Zhou et al. [9] that the conditional density function in Equation (2) owns important and desirable properties: (1) the density is log-concave and uni-modal, and (2) the density has a symmetry with respect to the origin in a simplex. These properties make this density analogous to a multivariate normal distribution, which is ideal for a center-outward rank on the observed data. Our new rank-based model on the point process will be based on this conditional density function.

In this subsection, we will formally define the joint likelihood of the ILR-transformed ISI for an HPP. This new likelihood is referred to as the rank-based likelihood.

To compute the joint density function of the ILR-transformed homogeneous Poisson process, we can decompose it as a product of the conditional density function and the marginal density function. That is, for any given spike train s = (s₁, s₂, ⋯ , s_k) in the time domain [0, T], and u^* ∈ ℝ^k as the ILR transformation of the ISI of s, the joint likelihood function of u^* can be written as:

By using the one-to-one mapping between the spike train s and the ILR-transformed vector u^*, the conditional density function in Equation (2) can be re-described using s as: fu*|k(u*|k)=ckTk+1∏i=1k+1(si-si-1) (see detailed derivation in Zhou et al. [9]). In addition, the spike train s and the ILR-transformed vector u^* must have the same cardinality. Hence P(|u^*| = k) = P(|s| = k), which is the probability of observing k events in the original HPP space. It is well known that in an HPP with constant rate λ, the number of events in any interval of length T would be a Poisson random variable with mean λT, and then the probability of this HPP has k events is P(|s|=k)=e-λT(λT)kk!. Therefore, the likelihood of u^* can be expressed as:

In this paper, we call this likelihood as the rank-based likelihood of the original spike train s. Its formal definition is given as follows:

Definition 2.1. Let s = (s₁, s₂, ⋯ , s_k) be a realization of a homogeneous Poisson process on [0, T] with constant rate λ > 0. Denote s₀ = 0 and s_k+1 = T. Then the rank-based likelihood of s is defined as:

where c_k is the normalizing constant given in Equation (2).

As opposed to the traditional likelihood of a spike train in a homogeneous Poisson process where the number of events is the only factor considered, the temporal distribution of the events plays a crucial role in determining the rank-based likelihood of the spike train.

For each given cardinality k, it is easy to see that the likelihood can attain its maximum value when the ISI is a constant vector u=(1k+1T,⋯,1k+1T). In this case, all spike times are 1k+1T,2k+1T,⋯, and kk+1T, which uniformly distribute on [0, T] and properly indicate the homogeneous pattern in the data. This result shows that the newly derived likelihood function can be used to rank a homogeneous Poisson process in any given cardinality.

HPP Outlier Detection: For outlier detection in a homogeneous Poisson process, we will need to incorporate both the cardinality and the distribution of events into the identification. We can slightly modify the newly defined likelihood calculation by normalizing the conditional density by its maximum value within each respective cardinality. This modified likelihood function on a process s with k events is given by:

In this subsection, we will extend the rank-based likelihood for homogeneous Poisson process to general point process spike trains. In this paper, we propose to achieve this goal by transforming the general point process into a homogeneous Poisson process using the classical Time-Rescaling (TR) method [3]. The TR method allows any point process to be transformed into a homogeneous Poisson process as long as the conditional intensity function is known (or can be estimated). Based on this method, we propose to define the rank-based likelihood for a general point process that utilizes such transformation. This definition allows for the computation of the rank-based likelihood in a feasible way by taking advantage of the simplified structure of the homogeneous Poisson process. Our definition is formally given as follows:

Definition 2.2. Let s = (s₁, s₂, ⋯ , s_k) be a realization of a general point process on [0, T] with k time events. Denote s₀ = 0 and s_k+1 = T. If the conditional intensity function λ^* is given and the cumulative function Λ*(t)=∫0tλ*(u)du,t∈[0,T], then we can transform s to Λ*(s)=(Λ*(s0),Λ*(s1),⋯,Λ*(sk+1)) via the time-rescaling transformation. Denote u*=(u1*,u2*,⋯,uk+1*)T as the ILR transformation of the ISI u=(u1,u2,⋯,uk+1)T=(Λ*(s1)−Λ*(s0),Λ*(s2)−Λ*(s1),⋯,Λ*(sk+1)−Λ*(sk))T. Then the rank-based likelihood of s conditioned on its cardinality is defined as:

where c_k is the normalizing constant given in Equation (2).

We point out that the above definition is only for the conditional likelihood because a general form of the marginal distribution on the number of events for a general point process is not available. However, if we only focus on the center-outward rank within the same cardinality, there is no need to calculate such marginal likelihood.

One commonly used special case is the inhomogeneous Poisson process (IPP). In this case, the conditional intensity is deterministic and can be written as λ^*(t) = λ(t). Let the cumulative intensity Λ(t)=∫0tλ(u)du,t∈[0,T]. Then the number of events in [0, T] follows a Poisson distribution with mean Λ(T). The probability of this process having k events is given by P(|s|=k)=e-Λ(T)(Λ(T))kk!. Therefore, the likelihood of an inhomogeneous Poisson process s with k spikes can be expressed as:

IPP outlier detection: Similar to the HPP case, for outlier detection in an IPP, we also need to incorporate both the cardinality and the distribution of events into the identification. This is done by normalizing the conditional density with its maximum value for each cardinality. The modified likelihood function for a process s with k events is given by:

In this section, we will illustrate the usage of the proposed likelihood for ranking and outlier detection in simulated spike trains. We at first generate 5000 realizations on [0, 1] from a homogeneous Poisson process with a constant rate of 5. For illustrative purposes, the ranking results with a cardinality of 4 are shown in Figure 1. There are three rows in this panel: the first row shows the constant intensity, and the second and third rows show observations with top and bottom 5 likelihoods in Equation (5). It is apparent that the spikes in the top observations are uniformly distributed, indicating a typical homogeneous pattern, whereas the spikes in the bottom observations appear outliers from the homogeneity. We emphasize that such ranking cannot be obtained with the classical likelihood method as all observations with cardinality 4 will have equal likelihood values.

Based on Equation (6), we display the outliers (observations with 5 lowest modified likelihood values) in the first row of Figure 1B. These spike trains clearly do not exhibit a homogeneous pattern. In contrast, we display spike trains with 5 lowest classical likelihoods in the second row. These trains all have only 1 spike. These single-spike trains properly address the outlier type on the number of events, whereas they cannot capture the outliers with respect to the distribution of events.

We then generate 5,000 inhomogeneous Poisson process realizations on [0, 1] from the intensity function λ(t) = 3e^3t (exponentially increasing from 0 to 1), and the result for the cardinality being 12 is shown in Figure 1C. There are also 3 rows in this panel: the first row shows the intensity function, and the second row shows observations with the top five new likelihoods (left) in Equation (8) and classical likelihoods (right), respectively. We can see that the top spike trains using the new likelihood better represent intensity function than those using the classical likelihood. The third row in the panel shows observations with the bottom 5 new likelihoods (left) and classical likelihoods (right), respectively. In this case, both outlier groups exhibit different patterns from the intensity function.

We finally show the result across different cardinalities in the same inhomogeneous Poisson process realizations. In Figure 1D, we find the spike trains with top (upper row) and bottom (lower row) 5 modified likelihood estimations using Equation (9). We can also see that the typical spike trains (with high likelihoods) well capture the variability in intensity function, and the outliers (with low likelihoods) clearly exhibit discrepancy from the intensity.