Efficient Inference of Flexible Interaction in Spiking-neuron Networks
Authors: Feng Zhou
Yixuan Zhang
Jun Zhu
Venue: ICLR 2021 Poster (2021)
Keywords: neural spike train
nonlinear Hawkes process
auxiliary latent variable
conjugacy

Abstract: Hawkes process provides an effective statistical framework for analyzing the time-dependent interaction of neuronal spiking activities. Although utilized in many real applications, the classic Hawkes process is incapable of modelling inhibitory interactions among neurons. Instead, the nonlinear Hawkes process allows for a more flexible influence pattern with excitatory or inhibitory interactions. In this paper, three sets of auxiliary latent variables (Polya-Gamma variables, latent marked Poisson processes and sparsity variables) are augmented to make functional connection weights in a Gaussian form, which allows for a simple iterative algorithm with analytical updates. As a result, an efficient expectation-maximization (EM) algorithm is derived to obtain the maximum a posteriori (MAP) estimate. We demonstrate the accuracy and efficiency performance of our algorithm on synthetic and real data. For real neural recordings, we show our algorithm can estimate the temporal dynamics of interaction and reveal the interpretable functional connectivity underlying neural spike trains. 

Published as a conference paper at ICLR 2021
EFFICIENT INFERENCE OF FLEXIBLE INTERACTION IN SPIKING-NEURON NETWORKS
Feng Zhou', Yixuan Zhang', Jun Zhu
'Dept. of Comp. Sci. & Tech., BNRist Center, THU-Bosch Joint ML Center, Tsinghua University Data Science Institute, University of Technology Sydney
(2h0ufeng6288, desz (tsz)))tsiinghtne edu. en, yixuan zhangeuts.
ABSTRACT
Hawkes process providss an effective statistical framework for analyzing the time- dependent interaction of neuronal spiking activities. Although utilized in many real applications, the classic Hawkes process is incapable of modelling inhibitory interactions among neurons. Instead, the nonlinear Hawkes process allows for a more flexible influence pattern with excitatory or inhibitory interactions. In this paper, three sets of auxiliary latent variables (Pólya-Gamma variables, la- tent marked Poisson processes and sparsity variables) are augmented to make functional connection weights in Gaussian form, which allows for simple it- erative algorithm with analytical updates. As result, efficient expectation- maximization (EM) algorithm derived to obtain the maximum a posteriori (MAP) estimate. We demonstrate the accuracy and efficiency performance of our algorithm on synthetic and real data. For real neural recordings, we show algorithm estimate the temporal dynamics of interaction and reveal the interpretable functional connectivity underlying neural spike trains.
INTRODUCTION
One of the most important tracks in neuroscience is to examine the neuronal activity in the cere- bral cortex under varying experimental conditions. Recordings of neuronal activity are represented through series of action potentials spike trains. The transmitted information and functional connection between neurons are considered to be primarily represented by spike trains (Kass et al., 2014; Kass Ventura, 2001; Brown et al., 2004; 2002). A spike train is a sequence of recorded times at which a neuron fires an action potential and each spike may be considered to be times- tamp. Spikes occur irregularly both within and across multiple trials, so it reasonable to consider spike train as point process with the instantaneous firing rate being the intensity function of point processes (Perkel et al., 1967; Paninski, 2004; Eden et al., 2004). An example of spike trains for multiple neurons is shown in Fig. 2a in the real data experiment.
Despite many existing applications, the classic point process models, e.g., Poisson processes, ne- glect the time-dependent interaction within one neuron and between multiple neurons, so fail to capture the complex temporal dynamics type of point processes which is able model the self-exciting interaction between past and future events. Existing applications cover a wide range of domains including seismology (Ogata, 1998; 1999), criminology (Mohler et al., 2011; Lewis et al., 2012), financial engineering (Bacry et al., 2015; Filimonov & Sornette, 2015) and epidemics (Saichev & Sornette, 2011; Rizoiu et al., 2018). Unfortunately, due to the linearly additive intensity, the vanilla Hawkes process can only represent the purely excitatory interaction because a negative firing rate may exist with inhibitory interac- tion. This makes the vanilla version inappropriate in the neuroscience domain where the influence between neurons is mixture of excitation and inhibition (Maffei et al., 2004; Mongillo et al., 2018).
neural population. In contrast, Hawkes process is one
In order reconcile Hawkes process with inhibition, various nonlinear Hawkes process variants are proposed to allow for both excitatory and inhibitory interactions. The core point of nonlinear Hawkes process is a nonlinearity which maps the convolution of the spike train with a causal influ- ential kernel to nonnegative conditional intensity, such as rectifier (Reynaud- RReynaud-Bouret et al., 2013),
"orresponding author.

---

Published as a conference paper at ICLR 2021
exponential (Gerhard et al., 2017) and sigmoid (Linderman, 2016; Apostolopoulou et al., 2019). The sigmoid mapping function has the advantage that the Pólya-Gammn Gamma augmentation scheme can be utilized to convert the likelihood into Gaussian form, which makes the inference tractable. In Lin- derman (2016), a discrete-time model is proposed to convert the likelihood from a Poisson process to Poisson distribution. Then Pólya-Gamma random variables are augmented on discrete obser- vations to propose a Gibbs sampler. This method is further extended to a continuous-time regime in Apostolopoulou al. (2019) by augmenting thinned points and Bolya-Gamma random variables to propose a Gibbs sampler. However, the influence function is limited to be purely exciting or inhibitive exponential decay. Besides, due to the nonconjugacy of the excitation parameter of expo- nential decay influence function, a Metropolis- Gibbs sampler making the Markov chain Monte Carlo (MCMC) algorithm further inefficient.
sampling step has to be embedded into the
To address the parametric and inefficient problems in aforementioned existing works, we develop flexible sigmoid nonlinear multivariate Hawkes processes (SNMHP) model the continuous-time regime, (1) which can represent the flexible
temporal dynamics among the neural population, (2) with the efficient conjugate inference. An EM inference algorithm is proposed to fit neural spike trains. Inspired by Donner & Opper (2017; 2018), three auxiliary latent variable sets: P6lya-Gamma variables, latent marked Poisson processes and sparsity variables are augmented make functional connection weights in a Gaussian form. As result, the EM algorithm has analytical updates with drastically improved efficiency. As shown in experiments, it is even more efficient than the maximum likelihood estimation (MLE) for the parametric Hawkes process in high dimensional cases.
2
OUR MODEL
Neurons communicate with each other by action potentials (spikes) and chemical neurotransmitters. A spike causes the pre-synaptic neuron to release chemical neurotransmitter that induces impulse responses, either exciting or inhibiting the post-synaptic neuron from firing its own spikes. The ad- dition of excitatory and inhibitory influence neuron determines whether a spike will occur. At the same time, the impulse response characterizes the temporal dynamics of the exciting or inhibit- ing influence which can be complex and flexible (Purves al., 2014; Squire al., 2012: Bassett & Sporns, 2017). Arguably, the flexible nonlinear multivariate Hawkes processes are suitable choice for representing the temporal dynamics of mutually excitatory or inhibitory interactions and functional connectivity of neuron networks.
2.1 MULTIVARIATE HAWKES PROCESSES
The vanilla multivariate Hawkes processes (Hawkes, 1971) are sequences of timestamps D (L W J. E [0, 7) where 12 is the timestamp of n-th event on i-th dimension with N; being the number points on i-th dimension, M the number of dimensions, the observation window. The i-th dimensional conditional intensity, the probability of an event occurring on i-th dimension in [t, t + dt) given all dimensional history before f, is designed in a linear superposition form:
M
Ax(1) - +
Pig(t - (H)),
j=1d
where μ4 > 0 the baseline rate of i-th dimension and Φ: 2 is the causal influence function (impulse response) from j-th dimension to i-th dimension which is normally a parameterized func- tion, e.g., exponential decay. The summation explains the self- and mutual-excitation phenomenon, i.e., the occurrence of previous events increases the intensity of events in the future. Unfortunately, one blemish is the vanilla multivariate Hawkes processes allow only nonnegative (excitatory) in- fluence functions because negative (inhibitory) influence functions may yield negative intensity which is meaningless. To reconcile the vanilla version with inhibitory effect and flexible influence function, we propose the SNMHP.
2

---

Published as a conference paper at ICLR 2021
2.2 SIGMOID NONLINEAR MULTIVARIATE HAWKES PROCESSES
Similar to the classic nonlinear multivariate Hawkes processes (Brémaud & Massoulié, 1996), the i-th dimensional conditional intensity of SNMHP is defined
(x(() (h)())())) ha(() p4
Σ EE (t- 1 j=
(2)
where b is the base activation of neuron ht(() is real-valued activation and (()) is the logistic (sigmoid) function which maps the activation into positive real value in 1) with A; being upper-bound to scale (0, A²). The sigmoid function is chosen because as seen later, the Pólya- Gamma augmentation scheme can be utilized to make the inference tractable. After incorporating see the influence functions, (15(-), can be positive neg-
the nonlinearity, is straightforward ative. If dig (-) is negative, the superposition of (-)(-) will lead to a negative activation (t) that renders the intensity to O; instead, the intensity tends to X; with positive (1)(-))
To achieve a flexible impulse response, the influence function is assumed to be a weighted sum of basis functions
() Σ wigsdn())
where { are predefined basis functions and Wiyb is the weight capturing the influence from j-th dimension to i-th dimension by b-th basis function with positive indicating excitation and neg- ative indicating inhibition. The basis functions are nonnegative functions capturing the temporal dynamics of the interaction. Although basis functions can be in any form, in order for the weights to represent functional connection strength, basis functions are chosen to be probability densities with compact support that means they have bounded support [0, Tal and the integral is one. As a result, the i-th dimensional activation
M
M
h1(1) = μ4 +
Σ
wi Wiyrcholl 12,) = 14: +
Wiyh Σ
E66l1 (6(1-12) (1)) 12))
j=1 P. Pt b=1
j=1
(4)
M
B
+
po = = -t)),,
where (() is the convolution of j-th dimensional observation with b-th basis function and can be precomputed; w, [e,, Will...., "aMa)" and (t) [1, Dn (t)........ (t) both (MB + 1) x 1 vectors. similar model is used in Linderman (2016) where a binary variable is included to characterize the sparsity of functional connection. As shown later, the sparsity in model is guaranteed by utilizing Laplace prior on weight instead.
In this paper, the basis functions are scaled (shifted) Beta densities, but alternatives such Gaussian or Gamma also can be used. The reason we choose Beta distribution is the inference of weights will be subject edge effects with infinite support densities when close the endpoints of [o, To]. The weighted sum of Beta densities is natural choice. With appropriate mixing, it can be used to approximate functions on bounded intervals arbitrarily well (Kottas, 2006).
INFERENCE
The likelihood of a point process model provided in Daley & Vere-Jones (2003). Correspondingly, the probability density (likelihood) of SNMHP on the i-th dimension function of parameters in continuous time is
N: II 2=1
Ao (ooh)(())) exp (-f." Nco (co((aa())d
PDD W., Xx)
(5)
is worth noting that (t) depends on w, and observations on all dimensions. Our goal is to infer the parameters i.e., weights and intensity upper-bounds, from observations, e.g., neural spike trains,
3

---

Published as a conference paper at ICLR 2021
over time interval [0, T) The functional connectivity in cortical circuits is demonstrated to be sparse in neuroscience (Thomson & Bannister, 2003; Sjöström et al., 2001). To include sparsity, factorizing Laplace prior is applied on the weights which characterize the functional connection. With the likelihood Eq. and Laplace prior pr.(w;) = I1 ta corresponds to Ll penalized log-likelihoon The i-th dimensional MAP estimate can expressed as
2a exp (- ( ), the log-posterior
X = argmax ((og((D)wr X,) + log PL lo
(6)
where w and X MAP estimates. The dependency of the log-posterior on parameters is com- plex because the sigmoid function exists in the log-likelihood term and the absolute value function exists in the log-prior term. As result, we have closed-form solutions for the MAP estimates. Numerical optimization methods can be applied, but unfortunately, the efficiency is low due to the high dimensionality of parameters which is (MB 2) M. To circumvent this issue, three sets of auxiliary latent variables: Polya-Gamma variables, latent marked Poisson processes and sparsity variables are augmented make the weights appear in a Gaussian form in the posterior. As a result, an efficient EM algorithm with analytical updates is derived to obtain the MAP estimate.
3.1 AUGMENTATION OF PÓLYA-GAMIA VARIABLES
Following Polson et al. (2013), the binomial likelihoods parametrized by log odds can be represented as mixtures of Gaussians w.r.t. representation of the sigmoid function
Gamma distribution. Therefore, we can define a Gaussian
((2) = 1 (1( (p,
where fw,,2) =2/2-2w/2--0 and pro ( 0) is the Pólya-Gamma Gamma distribution with E R+ Substituting Eq. 7 into the likelihood Eq. 5, the products of observations (f,) are transformed into a Gaussian form.
3.2 AUGMENTATION OF MARKED POISSON PROCESSES
Inspired by Donner & Opper (2018). a latent marked Poisson process is augmented to linearize the exponential integral term in the likelihood. Applying the property of sigmoid function ((2) 1 - -- and Eq. the exponential integral term is transformed to
exp (-1"" (( exp (- (1-e((lwu- (wp)(()11,)ddudl) (8)
The right hand side characteristic functional of marked Poisson process. According to the Campbell's theorem (Kingman, 2005) (App. I), the exponential integral term can be rewritten as
exp (- (oo = Eps;
fw. (w.t)E..
II
where IT, = ((wt K= denotes a realization of a marked Poisson process and px, is the proba- bility measure of the marked Poisson process II, with intensity \x(t, = ( 0). The events {tt KKK follow a Poisson process with rate X; and the latent Pólya-Gamma variable wf denotes the independent mark at each location t.. We can see that, after substituting Eq. 9 into the likelihood Eq. 5, the exponential integral term is also transformed into Gaussian form.
3.3 AUGMENTATION OF SPARSITY VARIABLES
The augmentation of two auxiliary latent variables above makes the augmented likelihood become Gaussian form w.r.t the weights. However, the absolute value in the exponent of the Laplace prior hampers the Gaussian forn of weights the posterior. To circumvent this issue, we augment the third set of auxiliary latent variables: sparsity variables. It has been proved that Laplace
4

---

Published as a conference paper at ICLR 2021
distribution can be represented infinite mixture of Gaussians (Donner & Opper, 2017; Pontil et al., 2000)
1.00
Bigb 2TTa
Bigh 2a Wish p(3rss)dBej.
(((Wish) = 2a exp
exp
(10)
where p((35b)) (Bizb/2)= exp (-1/(28.)6)) It is straightforward to see the weights are trans-
formed into a Gaussian form the prior after the augmentation of latent sparsity variables B.
3.4 AUGMENTED LIKELIHOOD AND PRIOR
After the augmentation of three sets of latent variables, we obtain the augmented joint likelihood and prior (derivation in App. II)
p(D, 1T WilW;, X.) II
(o (ed (ee Ps, (TT.).)
II (w,t)Ell;
(11a)
MB+1
Bagb 2
[-Bygb Bijb 2a² (-aa²
2
p(w,, Ba)
II
V
exp
Bijb
exp
2/zzbb
(11b)
j.b
where wE is the vector of w each t Bi is (MB + 1) vector of [3,00, Ba1. B:MBl Xa wn) (w,, 0). The motivation of augmenting auxiliary latent variables should now be clear: the augmented likelihood and prior contain the weights in Gaussian form, which corre- sponds quadratic expression for the log-posterior (L.1 penalized log-likelihood).
3.5 EM ALGORITHM
The original MAP estimate has been represented by Eq. 6. With the support of auxiliary latent variables, we propose an analytical EM algorithm obtain the MAP estimate instead of perform- ing numerical optimization. In the standard EM algorithm framework, the lower-bound (surrogate function) of the log- peg--ooterior can represented as
Q(w, A X = En,..,, [log p( D, II; wi/W; (o,/:);;)))). A.)] + Eg, [log (w,, 3:)],
(12)
2-1 "o:"",," .xx- and w 1- ,x, )
with expectation over posterior distributions p(II, w indicating parameters from last iteration.
E step: Based on joint distributions in Eq. 11, the posterior of latent variables can be derived. The detailed derivation is provided in App. III. The posterior distributions of V6lya-Gamma variables wi and sparsity variables B, and the posterior intensity of marked Poisson process II, are
N.
p(ws/w?-1)
II
PrG (w) ((,)),
(13a)
n=1
w/t, 1- .x,
) X; (-
(t))pro(w|1, M (f)),
(13b)
MB+1
p(B.lw
II PiG(Bijo j.b
w 1),
(13c)
where Wt( is the posterior intensity of I1 Pio the inverse Gaussian distribution. It worth noting that (t) depends on w The first order moments, E[wg]] 1/(2ng- tanh(h (t, (tt))2) and E[Biyb] = a a/w will be used in the M step.
M step: Substituting Eq. 13 into Eq. 12, we obtain the lower-bound Q(w,, X // - The updated parameters be obtained by maximizing the lower-bound The detailed derivation is provided in App. III. Due to the augmentation auxiliary latent variables, the update of parameters has closed-form solution
x = (N; + K2) KT)
(14a)
w E, [" (t)d),
(14b)
5

---

Published as a conference paper at ICLR 2021
where K; -. A,(t,wlw;-1,, ) /dwdt, E, = =;= [Io' (((t) (t)dt + diag (o (o-21((.1)
with diag(-) indicating the diagonal matrix of a vector, A.(1) = EW ( + Jo wt, (t) Σ 5(t 5(t
( -(, 1 T ( with o(-) being the Dirac delta function. It is worth noting that numerical quadrature methods, e.g., Gaussian quadrature, need to be applied to intractable integrals above.
3.6 COMPLEXITY AND HYPERPARAMETERS
The complexity of our proposed EM algorithm ONNNN, B (NN 1)3)) where N is the number of observations on all dimensions, Nr, is the the average num- ber of observations on the support of To on all dimensions and L is the number of iterations. The first term is due to the convolution nature of Hawkes process, the second and third term to the matrix multiplication and inversion in EM iterations. For one application, the number of dimensions M and basis functions B are fixed and much less than N. Therefore, the complex- ity can be simplified as ONN + L)).
Algorithm 1: EM inference for SNMHP
M (MB+ +
Result: {A:(1) Xioo {
Predefine basis functions (ob(-))=ii Initialize the hyperparameter a and (Ax, W.. w;. II, B.) for Iteration do for Dimension do
Update the posterior of w; by Eq. 13a; Update the posterior intensity of 11, by Eq. 13b;
Update the posterior of Be by Eq. 13c;
Update the intensity upper-hound X by Eq. 14a;
The hyperparameter in Laplace prior that en- codes the sparsity of weights and parameters of basis functions can be chosen by cross valida- tion or maximizing the lower-bound Q using numerical methods. For the number of basis functions: in essence, a large number leads to more flexible functional space while a small number results in a faster inference. experiments, we gradually increase it until no more significant improvement. Similarly, the number of quadra- ture nodes and EM iterations is also gradually increased until a suitable value. The pseudocode is provided in Alg. 1.
Update the weights w, by Eq. 14b.
end
Update the hyperparameter a.
end
4 EXPERIMENTS
We validate the EM algorithm for SNMHP in analyzing both synthetic and real-world spike data collected from the cat primary visual cortex. For comparison, the following most relevant baselines considered: (1) parametric linear multivariate Hawkes processes that are vanilla multivariate Hawkes processes with exponential decay influence functions, for which the inference is performed by MLE (Ozaki, 1979); (2) nonparametric linear multivariate Hawkes processes with flexible influ- ence functions, for which the inference is by majorization minimization Euler-L Lagrange (MMEL) (Zhou al., 2013); (3) parametric nonlinear multivariate Hawkes processes with exponential de- cay influence functions, for which the inference is by MCMC based on augmentation and Poisson thinning (MCMC-Aug) (Apostolopoulou et al., 2019), The implementation of our model is publicly ctt htt.aa/......oooooooooon////////
available at
4.1 SYNTHETIC DATA
We analyze spike trains obtained from the synthetic network model shown in Fig. la. The synthetic neural network contains four groups of two neurons each. In each group, the 2 neurons are self- exciting and (shifted) Beta distributions as basis functions with support [0. To = 6) in Fig. 1b. For the ground
while groups are independent of each other. We assume 4 scaled
truth, it is assumed that Φ1 Φ3 Øss Ôग7 Φ1. Φ22 Φ4ه Φns das 2.,
$12 031- Os6 ØTs - 1 Ø2. Φ21 ΦA3 des Øs7 _os with positive indicating excitation and negative indicating inhibition. With base activation
{4 0 and upper-bounds
(7)) 5, we use the thinning algorithm (Ogata, 1998) generate two sets of synthetic spike

---

Published as a conference paper at ICLR 2021
WATTOA TIE
YWMAL YIM!
(a)
(b)
(c)
(d)
Funcioral Comedrity Goand trun Fandoral Corned ety Caimit
Ligineissed
VLC MICMC--u
R
qakcatuuue tockss
Auraged polio cieeaioo
(c)
(h)
Figure 1: The synthetic network model and experimental results. (a): The synthetic neural popu- lation contains 4 independent groups. In each group, the interdependencies between 2 neurons are self-exciting and mutual-inhibitive with red arrows indicating excitation and blue arrows indicating inhibition. (b): Four scaled (shifted) Beta densities basis functions the support of [0, 6). (c): The intensities and spike times of neurons in the synthetic data. (d): The estimated influence functions of 1-st and 2-nd neurons where the estimated 611, d12, 621, 622 are close to the ground truth, the other ground truth 113....8 and Φ233...8 are not labeled since they all zero (GT=Ground Truth). (e): The heat map of functional connectivity among neural population with ground truth (left) and estimation (right). (f): The training and test log- likeeiiolooo curve w.r.t. EM iterations. (g): The trade-off' between accuracy and efficiency w.r.t. # of quadrature nodes and basis functions for synthetic data. (h): The running time of 2D data for EM algorithm and alternatives w.r.t. the average observation number on each dimension (the precomputation of (t) is included).
data on the time window [0, T = 1000] with one being the training dataset in Fig, Ic and the other one test dataset in App. IV. Each dataset contains sequences and each sequence consists of 3340 events on average. We aim to identify the functional connectivity of the neural population and the temporal dynamies of influence functions from statistically dependent spike trains. More experimental details, e.g., hyperparameters, are given in the App. IV.
The temporal dynamics of interactions among the neural population is shown Fig. 1d where we plot the estimated influence functions of 1-st and 2-nd neurons (other neurons are shown in the App. IV). The estimated Φ11 and du2 exhibit the self-exciting relation with d12 and du1 character- izing the coutual-inhibiie interactions. All estimated influence functions are in a flexible form and close the ground truth. Besides, shown in Fig. Ie, the estimated functional connectivity covers the ground-truth structure successfully. The functional connectivity is defined as J (ovs (d) meaning there no connection only if neither excitation nor inhibition exists.
The training and test log-likelihood (LogL) curves w.r.t. EM iterations are shown in Fig. If where our EM al- gorithm converges fast with only 50 iterations needed to obtain a plateau. The trade-off between accuracy (LogL) and efficiency (running time) w.r.t. the num- ber of quadrature nodes and basis functions is shown in Fig. 1g where we can see the accuracy is not sensitive to the number quadrature nodes over 100 and the optimal number of basis functions is 4. A larger number does not significantly improve the accuracy but leads to longer running time. Moreover, we compare the running time of our method with alternatives in Fig. 1h where the number of dimen- sions M is fixed to 2, basis functions B to 4, quadrature nodes to 200 and iterations of all methods to 200. We can observe that our EM algorithm is the most efficient, even superior to MLE for the classic parametric case, which verifies its efficiency. Also, compare our model's fitting and pre-
Table 1: Training/test Logl. (x10²) of different models for synthetic data.
MLE
MMEL.
MCMC-Aug
FM
Training LogL
2.051
1.993
2.199
Test Lagl.
1.866
1.843
2.278
2,373

---

Published as a conference paper at ICLR 2021
funden de
1 Kredin
Fundlisa Connecilty Eséeraton
[x100mal
Nawent
-1
Harce landennss irance Arddon 6
Schwiion
-11
trsiiing inplioal hood oT bigiialnsoee
www. RU
(x100mg
1x100mi
(a)
(b)
(c)
(d)
Figure 2: The real data experimental results. (a): The training spike trains extracted from real data (test spike trains in App. IV). (b): The estimated influence functions between 8-th and 9-th neurons.
(c): The heat map estimated functional connectivity among 25 neurons. (d): The training and test Logl. curves w.r.t. EM iterations.
diction ability with baseline models for 1-st and 2-nd neurons. Training and test LogL are shown in Tab. where our SNMHP with EM inference is the champion due to its superior generalized expressive ability.
4.2 REAL DATA
In this section, we analyze our model performance on a real multi-neuron spike train dataset. We aim to draw some conclusions about the functional connectivity of cortical circuits and make inferences of the temporal dynamics of influence.
Spike Train Data (Blanche, 2005; Apostolopoulou et al., 2019) Several multi-chammn silicon elec- trode arrays are designed to record simultaneously spontaneous neural activity of multiple isolated single units in anesthetized paralyzed cat primary visual cortex areas 17. The spike train dataset contains spike times of 25 simultaneously recorded neurons.
We extract the spike times in the time window [0, 300] (time unit: 100ms, the same applies to the following) as the training data (Fig. 2a) and [300, 600] as the test data (App. IV). Both datasets contain approximate 7000 timestamps. All hyperparameters are fine tuned to obtain the maximum test LogL: the scaled (shifted) Beta distribution Beta(a = 50, 8 50, shift 55) with support [0, Ta = 10] is designed as the basis function; the number of quadrature nodes is set to 1000 and EM iterations to 100. More experimental details, e.g., hyperparameters, are given in the App. IV.
Preliminary Setup
Results 25 x 25 influence functions among the neuron population estimated in the application. An example the influence functions between 8-th and 9-th neurons are plotted in Fig. 2b where SNMHP model successfully captures the exciting or inhibitive interaction between neurons. Be- sides, the estimated functional connectivity is shown in Fig. 2c where we can see the functional connection structure among neural population is sparse. Unfortunately, because the ground-truth functional connectivity of cortical circuits is unknown, the estimated functional connectivity cannot be compared with the ground truth but here we resort to the test LogL verify whether the estima- tion is good. The training and test LogL curves are shown in Fig. 2d where they both reach a close plateau indicating the estimation is appropriate without overfitting or underfitting.
A significant advantage of our EM algorithm is the efficiency. The 25-dimensional observation in the real data is a challenge for the inference. For the running time, our EM algorithm costs 3 min- utes, the MCMC-Aug costs I hour and 45 minutes with the same number of iterations while MLE and MMEL cannot finish in 2 days due to the curse of dimensionality. Moreover, the fitting and prediction ability is compared in Tab. 2. The superior performance of SNMHP w.r.t. training and test Logl.
Table 2:
Training/test LogL (x103) and run- ning times of different models for real data.
MLE
MMEL
MCMC-Aug
EM
Training Logl.
-15.228
-5.519
Test Logl.
-6.135
-5.862
Running Time
2 days
1b 45m
3m
8

---

Published as a conference paper at ICLR 2021
demonstrates model capture the complex mixture of exciting and inhibitive interactions among neural population which leads to better goodness-of of-fit.
5 DISCUSSION AND CONCLUSION
Although propose point-estimation method (EM algorithm) in this work, straightforward extension to Gibbs sampler is already at hand. Based on the augmented likelihood and prior, we can obtain the conditional densities of latent variables and parameters in closed form, which constitutes a Gibbs sampler with better efficiency than MCMC-Aug since the mome-conuunig Metropolis- Hasting sampling in MCMC-Aug not needed. However, the proposed Gibbs sampler is less efficient than the proposed EM algorithm because the latent Poisson processes have to be sampled by thinning algorithm Gibbs sampler which time consuming. For the model Apostolopoulou et al. (2019). a tighter intensity upper-bound is used to reduce the number of thinned points to accelerate the sampler. Instead, our EM algorithm does not encounter this problem we compute the expectation rather than sampling. Moreover, Apostolopoulou et al. (2019) can only use one basis function, which limits influence functions be purely exciting or inhibitive exponential decay. Instead, model can utilize multiple basis functions to characterize an influence function that is mixture of excitation and inhibition.
In this paper, we develop SNMHP model in the continuous-time regime which can characterize mixture temporal dependencics among the neural population. Three auxiliary latent variables augmented to make the corresponding EM algorithm in closed form to improve efficiency. The synthetic and real data experimental results confirm that our model's accuracy and efficiency are superior to the state of the arts. From the application perspective, although proposed in the neuroscience domain, it can be applied to other applications where the inhibition vital factor, e.g., in the coronavirus (COVID-19) spread, the inhibitive effect may represent the medical treatment or cure, or the forced isolation by government. From the inference perspective, point-estimation method; other efficient distribution-estinee estimation methods can
model
our EM algorithm be developed, e.g., the Gibbs sampler mentioned above or the mean-fich. finnff variational inference.
ACKNOWLEDGMENTS
The authors would like to thank the anonymous reviewers for insightful comments which greatly im- proved the paper. This work was supported by NSFC Projects (Nos. 62061 136001, 61620106010), Beijing NSF Project (No. JQ19016). Beijing Academy of Artificial Intelligence (BAAI), Tsinghua- Huawei Joint Research Program, a grant from Tsinghua Institute for Guo Qiang, Tiangong Institute for Intelligent Computing, and the NVIDIA NVAIL Program with GPU/DGX Acceleration. F. Zhou was partially funded by China Postdoctoral Science Foundation.
REFERENCES
Ifigeneia Apostolopoulou, Scott Linderman, Kyle Miller, and Artur Dubrawski. Mutually regressive point processes. In Advances in Neural Information Processing Systems, pp. 5116-5127, 2019.
Emmanuel Bacry, lacopo Mastromatteo, and Jean-François Muzy. Hawkes processes in finance. Market Microstructure and Liquidity, 1(01):1550055, 2015.
Danielle S Bassett and Olaf Sporns Network neuroscience. Nature neuroscience, 20(3):353, 2017.
Tim Blanche. The neural data was recorded by Tim Blanche in the laboratory of Nicholas Swindale, University of British Columbia, and downloaded from the NSF-funded CRCNS Data Sharing website., 2005.
Pierre Brémaud and Laurent Massoulié. Stability of nonlinear hawkes processes. The Annals of Probability, pp. 1563-1588, 1996.
Emery N Brown, Riccardo Barbieri, Valérie Ventura. Robert Kass, and Loren M Frank. The time- rescaling theorem and its application to neural spike train data analysis. Neural computation, 14 (2):325-346 2002.

---

Published as a conference paper at ICLR 2021
Emery N Brown, Robert E Kass, and Partha P Mitra. Multiple neural spike train data analysis: state-of-the-art and future challenges. Nature neuroscience, 7(5):456-461 2004.
Daryl J Daley and David Vere-Jones. An introduction to the theory of point processes. vol. prob- ability and its applications, 2003.
Christian Donner and Manfred Opper. Inverse Ising problem in continuous time: A latent variable approach. Physical Review E, 96(6):062104, 2017.
Christian Donner and Manfred Opper. Efficient Bayesian inference of sigmoidal Gaussian Cox processes. The Journal of Machine Learning Research, 19()):2110-2433, 2018.
Uri Eden, Loren M Frank, Riccardo Barbieri, Victor Solo, and Emery N Brown. Dynamic analysis of neural encoding by point process adaptive filtering. Neural compulation, 16(5):771-99) 2004.
Vladimir Filimonov and Didier Sornette. Apparent criticality and calibration issues in the Hawkes self-excited point process model: application to fegefrreyenccy financial data. Quantitative Fi- nance, 15(8):1293-1333 2015.
Felipe Gerhard, Moritz Deger, and Wilson Truccolo. On the stability and dynamics of stochastic spiking neuron models: Nonlinear Hawkes process and point process GLMs. PLoS computational biology, 13(2), 2017.
Alan G Hawkes. Spectra of some self-exciting and mutually exciting point processes. Biometrika, 58(1))83-90, 1971.
Robert E Kass and Valérie Ventura. A spike-train probability model. Neural computation, 13(8): 1713-1720, 2001.
Robert E Kass, Uri Eden, and Emery N Brown. Analysis of neural data, volume 491. Springer, 2014.
John Frank Charles Kingman. Poisson processes. Encyclopedia of biostatistics, 6, 2005.
Athanasios Kottas. Dirichlet process mixtures of Beta distributions, with applications density and intensity estimation. In Workshop on Learning with Nonparametric Bayesian Methods, 23rd International Conference on Machine Learning (ICML), volume 47, 2006.
Erik Lewis, George Mohler, P Jeffrey Brantingham, and Andrea L Bertozzi. Self-exciting point process models of civilian deaths in Iraq. Security Journal, 25(3):244 264, 2012.
Scott Warren Linderman. Bayesian Methods for Discovering Structure in Neural Spike Trains. PhD thesis, Harvard University, 2016.
Arianna Maffei, Sacha Nelson, and Gina Turrigiano. Selective reconfiguration of layer 4 visual cortical circuitry by visual deprivation. Nature neuroscience, 7(12):335-1135) 2004.
George O Mohler, Martin B Short, P Jeffrey Brantingham, Frederic Paik Schoenberg, and George E Tita. Self-exciting point process modeling of crime. Journal of the American Statistical Associa- fion, 106(493):100-100 2011.
Gianluigi Mongillo, Simon Rumpel, and Yonatan Loewenstein. Inhibitory connectivity defines the realm of excitatory plasticity. Nature neuroscience, 21(10):1463-1470 2018.
Yosihiko Ogata. Space-time point- proness models for earthquake occurrences. Annals of the Insti- nute of Statistical Mathematics, 50(2):379 402, 1998.
Yosihiko Ogata. Seismicity analysis through point-process modeling: A review. In Seismicity pat- terns, their statistical significance and physical meaning, pp. 471-507. Springer, 1999.
Tohru Ozaki. Maximum likelihood estimation of Hawkes' self-exciting point processes. Annals of the Institute of Statistical Mathematics, 31(1):145-155, 1979.
Liam Paninski. Maximum likelihood estimation of cascade point-process neural encoding models. Network: Computation in Neural Systems, 15(4):243-262 2004.
10

---

Published as a conference paper at ICLR 2021
Donald H Perkel, George L Gerstein, and George P Moore. Neuronal spike trains and stochastic point processes: II. Simultaneous spike trains. Biophysical journal, 7(4):419 440, 1967.
Nicholas G Polson, James G Scott, and Jesse Windle. Bayesian inference for logistic models using Pólya-Gamma latent variables. Journal of the American statistical Association, 1336 1349, 2013.
Massimiliano Pontil, Sayan Mukherjee, and Federico Girosi. On the noise model of support vector machines regression. In International Conference on Algorithmic Learning Theory, pp. 316-324. Springer, 2000.
Dale Purves, George J Augustine, David Fitzpatrick, WC Hall, AS LaMantia, JO McNamara, and White. Neuroscience, 2008. De Boeck, Sinauer, Sunderland, Mass, pp. 15-16, 2014.
Patricia Saynaud-BBurde Vincent Rivoirard, and Christine Tulcau-Malot, Inference of functional connectivity neurosciences via Hawkes processes. In 2013 IEEE Global Conference on Signal and Information Processing, pp. 317-320. IEEE, 2013.
Marian-Andrei Rizoiu, Swapnil Mishra, Quyu Kong, Mark Carman, and Lexing Xie. SIR-Hawkes: linking epidemic models and Hawkes processes to model diffusions in finite populations. In Proceedings of the 2018 World Wide Web Conference, pp. 419-428, 2018.
Al Saichev and Didier Somette. Generating functions and stability study of multivariate self-excited epidemic processes. The European Physical Journal B, 83(2):271, 2011.
Per Jesper Sjöström, Gina G Turrigiano, and Sacha B Nelson. Rate, timing, and cooperativity jointly determine cortical synaptic plasticity. Neuron, 32(6): 1499-1164 2001.
Larry Squire, Darwin Berg, Floyd E Bloom, Sascha Du Lac, Anirvan Ghosh, and Nicholas C Spitzer. Fundamental neuroscience. Academic Press, 2012.
Alex M Thomson and A Peter Bannister. Interlaminar connections in the neocortex. Cerebral cortex, 13(1):5-14, 2003.
Ke Zhou, Hongyuan Zha, and Le Song. Learning triggering kemels for multi dimensional Hawkes processes. In International Conference Machine Learning, pp. 1301-1309, 2013.

---

Published as a conference paper at ICLR 2021
APPENDIX
I CAMPBELL'S THEOREM
Let I = ((z, wn) In be a marked Poisson process on the product space 2 = x with intensity A(z, w) = A(z)p(wlzz) A(z) is the intensity for the unmarked Poisson process {znn V.V. with Wn p(w,, knn) being an independent mark drawn each Znn Furthermore, define function ((,ww) 2 x Ω + and the sum H)) = -(awe) w(z,w). If A(2, < 00, then
exp exp V/ (cchlew) wdwod
Bng [exp (EH(
for any & E (C. The above equation defines the characteristic functional of marked Poisson process. This proves Eq.9 in the main paper. The mean is
Em [111 = / l(zz, w) W) w)dwdz,
which is used when substituting Eq. 13 into Eq. 12.
II DERIVATION OF AUGMENTED LIKELIHOOD AND PRIOR
Substituting Eq.7 and into Eq.5 in the main paper, the augmented likelihood obtained
N: II dio
exp (-1² (xo(h(())d)))
D(D]w: A.)
(ho((s)())))
No
(o Xyed prolws /1, O)cck,
II a=1
Epss
II
fww,))(1))
(w.t)ETI
//
II
Paitth w/, Je Ph: (T,;
II (w,t)ell,
where wi is the vector of w and wh) = \ppro(wll1,0) It is straightforward to see the augmented likelihood is
N2
p(D, wi = II
Aeet w (e
PA, (TI,AA) (wit)eT.
II
ef(e))
n=1
which is Eq.
Similarly, the integrand in Eq. 10 is just the augmented prior in Eq. 11b.
III DERIVATION OF EM ALGORITHM
In the standard EM algorithm framework, the lower-bound of log-posterior has been provided in Eq. 12. The posterior of latent variables can be derived from the joint distribution in Eq. 11. The derivation is relatively easy for w; and B3, while II; is difficult. In the following, and mean the last and current iteration in the EM algorithm.
STEP
1. The posterior of Polya-Gamma variables w, is dependent on the activation (6-(()) {t) (nnnee which is further dependent w;" 'troough Eq.
(t)
N: PPG (wi 1,- (t,))
p(w:|w?-1 = II
12

---

Published as a conference paper at ICLR 2021
where utilize the tilted Pólya- Gamm- density (w)b,)) x e- pro(web 0) (Polson et al., 2013).
2. The posterior of sparsity variables Gg is an inverse Gaussian distribution which is dependent weights w."
MB+1
II
(
2-11 1).
5,b
w
15b
The posterior of II; is dependent on both h and X
PA, (1I I[w.tenm ew 2(0)) PX. (T ') I[(t)) ef -1- (t),
w
The Campbell's theorem can be applied to convert the denominator, the equation above can be transformed as
Ps. (TI, 1I I[ ef(w.-n;- (0))
p(T,, W
exp (1 - -A, (0)))
- )T pro(w/1, pro(w(1,O)ddul)) 0)dwvlt)
pro(w/1,0) ) exp (-1/ c/w 100x; PPo((Wl, O)dwdt
= II (c
(w,f)
The above posterior distribution is in the likelihood form of a marked Poisson process with intensity function
A xx = ff (00x (pa(w/1,0) = =x x,- (1))p((()))
M STEP
Substituting posterior distributions of latent variables into Eq. 12, obtain the lower-bound Q. The first term of Eq. 12 is
En, [logp(D, w., = -2222 w L'" (.(t))(() (t)dt w. w. w.+w! w/ // (t))
X (N++ Ad(t: log I, + "ft+/A/(t:)aaaan
where we utilize the mean rule in Campbell's theorem, C' is a constant and
N; Σ
Elwijo(t - (2) f wW,(t,
A,
N
(t
6(t - (2,) - 2 Σ
/ A)(t,))dw,
2 ==1
with 8(-) being the Dirac delta function and B[w 1/(2h=-((()) (anh( (Polson et al., 2013). The integral of intensity function has no closed-form solution but can be solved by numerical quadrature methods.
The second term of Eq. 12 is
E(G.) a²
Eg, [log(ww; = - w. diag
w;
+C,
where C constant, EJG.] (E|3igu]) MB+I {a/w diagonal matrix of vector.
5-1 MB+1 and diag(-) indicates the
13

---

Published as a conference paper at ICLR 2021
The updated parameters X and w; can be obtained by setting the gradient of Q to zero. Due to auxiliary variables augmentation, we can see the weights are in a quadratic form in the lower-bound, which leads to analytical expression
x = (N: + K.) /T,
w E, /' (t(t)
where -,= - _o¹ A,(t,w/w;" "dwdtt )dwdt, Σ, is worth noting that numerical quadrature methods need to be applied to intractable integrals above.
Sh (,(t) (t) + diag (a-2))
IV EXPERIMENTAL DETAILS
In this appendix, we elaborate on some experimental details.
SYNTHETIC DATA EXPERIMENTS
For the synthetic data, the intensities and spike times of simulated training and test data shown in Fig. 1. As shown in the experiment of fog-likelihood and running time w.r.. the number of basis functions, the optimal number of basis functions is 4, which are chosen the ground truth: (11,2,3,4) Beta(ā 50, 33 50, scale 6, shift hyperparameter is chosen to be 0.05. As shown in the experiment of log- log-likelihood and running time w.r.t. the number of quadrature nodes, the accuracy is not sensitive to the number of quadrature nodes over 100, so the number of quadrature nodes is set to 2000. The number of EM iterations is set to 200 which large enough for convergence. We plot the estimated influence functions of neurons in Fig. 2. For comparison, we also plot the estimated influence functions of 8 neurons from vanilla multivariate Hawkes processes using the MLE algorithm in Fig. 3 and the functional connectivity graph in Fig. 4. We can both estimated influence functions and functional connectivity graph are far from the ground truth. This demonstrates the necessity of incorporating inhibitive interaction into the model when the Hawkes process is applied in the neuroscience domain. The running time experiment and the fitting and prediction experiment are both conducted for neurons because the baseline models cannot finish in days with 8 neurons because of the curse of dimensionality.
13). By cross validation, the
Intora
Intersitios
TT
intermiity spikes
internsity
Training
Datissel
Test Datasel
Figure 1: The intensities and spike times of neurons in our synthetic training dataset (left) and test dataset (right).
14

---

Published as a conference paper at ICLR 2021
Neaton-1
Hilzeres Neuron2
Iriuernc Neuron
IrBuenco kinctions of Neuioed
InEuence fincrions dhuuss
Irfuercce funollion Newron??
Informnce functions of Neurovd
Actastis
43
Figure 2: The estimated influence functions of all neurons where the estimated d's are close to the ground truth and some ground truth are not labeled since they are all zero (GT=Ground Truth).
Infucnce funcioons of houron-1
infucnce Sunciions of houron-2
Infucnce Suncions houron-3
Infiuance functions houron
Acchation
Actratton
Acthott on
Infuence funciions of heuron-5
funciions heuron-5
ifuuene Suncions houron-7
funciions houren-8
Actrokton
Actration
Acthvation
02
02-
Figure 3: The estimated influence functions of all neurons from vanilla multivariate Hawkes pro- cesses using MLE; some influence functions are not labelled since they are all zero (GT=Ground Truth).
REAL DATA EXPERIMENTS
For the real spike data in cat primary visual cortex areas 17, it contains spike times of 25 simultane- ously recorded neurons. We extract the spike times in the time window [o, 300] (time unit: 100ms) as the training data and [300, 600) as the test data. Both datasets contain approximate 7000 timestamps. The training and test spike trains are plotted in Fig. 5 below.
All hyperparameters are fine tuned in real data experiments. Specifically, the optimal basis function chosen as: Ō Beta(ā 50, 50, scale 10, shift -5.. The hyperparameter a is optimised to be 0.1 by cross validation. The number of quadrature nodes is chosen to be 1000 for which the running time acceptable. The number of EM iterations is set 100 which is large enough for convergence.
15

---

Published as a conference paper at ICLR 2021
Functional Connectivity: Ground Truth Functional Connectivity: MLE Estimation
1.0
Neuron
Neuron
0.2
0.00
Neuron
Neuron
Figure 4: The heat map of functional connectivity among neural population with ground truth (left)
and estimation from vanilla multivariate Hawkes processes (right).
Training Spike Trains
Test Spike Trains
Neurons
50
Figure 5: The training and test spike trains the dataset cat primary visual cortex areas 17.
16


Reviews:
[{"review_id": "i91wjsJwZyM", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The manuscript 'EFFICIENT INFERENCE OF NONPARAMETRIC INTERACTION IN SPIKING-NEURON NETWORKS' develops an expectation-maximization (EM) algorithm, based on an existing trick of introducing additional  latent variables, to obtain the maximum a posteriori (MAP) estimate of parameters of the model. In addition, the authors use a basis function set to express the influence function. Estimating functional interaction among neurons is an important practical problem.\n\n\nPros:\nThey seem to achieve great runtime performance compared to some existing methods, one of which involves MCMC.\nThe basis set allows the authors to express more complex influence functions.\nThe writing is clear.\n\nCons:\nThis is a special application, useful for analyzing neural data.\nThe work is somewhat incremental, although practically useful.\nNot sure that the work provides any insight for a broader class of machine learning problem.\n\nSmaller issues:\nFigure 1 subplots are too small, with (d) and (g) particularly hard to see.\n\nIt might help to explain how to read log likelihood in table 1 and, especially, table 2. Is it obvious that we can compare the columns to each other?", "extra_scores": {}, "extra_text": {}}, {"review_id": "bHtz6o5VLNC", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a novel multivariate nonlinear Hawkes model by modeling the intensity as the product of a upper-bounding intensity, and a sigmoidal function that maps the real-valued latent intensities into a nonnegative real value in [0,1]. Hence, the positive exciting and negative inhibiting interactions among dimensions can be captured by the latent intensities and real-valued influence functions. To further induce sparsity over the influence networks among dimensions, the Laplace prior is used to model the weight parameters of influence functions. A standard EM algorithm is developed to perform inference. Experiments conducted on simulated data demonstrate the novel model can correctly recover the interaction network among neurons, and the influence functions, with less computation time compared with the baselines. They also compared the model with the baselines on real neuron spike data in terms of both train- and test-loglikelihood.\n\n-Quality:  The novel Hawkes model that capture complicated interactions among dimensions, is reasonable to me. The math derivations in Sec 3, appear to be solid. \n\n-Clarity: The motivation and structure of this paper is pretty clear. The authors clearly introduces the main techniques they used from [Adams2009,Donner2017], and their novel Hawkes model.\n\n-Originality & Significance: Although many techniques (sigmoidal function framework, Laplace prior) have been well studied before, they are first reflected in this context of capturing both exciting and inhibiting interactions by a novel multivariate nonlinear Hawkes model. Although the new contributions seem to be incremental, it may has some impacts for modeling complicated neuron interactions. Specifically, the low computational complexity of their novel model is very clear compared with those baselines.\n\n-Pros: A sensible and useful new Hawkes model can capture complicated interactions among dimension, outperform the previous canonical Hawkes models in terms of computation efficiency and test log-likelihood. \n\n-Cons: \nI think the experiments could be more convincing. For synthetic data, the authors only consider one example of neuron interaction networks (Fig.1 a). Although the results show their model correctly recover the ground truth, I am wondering that the authors can provide more examples in supplements, in which the functional connectivity among neurons inferred by the canonical multivariate Hawkes process can be shown to clearly demonstrate why both exciting and inhibiting need to be captured.\n\nIn addition to the train- and test- loglikelihood, can you compare these methods using some other evaluation metrics? In terms of recovering the underlying functional connectivity, does it make sense to compare the estimate functional connectivity graph with the ground truth using AUROC, PR, or F1 scores?\n\nTo my knowledge, it seems to be the first attempt to use the sigmoidal function framework to capture both exciting and inhibiting interactions among dimensions in the multivariate Hawkes model. Nevertheless, I think there also exist some other related multivariate nonlinear point processes can also capture self- and mutual-inhibiting behaviors, e.g., [1]. The authors should consider more strongly related baselines for the comparisons.  \n\n[1] The Multivariate Hawkes Process in High Dimensions: Beyond Mutual Excitation. https://arxiv.org/abs/1707.04928", "extra_scores": {}, "extra_text": {}}, {"review_id": "Hp3gyHjFEu", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper introduces a novel variation on the multivariate Hawkes process that can be fitted efficiently and that allows for interactions to have both excitatory and inhibitory effects. An EM algorithm is provided for fitting the model to empirical data (other methods are discussed in the supplement, but I didn't look at those), and application examples are given on synthetic data and on neuronal spiking data recorded in visual cortex in cat.\n\nThe core contribution of the paper is novel and potentially very useful and relevant to the ICLR community. Unfortunately, the content is not done justice to by the presentation and the quality of the writing, which could be improved or made clearer in several passages. Moreover, a couple of conceptual issues temper my enthusiasm for this work (although they should be easily addressable). All things considered, my initial recommendation is to accept. I have three main points to raise to the attention of the authors (one of which is more of a question/suggestion), and several minor issues to highlight from the point of view of clarity. \n\n###  Main points: \n\n1. From a conceptual standpoint, I am confused by the way in which the neuroscience term \"functional connectivity\" is used. Although this term is very broadly used in different contexts, it is generally meant to convey the notion of effective or statistical connectivity as opposed to direct anatomical or synaptic connectivity. For instance (to cite an influential paper in a subfield not too distant from that of the present work), when describing its model of functional connectivity among retinal cells, Pillow (2008) writes: \"Although it provides an accurate functional description of correlated spike responses, the generalized linear model does not reveal the biophysical mechanisms underlying the statistical dependencies between neurons: coupling does not necessarily imply anatomical connections between cells, but could (for example) reflect dependencies due to shared input noise\". This is, indeed, what the proposed method can measure. But the paper appears to use this concept interchangeably with that of synaptic connectivity (for instance, in section 4.1, first paragraph: \"We aim to identify the synaptic connectivity (functional connectivity) of the neural population…\"), which is problematic. Functional and synaptic connectivity may coincide in carefully constructed synthetic examples, but this is not true in general and certainly not for empirical recordings of neuronal activity. \nThis not just a terminology issue, as it seems to inform at least one technical choice. Indeed, when motivating the choice of basis function sets, the author state \"Although basis functions can be in any form, in order for the weights to represent synapses connection strength, basis functions are chosen to be probability densities with compact support that means they have bounded support [0, $T_\\phi$] and the integral is one\". But by the discussion above the proposed method doesn't tell us anything about the synaptic connectivity between the measured neurons, so the interpretation of the weights as synaptic weights is unfounded. One should rather talk of functional connectivity weights. To see this in a different way, consider the discussion of the real data fit in section 4.2: for a given pair of cells, i and j, the influence function of i on j can take on positive or negative values at different lags (this is indeed one of the main points of the method). But by the argument just discussed on synaptic weights, are we supposed to interpret this as evidence for the fact that neuron i makes multiple synapses onto neuron j (with different time courses), some of which are excitatory and some inhibitory? This makes little sense from a neurobiological/synaptic point of view, while of course still being perfectly fine from the point of view of functional connectivity. I suggest that the author rephrase all their discussions of synaptic connectivity in terms of functional connectivity. \n2. The point above also leads into another question/suggestion for the authors. If you were to drop the requirement that the basis functions used in equations 3 were area-normalized, other basis sets could present themselves as natural choices (I'm bringing this up as it seems that the only justification for the area-normalization requirement was the synaptic weight interpretation, which as I've explained above I believe to be misleading). In particular, B-splines come to mind: like the pdfs used in the paper, they are positive and have finite support. But additionally they offer a very precise characterisation of the class of function they allow to model (the set of splines with the same order, defined over the same knots). This allows one to immediately get an intuitive idea of the complexity of the set of functions over which inference is performed. I don't think this is easily feasible by using, say, Beta pdfs. \n3. I am not sure I understand why the proposed method is called \"nonparametric\". It seems to me that there is a finite set of parameters ($\\mu_i$ and $w_{ijb}$), whose size is fixed from the start and over which inference is performed. The mere fact of looking for a solution over a finite-dimensional space of functions (via eq 3) does not make a model \"nonparametric\". In my understanding, methods are called nonparametric when they involve looking for a solution in an infinite dimensional space (for instance a space of functions constrained only by, say, some smoothness and integrability requirement; but not the finite-dimensional vector space comprising all the linear combinations of a fixed, finite set of basis functions). The authors cite Zhou et al 2013 and state that \"a similar nonparametric scheme\" is implemented in that paper. But if I'm reading that right, that paper also infers the shape of the basis functions ($\\phi$  here, $g$ in Zhou) from data, which involves solving some differential equations during parameter estimation. I see nothing of a similar nature here. And even the authors themselves write things like \"Our goal is to infer the **parameters** i.e. weights and intensity upper bounds…\" (just below equation 5; emphasis mine). So either I am missing something important about the proposed method (and in that case I hope the authors could clarify that in a revised version of the manuscript), or the term \"nonparametric\" is a misnomer here and should be removed from the name of the method and from the paper. As a side effect, it would also make the SNNMHP acronym slightly less intimidating.\n\n### Minor points / clarity \n\n1. Section 2, first paragraph: \"in the meantime\": you probably mean \"at the same time\"? \n2. Same paragraph: \"Apparently, the nonparametric nonlinear multivariate Hawkes processes are a suitable choice for representing the temporal dynamics of mutually excitatory or inhibitory interactions and synapses connectivity of neuron networks.\" I am very confused by that \"apparently\". Do you mean something like \"clearly\", \"evidently\" or \"arguably\" (probably the better option)? \n3. Between eqs 5 and 6: \"As demonstrated in neuron science (Thomson & Bannister, 2003; Sjöström et al., 2001), the synaptic connectivity in cortical circuits is unraveled to be sparse.\". This sentence needs some editing. You probably meant something else than \"unraveled\", and \"neuron science\" is probably meant to be \"neuroscience\". Also the structure of the sentence itself is confusing — it reads as if \"neuron science\" could be only one of many branches of science that could try to have something factual to say about the sparseness of cortical circuits. \n4. Page 4, top: \"exits\" repeated twice in the sentence, is probably meant to be \"exists\"? \n5. Beginning of section 4.1: \"The synthetic neural network contains 4 groups of neurons with a couple of ones in each group\". How about \"The synthetic neural network contains four groups of two neurons each\"? \n6. Section 4.2., paragraph \"Spike train data\": \"Several multi-channel silicon electrode arrays are designed to record simultaneously spontaneous neural activity of multiple isolated single units in anesthetized paralyzed cat primary visual cortex areas 17. It contains spike times of 25 simultaneously recorded neurons.\" Another sentence that needs some editing. I suppose \"it\" is meant to refer to a dataset, but this is never mentioned as the subject of the previous sentence is \"electrode arrays\". \n7. Same section, paragraph \"results\": \"…the real functional connectivity of cortical circuits is unknown…\". This is another point at which the general misuse of the concept of functional connectivity surfaces. Is \"the real functional connectivity\" even a well defined concept? \n8. Similarly, a few lines below: \"most synapses of the 11th neuron are bidirectional\" - I don't think this analysis allows one to conclude anything about direct synaptic contacts between the recorded neurons, given that the 25 neurons are embedded in a vast network. This is an example of a statement that would be better if rephrased in terms of functional connectivity. \n\n\n\n### References\n\nPillow, J. W., Shlens, J., Paninski, L., Sher, A., Litke, A. M., Chichilnisky, E. J., & Simoncelli, E. P. (2008). Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature, 454(7207), 995-999.\n\n----\n\nReview update after revision: my concerns have been addressed.", "extra_scores": {}, "extra_text": {}}, {"review_id": "G1sUwyVW3w", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors here present an extension the Hawkes process to incorporate negative interactions. This allows for inference of excitatory and inhibitory interactions among point-processes using the Hawkes process framework. They present a novel inference procedure for this model using three augmentations to pieces of the model 1) representing the sigmoidal nonlinearity in the model (which is needed to keep the Poisson rates positive) with a mixture of Gaussians wrt a polya-gamma distribution. 2) Using Cambell's theorem to transform the marked Poisson process into gaussian form and 3) using a mixture of gaussians to represent the laplace prior over the weights. The authors fit the model to real and simulated data.\n\nI think this work is an interesting extension to Hawkes processes, the novel EM inference procedure is impressive, and the use on real and simulated data is laudable. \n\nHowever, I think more work needs to be done comparing this model to similar models used to infer E and I connectivity neural data. Since the application here seems to be primarily neural data, I wonder if this model can infer anything that existing models cannot. The first comparison that comes to mind is the GLM. There are recent non-parametric versions (see Dowling et al. 2020) that seem to me to be inferring similar time-varying interactions from the same kind of data as this extended Hawkes process model. I would also like to see more of a discussion about how this model compares to  Apostolopoulou 2019, which, to my knowledge, also uses a Hawkes process extension to capture inhibitory and excitatory connections. The authors should more directly address why their model should be used in place of these other methods.\n\nIt is possible that the contribution here is primarily the speed of the inference procedure. If so, I think more attention should be devoted to this, and more discussion as to why this is much faster than Apostolopoulou 2019. Is it also possible also to compare to an out-of-the-box variational inference method? Some further discussion of how each method scales with number of neurons or length of time series would also be useful.\n\nI'm primarily trying to figure out what the authors want us to take away -- a new kind of model for neuroscience, an inference improvement to Apostolopoulou 2019, or something else. I find the work an interesting extension to the existing body of literature, but a clearer outline of the primary contribution is important, especially given the fact that this seems to be primarily of interest to the neuroscience community, who will want to decide whether or not to use this model on their data.\n\n####EDIT####\nI am satisfied with the reviewers response and will raise my score to a 6.\n\nI think the authors should take care though in speaking practically about their work. The authors highlight the difference of 'excitation-inhibition mixture\" in their response, which is indeed a mathematical difference to the existing work but, at least in the context of neuroscience, I am not sure if this serves any scientific purpose. Do the authors believe this has biophysical meaning or bears a relationship to neural coding? \n\nI do still think some mention of GLMs will be useful as context for models that also infer neuron-neuron E and I interactions.\n\nI appreciate the authors efforts to explicitly highlight the advantages of this work as compared to Apostolopoulou 2019.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This article proposes latent variable augmentation scheme for inference in nonlinear multivariate Hawkes processes. It combines existing approaches (Polya-gamma and sparsity-inducing variables) in a sensible way and is clearly written. Concerns were raised with respect to the comparison to alternative baselines, and answered by the authors. As a result, some reviewers have increased their score, and I recommend acceptance.

Decision: Accept (Poster)