Modelling spine locations on dendrite trees using inhomogeneous Cox point processes

TL;DR

This paper models the distribution of dendritic spine locations on neuron dendrites using inhomogeneous Cox point processes, revealing clustering patterns and proposing a new statistical model to describe spine arrangements.

Contribution

It introduces an inhomogeneous Cox process model for spine locations on dendrites, capturing clustering and spatial patterns not explained by simpler models.

Findings

Some spine data fit inhomogeneous Poisson models

Other data show clustering at larger scales

The Cox process model captures clustering but has high variance in point counts

Abstract

Dendritic spines, which are small protrusions on the dendrites of a neuron, are of interest in neuroscience as they are related to cognitive processes such as learning and memory. We analyse the distribution of spine locations on six different dendrite trees from mouse neurons using point process theory for linear networks. Besides some possible small-scale repulsion, { we find that two of the spine point pattern data sets may be described by inhomogeneous Poisson process models}, while the other point pattern data sets exhibit clustering between spines at a larger scale. To model this we propose an inhomogeneous Cox process model constructed by thinning a Poisson process on a linear network with retention probabilities determined by a spatially correlated random field. For model checking we consider network analogues of the empirical $F$-, $G$-, and $J$-functions originally introduced for inhomogeneous point processes on a Euclidean space. The fitted Cox process models seem to catch the clustering of spine locations between spines, but also posses a large variance in the number of points for some of the data sets causing large confidence regions for the empirical $F$- and $G$-functions.

Figures

Tables

Table 1. Table 1: Number of spines, length, and intensity estimates for the main and side branches separately, where the intensity estimates are for the parametric model ( 4.1 ) given in Section  4.1 .

Dendrite	$n_m$	$n_s$	$\left|L_m\right|$	$\left|L_s\right|$	${\hat{\rho }}_m$	${\hat{\rho }}_s$
1	51	72	$212\mu m$	$202\mu m$	0.240	0.356
2	36	145	$204\mu m$	$430\mu m$	0.176	0.337
3	69	63	$211\mu m$	$202\mu m$	0.328	0.312
4	33	308	$225\mu m$	$652\mu m$	0.134	0.477
5	34	83	$286\mu m$	$450\mu m$	0.119	0.184
6	30	62	$178\mu m$	$250\mu m$	0.168	0.248

Table 2. Table 2: Estimates of ρ Y , m subscript 𝜌 𝑌 𝑚 \rho_{Y,m} , ρ Y , s subscript 𝜌 𝑌 𝑠 \rho_{Y,s} , σ 2 superscript 𝜎 2 \sigma^{2} , and β 𝛽 \beta for each spine data set.

Dendrite	${\hat{\rho }}_{Y,m}$	${\hat{\rho }}_{Y,s}$	${\hat{\sigma }}^2$	$\hat{\beta }$
1	$0.312$	$0.463$	$0.686$	$0.037$
2	$0.275$	$0.525$	$1.427$	$0.020$
3	$0.328$	$0.312$	$5.170\times {10}^{-8}$	$30.178$
4	$0.197$	$0.701$	$1.159$	$0.010$
5	$0.266$	$0.413$	$4.023$	$0.030$
6	$0.322$	$0.474$	$2.662$	$0.013$

Table 3. Table 3: Overview of runs made in the simulation study for the minimum contrast procedures. Here ( σ ∗ 2 , β ∗ ) subscript superscript 𝜎 2 subscript 𝛽 (\sigma^{2}_{*},\beta_{*}) denote the start parameters for the optimisation algorithm and b ​ w 𝑏 𝑤 bw is the automatically selected bandwidth used to calculate g ^ ^ 𝑔 \hat{g} in spatstat .

Run no.	${\sigma }^2$	$\beta$	${\rho }_{Y,m}$	${\rho }_{Y,s}$	$p$ for MCE-$g$ (MCE-$K$)	$r_l$	$r_u$	$({\sigma }_{\ast }^2,{\beta }_{\ast })$
1	5	0.1	0.8	1.2	1 (0.25)	0	30	(0.5, 0.5)
2	5	0.1	0.8	1.2	0.5 (0.5)	0	30	(0.5, 0.5)
3	5	0.1	0.8	1.2	1 (0.25)	0	50	(0.5, 0.5)
4	5	0.1	0.8	1.2	1 (0.25)	0	20	(0.5, 0.5)
5	5	0.1	0.8	1.2	1 (0.25)	0	30	(3, 0.2)
6	5	0.1	0.8	1.2	1 (0.25)	0	30	(0.2, 3)
7	5	0.1	0.3	0.7	1 (0.25)	0	30	(0.5, 0.5)
8	5	0.1	1	1	1 (0.25)	0	30	(0.5, 0.5)
9	5	0.5	0.8	1.2	1 (0.25)	0	30	(0.5, 0.5)
10	5	1	0.8	1.2	1 (0.25)	0	30	(0.5, 0.5)
11	1	0.1	0.8	1.2	1 (0.25)	0	30	(0.5, 0.5)
12	1	0.5	0.8	1.2	1 (0.25)	0	30	(0.5, 0.5)
13	5	0.1	0.8	1.2	1 (0.25)	$2\ast bw$	30	(0.5, 0.5)
14	5	0.1	0.8	1.2	1 (0.25)	$0.5\ast bw$	30	(0.5, 0.5)
15	5	0.1	0.8	1.2	1 (0.25)	$2$	30	(0.5, 0.5)
16	5	0.1	0.8	1.2	1 (0.25)	$0.5$	30	(0.5, 0.5)

Equations

$$\mathrm{E}n(X\cap B)={\int }_B\rho (u){\mathrm{d}}_{L}u<\infty ,$$

$$\mathrm{E}\{n(X\cap A)n(X\cap B)\}={\int }_A{\int }_{B}g(u,v)\rho (u)\rho (v){\mathrm{d}}_{L}u{\mathrm{d}}_{L}v<\infty .$$

$$K(r)={\int }_0^r{g}_0(s)\mathrm{d}s,r\ge 0.$$

$$\hat{F}(r)$$

$$\hat{G}(r)$$

$$\hat{J}(r)$$

$$\rho (u)={\rho }_m^{\mathbb{I}(u\in L_m)}{\rho }_s^{\mathbb{I}(u\in L_s)},u\in L,$$

$${\hat{\rho }}_m=\frac{n_m}{|L_m|},{\hat{\rho }}_s=\frac{n_s}{|L_s|},$$

$$\Pi(u)=\exp\Bigl{\{}-\frac{\sigma^{2}}{2}\sum_{j=1}^{k}Z_{j}^{2}(u)\Big{\}},\qquad u\in L,$$

$$c(u,v)=\exp \{-\beta d_L(u,v)\},u,v\in L,$$

$$\mathrm{E}\{\Pi (u)\Pi (v)\}=\{(1+{\sigma }^2{)}^2-({\sigma }^2{)}^{2}c(u,v{)}^2{\}}^{-k/2},$$

$$\rho (u)=(1+{\sigma }^2{)}^{-k/2}{\rho }_Y(u),u\in L,$$

$$g(u,v)={\left\{\frac{(1+{\sigma }^2{)}^2}{(1+{\sigma }^2{)}^2-({\sigma }^2{)}^{2}c(u,v{)}^2}\right\}}^{k/2},u,v\in L,u\ne v.$$

$${\rho }_Y(u)={\rho }_{Y,m}^{\mathbb{I}(u\in L_m)}{\rho }_{Y,s}^{\mathbb{I}(u\in L_s)},u\in L,$$

$$(\hat{\sigma}^{2},\hat{\beta})=\underset{(\sigma^{2},\beta)}{\arg\min}\int_{r_{l}}^{r_{u}}\bigl{\{}\hat{T}(r)^{p}-T_{(\sigma^{2},\beta)}(r)^{p}\bigr{\}}^{2}\,\mathrm{d}r,$$

$$Y=\{u\in \tilde{Y}:d_L(u,v)>h\text{ for all }v\in \tilde{Y}\backslash \{u\}\},$$

$$K(r)={\int }_0^r{\left\{1-\frac{\exp (-2\beta t)}{{\left(1+1/{\sigma }^2\right)}^2}\right\}}^{-k/2}\mathrm{d}t.$$

$$\frac{1}{\beta }\left\{\log \left(\sqrt{{\mathrm{e}}^{2\beta r}-\alpha }+{\mathrm{e}}^{\beta r}\right)-\log \left(\sqrt{1-\alpha }+1\right)\right\}$$

$$\frac{1}{2\beta }\log \left(\frac{{\mathrm{e}}^{2\beta r}-\alpha }{1-\alpha }\right)$$

$$\displaystyle\frac{1}{\beta}\biggl{\{}$$

$$\displaystyle\frac{1}{2\beta}\Biggl{\{}\frac{\alpha}{\alpha-\mathrm{e}^{2\beta r}}-\frac{\alpha}{\alpha-1}+\log\Bigl{(}\frac{\mathrm{e}^{2\beta r}-\alpha}{1-\alpha}\Bigr{)}\Biggr{\}}$$

$$\displaystyle\frac{1}{\beta}\biggl{\{}\frac{\alpha+\mathrm{e}^{-\beta r}(\mathrm{e}^{2\beta r}-\alpha)^{3/2}\log(\sqrt{\mathrm{e}^{2\beta r}-\alpha}+\mathrm{e}^{\beta r})-\frac{4}{3}\mathrm{e}^{2\beta r}}{\sqrt{1-\alpha\mathrm{e}^{-2\beta r}}(\mathrm{e}^{2\beta r}-\alpha)}-\log(\sqrt{1-\alpha}+1)+\frac{\frac{4}{3}-\alpha}{(1-\alpha)^{3/2}}\biggr{\}}$$

$$CL({\sigma }^2,\beta )=$$

$$w(u,v)=\mathbb{I}(d_L(u,v)\le r_0)$$

$$\nabla CL($$

$$=\sum _{u,v\in x}^{\ne }w(u,v)\frac{\nabla g(u,v)}{g(u,v)}-{\int }_L{\int }_{L}w(u,v)\hat{\rho }(u)\hat{\rho }(v)\nabla g(u,v){\mathrm{d}}_{L}u{\mathrm{d}}_{L}v.$$

$$\displaystyle w(u,v)=\mathbb{I}\biggl{(}\frac{|g(u,v)-1|}{M(u,v)}>\epsilon\biggr{)},$$

$$w(u,v)=\{\begin{array}{ll}\exp \left[1/\{h(u,v{)}^2-1\}\right]& \text{ for }-1\le h(u,v)\le 1,\\ 0& \text{ else,}\end{array}$$

$${\int }_L{\int }_{L}f(u,v){\mathrm{d}}_{L}u{\mathrm{d}}_{L}v=\sum _{i,j}{\int }_{L_i}{\int }_{L_j}f(u,v){\mathrm{d}}_{L}u{\mathrm{d}}_{L}v.$$

$${\int }_{L_i}{\int }_{L_j}f(u,v){\mathrm{d}}_{L}u{\mathrm{d}}_{L}v$$

Full text

\keywords

Empirical summary functions; linear networks; random fields; thinned point process.

1 Introduction

Point patterns on linear networks arise in a broad range of fields, where the network for example represents roads, a river network, or a dendrite tree. This paper focuses on the latter type of data: the left panel in Figures 1–2 shows six linear networks each representing a dendrite tree from a mouse neuron grown in vivo. On the dendrites small protrusions called spines are found that among other things help transmitting electrical signals to the soma. In neuroscience, the behaviour of spines is of interest as changes can be linked to changes in cognitive processes. The spine locations can be viewed as a point pattern on the dendrite tree and thus analysed using point process theory for linear networks. These six spine point pattern data sets will hereafter be referred to as the ‘spine data’.

Over the last two decades, methods for analysing point patterns on linear networks have been developed. Particularly, a network analogue of Ripley’s $K$-function was first presented in Okabe and Yamada, (2001) and later modified and extended to the inhomogeneous case in Ang et al., (2012). When defining the $K$-function, Ang et al., (2012) required that the underlying point process model fulfils an invariance property called second-order pseudo-stationarity (an analogue to second-order intensity-reweighted stationarity as introduced in Baddeley et al., , 2000). This property is fulfilled whenever the pair correlation function is isotropic, i.e. when it only depends on the shortest path distance. Baddeley et al., (2017) showed that certain constructions, e.g. special types of Cox point processes that lead to point processes in the Euclidean space with an isotropic pair correlation function rarely result in second-order pseudo-stationary point processes when adapted to linear networks. Even without the requirement of pseudo-stationarity, there are currently only a limited number of point process models available for linear networks. For point processes on directed acyclic linear networks, Rasmussen and Christensen, (2019) presented both regular and clustered models defined by a generalisation of the conditional intensity function for temporal point processes. Anderes et al., (2020) supplied a list of valid isotropic covariance functions for connected linear networks that can be used to construct Cox point processes, particularly log Gaussian Cox processes (LGCPs; see also Møller et al., , 1998).

Only few studies use point process theory to analyse the behaviour of spines: treating the dendrite tree as a directed tree, Rasmussen and Christensen, (2019) analysed the spine data for dendrite 4 (visualised in the first row of Figure 2). The distribution of spines (and their shape) have further been investigated using point process theory in Jammalamadaka et al., (2013) (testing a homogeneous Poisson process model) and Baddeley et al., (2014) (using multitype Poisson process models to account for the shape classification) for in vitro grown neurons. Based on the network $K$-function, Jammalamadaka et al., (2013) concluded that a homogeneous Poisson process model seems adequate to describe the spine locations. However, Jammalamadaka et al., (2013) also stated that their results for the in vitro setting are unlikely to hold in an in vivo setting.

Instead of Poisson process models, this paper suggests a new class of Cox process models on a linear network. Such a model applies for an undirected graph and is not a LGCP, but its construction still exploits a Gaussian random field so that the covariance functions from Anderes et al., (2020) become useful. Moreover, seemingly for the first time in connection to point process model fitting on linear networks, we demonstrate the use of minimum contrast and composite likelihood estimation procedures. Finally, we introduce new empirical summary functions and demonstrate their usefulness for model checking.

Our new empirical summary functions are analogous to the empirical $F$-, $G$-, and $J$-functions in van Lieshout, (2011) for inhomogeneous point processes on a Euclidean space (Okabe and Sugihara, , 2012, introduced a version of $F$- and $G$-functions but for homogeneous point processes on a linear network). We also consider empirical intensity functions, and after completing the first version of our paper we became aware of later appearing papers on intensity estimation (Rakshit et al., , 2019) and summary statistics (Cronie et al., , 2019) for inhomogeneous point processes on a linear network. In particular, Cronie et al., (2019) introduced other kinds of empirical $F$-, $G$-, and $J$-functions which may be preferred from ours since their functions account for the geometry of the linear network.

The paper is organised as follows. The spine data is described in more detail in Section 2 along with the general notion of a linear network. In Section 3 we discuss existing as well as our new summary functions for point processes on linear networks; these are used for analysing the spine locations in Section 4. We initially suggest to model the spine locations by an inhomogeneous Poisson process model in Section 4.1, but due to clustering between spines we propose in Section 4.2 an inhomogeneous Cox process model. Lastly, we discuss in Section 5 possible extensions and future research directions.

2 Dendritic spine data

The spine data origin from six apical dendrite trees corresponding to two neurons from each of three different mice, with each mouse neuron grown in vivo. The numbering of the dendrites is as follows: dendrite 1 and 2 come from mouse no. 1; dendrite 4 and 5 from mouse no. 2; and dendrite 3 and 6 from mouse no. 3. The first and middle columns of Figures 1–2 show different ways of viewing each dendrite as a linear network (more details are given below). Specifically, a linear network is a union $L=\bigcup_{i=1}^{N}L_{i}$ of a finite number $N$ of line segments $L_{i}\subset\mathbb{R}^{d}$, $d\geq 2$, with finite length and intersecting only at the end points. A linear network may also be viewed as a graph consisting of a set of vertices and a set of weighted edges, where the edges coincide with the line segments $L_{1},\dots,L_{N}$, the vertices correspond to the end points of these line segments, and the weight of an edge is the length of the corresponding line segment.

Throughout this paper we assume that $L$ is a connected set and the distance between two points $u,v\in L$ is measured by the shortest path distance and is denoted by $d_{L}(u,v)$. For the spine data, the linear network $L$ of a dendrite is a tree, meaning that there is only one path between any pair of points in $L$. Naturally, in other applications more complicated networks than a tree occur in which case we may need to take more care when choosing the distance metric $d_{L}$ (see Section 5 for details).

The linear networks visualised in the left column of Figures 1–2 are approximations of the underlying apical dendrite trees. The vertices of each of the linear networks are described by three-dimensional coordinates which represent a spine location or another point chosen to obtain the approximation. For each dendrite tree, we talk about two subsets: the main branch and the side branches. Main branch refers to the tree’s stem, while side branches constitute the rest of the tree. The left column of Figures 1–2 shows which parts of the approximations of the underlying apical dendrite trees belong to the main branch and which to the side branches.

To utilise the functionalities of the R-package spatstat (Baddeley et al., , 2015), each of the linear networks has been transformed into a simpler network embedded in $\mathbb{R}^{2}$ in the following way. First, we reposition any sequence of edges which are consecutively connected by vertices of degree two (such edges will either be contained in the main branch or in the same side branch) with a single edge/line segment whose weight is the sum of the old edge weights. This procedure straightens out ‘kinks’ in any connected subset of the network containing only such line segments, as the angle between a pair of neighbouring edges becomes $180^{\circ}$. Next, any line segments meeting at a vertex of degree higher than two can be repositioned such that the entire network is contained in a plane with no overlap between line segments (except at the end points). Thereby the network has naturally been embedded in $\mathbb{R}^{2}$. Note that this transformation relies on the fact that the linear network forms a tree, and that it preserves the distance between any pair of points in the network. Most importantly, the transformation allows us without any loss of information to consider the spine locations as a point pattern on the simplified and embedded linear network when in Section 4 we are dealing with

• •

non-parametric kernel estimates of the intensity, since they will be unchanged under the performed transformation,

• •

specific Poisson and Cox process models and statistical tools, since as we shall see they will not directly depend on the three-dimensional coordinates but only on

• –

how many points in a spine point pattern data set will be on the main branch and on the side branches, respectively,

• –

the distances between pairs of points.

The simplified and embedded versions of the networks are shown in the middle column of Figures 1–2 along with the spine locations. It is also shown which parts of the simplified embedded trees belong to the main branch and which to the side branches. In the following, $L=L_{m}\cup L_{s}$ refers to one of the six simplified and embedded linear networks in Figures 1–2, where $L_{m}$ is the main branch, and $L_{s}$ is the union of the side branches. Further, $n_{m}$ and $n_{s}$ denote the number of spines on $L_{m}$ and $L_{s}$, respectively. Lastly, we let $|B|$ denote the size of $B\subseteq L$ or more precisely the total length of the (partial) line segments constituting $B\subseteq L$; note that $|L|=|L_{m}|+|L_{s}|$. Table 1 summarises the number of spines and sizes for each dendrite tree.

3 Point processes on linear networks

The six spine point pattern data sets are modelled as realisations of six point processes defined on the six dendrite trees. In general, by a point process $X$ on a linear network $L$ we mean a random finite subset of $L$; we use this generic notation throughout this paper. In this section we consider summary functions useful for analysing point processes on linear networks, including the introduction of new empirical summary functions.

3.1 Summary functions for first and second-order moment properties

We assume that $X$ has intensity $\rho$, that is, for $B\subseteq L$,

[TABLE]

where $n(X\cap B)$ is the number of points from $X$ falling in $B$ and $\mathrm{d}_{L}$ denotes integration with respect to one-dimensional arc-length along $L$. Intuitively, $\rho(u)\,\mathrm{d}_{L}u$ is the probability of $X$ having a point in an infinitesimal small subset of $L$ that contains $u$ and has size $\mathrm{d}_{L}u$. If the intensity $\rho(\cdot)\equiv\rho$ is constant, we say that $X$ is homogeneous; otherwise $X$ is said to be inhomogeneous. In case of homogeneity, $\rho$ is the expected number of points per unit length.

We also assume that $X$ has pair correlation function $g$, that is, for disjoint $A,B\subset L$,

[TABLE]

We can interpret $g(u,v)\rho(u)\rho(v)\,\mathrm{d}_{L}u\,\mathrm{d}_{L}v$ as the joint probability that two infinitesimal small regions around $u$ and $v$ of size $\mathrm{d}_{L}u$ and $\mathrm{d}_{L}u$, respectively, each contains a point from $X$.

We say that the pair correlation function is isotropic if it only depends on the shortest path distance, that is, $g(u,v)=g_{0}\{d_{L}(u,v)\}$ where $g_{0}$ is a non-negative function. When $X$ has an isotropic pair correlation function, the (geometrically corrected network) $K$-function introduced by Ang et al., (2012) can be expressed as

[TABLE]

Alternatively, $K(r)$ may be written as an expectation with respect to a Palm distribution, see (Ang et al., , 2012, Theorem 3). If the $K$-function or the pair correlation function is expressible on closed form, we can use a minimum contrast or composite likelihood procedure to estimate the model parameters; this is described further in Section 4.2.2.

The simplest point process model is a Poisson process, which is characterised by that $n(X)$ follows a Poisson distribution with mean given by (3.1) with $B=L$ and further that the points of $X$ conditioned on $n(X)$ are independent and identically distributed, with density proportional to $\rho$. For a Poisson process, $g\equiv 1$ and $K(r)=\nobreak r$.

3.2 New empirical summary functions

For estimating the pair correlation function and the $K$-function we follow Ang et al., (2012). These empirical summary functions can be used in minimum contrast or composite likelihood estimation procedures as well as for model checking. Obviously, if the $K$-function or pair correlation function have been used to fit the model, neither should be used to check the adequacy of the model. Due to the shortage of summary functions for point processes on linear networks, we may let a simple visual comparison of the observed point pattern and simulations from the fitted model serve as a model check. It is needless to say that a more rigorous model checking would be preferred.

Therefore, we now introduce three purely empirical summary functions. These are obtained by modifying the empirical $F$-, $G$-, and $J$-functions for inhomogeneous point patterns on a Euclidean space (introduced by van Lieshout, , 2011) to linear networks. The modification simply consists of replacing the Euclidean space with the linear network, introducing the shortest path distance instead of the Euclidean distance, and adapting the notion of an eroded set to linear networks. The functions are then defined as follows. Assume that the intensity $\rho$ is known or has been estimated by $\hat{\rho}$ and that $\bar{\rho}=\inf_{u\in L}\tilde{\rho}(u)>0$, where either $\tilde{\rho}=\rho$ or $\tilde{\rho}=\hat{\rho}$. For $r\geq 0$, let $L_{\ominus r}$ consist of the points in $L$ with distance greater than $r$ to any vertex of $L$ with degree one. Furthermore, let $H\subset L$ be a finite ‘lattice’. For a tree network the points of $H$ may be chosen equidistant, while on more complicated networks we can e.g. choose $H$ to be a collection of points that are equidistant within each line segment. For an observed point pattern $X=x$, the empirical summary functions $\hat{F}$, $\hat{G}$, and $\hat{J}$ are then defined for $r\geq 0$ by

[TABLE]

where we restrict attention to $r$-values small enough to ensure that $\#(H\cap L_{\ominus r})>0$ for $\hat{F}(r)$, $\#(x\cap L_{\ominus r})>0$ for $\hat{G}(r)$, and $\hat{F}(r)<1$ for $\hat{J}(r)$. Further, $r$ should not be chosen larger than a reasonable amount of the network remains when considering $L_{\ominus r}$.

In Ang et al., (2012), the $K$-function and its empirical estimate include a factor that corrects for the network geometry, such that its shape can be compared for point patterns on different networks. As it was not obvious to us how to extend such a correction to $\hat{F}$, $\hat{G}$, and $\hat{J}$ in a meaningful way, our definitions in (3.3)–(3.5) do not correct for the network geometry. Further, we do not have any theoretical counterpart to $\hat{F}$, $\hat{G}$, and $\hat{J}$ and therefore their shapes alone can in general not be used to conclude anything about e.g. the presence of regularity or clustering. However, $\hat{F}$, $\hat{G}$, and $\hat{J}$ are still useful tools for providing a so-called global rank envelope; this is a confidence region for a given test function obtained from simulations under a fitted model (for details, see Myllymäki et al., , 2017). In a global rank envelope procedure, the shape of the test function for the data is compared to that of the simulations and Myllymäki et al., (2017) discussed how this provides a test and an interval with lower and upper bounds given by a liberal and a conservative $p$-value, respectively.

4 Modelling spine locations

In this section each of the six data sets is analysed with the aim of finding a model that adequately describe the spine locations.

4.1 Fitting a Poisson process model

The simplest model we can propose is a Poisson process. To investigate the behaviour of the spine intensity, we calculated the non-parametric intensity estimate suggested by McSwiggan et al., (2016) using the density.lpp function from the spatstat-package with Scott’s rule of thumb modified to linear networks for choosing the bandwidth (Rakshit et al., , 2019); the resulting estimates are seen in the right panel of Figures 1–2, where the bandwidth is 18, 20, 18, 22, 23, 14 for dendrite $1,\dots,6$, respectively. Note that the estimated spine intensity tends to be higher on the side branches than on the main branch (except perhaps for dendrite 3) and it seems plausible to assume a constant intensity value $\rho_{m}$ and $\rho_{s}$ on the main and side branches, respectively (with different $(\rho_{m},\rho_{s})$-values for the six dendrites). Therefore, recalling the notation in Section 2, we assume that

[TABLE]

where $\mathbb{I}(\cdot)$ denotes the indicator function. Naturally, we could also propose inhomogeneous Poisson process models involving other covariates (than $\mathbb{I}(u\in L_{m})$ and $\mathbb{I}(u\in L_{s})$) on the network. As no quantities beside the location of the spines have been recorded, the only covariates available are those that we can create within the network, such as distance to a vertex of degree higher than two. However, as we see no clear relation between such covariates and the intensity function we restrict ourselves to consider the Poisson process model specified by (4.1).

The maximum likelihood estimates of the intensity parameters in (4.1) are easily found and given by

[TABLE]

cf. the notation in Section 2. These estimates are shown in Table 1.

To test whether the proposed inhomogeneous Poisson process model adequately describes the spine locations, we performed global rank envelope tests using $K$ as test function, cf. Section 3.2. Results from these tests are shown in Figure 3. For all dendrites the conservative $p$-value is below $3\%$, suggesting that the fitted Poisson process models do not describe the spine locations adequately. Specifically, the empirical $K$-functions for dendrite 2, 4, 5, and 6 fall above the global rank envelopes for certain $r$-values, indicating that the spines tend to cluster at these distances. Further, for all six spine data sets the empirical $K$-function falls below the global rank envelope for small $r$-values, indicating a small-scale repulsion between spines. For dendrite 1 and 3, this small-scale repulsion is the only deviation from the Poisson process model revealed by the global rank envelope test. Disregarding the small distances ($r<$1\text{,}\mathrm{\SIUnitSymbolMicro m}$$) for the global rank envelope test with $K$ as test function, does not change the $p$-intervals significantly for dendrite 2, 4, 5, and 6. For dendrite 1 and 3 on the other hand, the $p$-intervals change from $(0.024,0.040)$ to $(0.048,0.060)$ and from $(0,0.019)$ to $(0.160,0.168)$, respectively, giving (most clearly for dendrite 3) no evidence against the proposed Poisson process model. Global rank envelopes with a concatenation of $\hat{F}$, $\hat{G}$, and $\hat{J}$ as test function, and where distances less than $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ were disregarded, are shown in Figure C.1 in Appendix C; these do not provide any evidence against the Poisson process model for dendrite 1 and 3 either.

As the physical scale of the spine data is quite small (the dendrites range in size from $412\text{\,}\mathrm{\SIUnitSymbolMicro m}$ to $876\text{\,}\mathrm{\SIUnitSymbolMicro m}$, cf. Table 1) and as there is uncertainty in the exact choice of the point representing a spine’s location, we must expect some degree of imprecision and therefore we may not want to put too much value into the observed small-scale repulsion. In the following we will not take the small-scale repulsion into account but rather focus on modelling the large scale clustering.

4.2 Fitting a Cox process model

In addition to inhomogeneity in the location of spines, Figures 1–2 also shows a tendency to large bare areas where no spines occur. To model this behaviour we considered a point process model introduced by Lavancier and Møller, (2016) for the Euclidean space, which is easily adapted to a linear network $L$. The point process is constructed as a thinning of a Poisson process $Y$ on $L$ with intensity function $\rho_{Y}$, where the retention probabilities are determined by a random field $\Pi=\{\Pi(u):u\in L\}$ that may be spatially correlated. That is, the point process is given by $X=\{u\in Y:\Pi(u)\geq R(u)\}$, where $R=\{R(u):u\in L\}$ consist of independent uniform random variables on $[0,1]$, and where $Y$, $\Pi$, and $R$ are independent. Thus $X$ is a Cox point process driven by the random field $\Lambda=\{\rho_{Y}(u)\Pi(u):u\in L\}$. We let

[TABLE]

where $k\in\{1,2,\dots\}$ and $\sigma^{2}>0$ are parameters, and $Z_{1},\dots,Z_{k}$ are IID zero-mean unit-variance Gaussian random fields with correlation function $c$. If $c(u,v)=c_{0}\{d_{L}(u,v)\}$ depends only on the shortest path distance, we say that $c$ is isotropic; see Anderes et al., (2020) for a list of isotropic correlation functions for linear networks. For the spine data, we considered the exponential correlation function, that is,

[TABLE]

which is a valid correlation function for any parameter value $\beta>0$ and any tree network but not necessarily for other kinds of linear networks (see Section 5 for a comment on this).

We have that $\mathrm{E}\,\Pi(u)=(1+\sigma^{2})^{-k/2}$ and

[TABLE]

implying that $X$ has intensity

[TABLE]

and pair correlation function

[TABLE]

If $c$ is isotropic, then $g$ is isotropic and the $K$-function can be expressed by (3.2). Closed form expressions of the $K$-function are given in Appendix A for $c$ equal to the exponential correlation function and $k=1,\dots,5$.

Note that for each $u\in L$, $\sigma^{2}\sum_{j=1}^{k}Z_{j}^{2}(u)$ in (4.3) follows a $\sigma^{2}\chi^{2}$-distribution with $k$ degrees of freedom. That is, the skewness decreases as $k$ increases, while the range is stretched/compressed depending on the value of $\sigma^{2}$. The larger $\sigma^{2}$ is, the more $Y$ is thinned to obtain $X$ and also the more variation in the thinning probabilities. The pair correlation function in (4.6) is an increasing function of both $\sigma^{2}$ and $k$. When $c$ is given by (4.4), $\beta$ controls the correlation of the retention probabilities: the smaller $\beta$, the longer range of correlation in $\Pi$ and thus larger coherent bare/populated areas in $X$. Finally, the pair correlation function decreases towards 1, when $\beta$ increases.

For the spine data, we still want to model the intensity function by (4.1), which by (4.5) requires $Y$ to have a similar intensity structure, that is,

[TABLE]

for non-negative parameters $\rho_{Y,m}$ and $\rho_{Y,s}$.

4.2.1 Simulation

To perform model checking with global rank envelopes or to carry out simulation studies, we need to be able to simulate point patterns from the model of interest. Fortunately, it is straightforward to simulate a point pattern on $L$ from the proposed Cox process model by the following three steps:

1. (a)

Simulate a discretised version of the random field $\Pi$ by first simulating the independent Gaussian random fields $Z_{j}$, $j=1,\dots,k$, at chosen grid locations along the network and then transforming the random fields according to (4.3) to obtain the retention probabilities. 2. (b)

Simulate a point pattern $y$ from a Poisson process on $L$ with intensity $\rho_{Y}$. 3. (c)

Thin $y$ using the retention probabilities simulated in (a).

4.2.2 Estimation procedure

In the following we describe a procedure for estimating the model parameters of the proposed Cox process model where. For specificity we consider the case where $c$ is given by (4.4) and where $\rho_{Y}$ is given by (4.7).

To begin we assume that $k$ is known, whereas the remaining parameters are estimated through a two-step procedure (Waagepetersen, , 2007; Waagepetersen and Guan, , 2009). In short, we first estimate $(\rho_{m},\rho_{s})$ and then plug in these estimates in a second-order procedure where $(\sigma^{2},\beta)$ is estimated. Lastly, an estimate of $(\rho_{Y,m},\rho_{Y,s})$ can be found by using (4.5).

First, to estimate $(\rho_{m},\rho_{s})$ we use the first order composite likelihood (Waagepetersen, , 2007) which simply corresponds to a Poisson likelihood yielding the estimates in (4.2).

Second, as we know explicit formulas for the pair correlation and $K$-function, we can estimate $(\sigma^{2},\beta)$ using a minimum contrast procedure (Guan, , 2009; Diggle, , 2014) or a second-order composite likelihood approach (Waagepetersen, , 2007; Lavancier et al., , 2018). The latter is not considered here but described in Appendix B, where results from a simulation study comparing the two approaches also can be found. The simulation study suggests that the minimum contrast procedure provides more meaningful estimates with less bias and variance than the second-order composite likelihood.

For a chosen summary function $T_{(\sigma^{2},\beta)}$ which depends on $(\sigma^{2},\beta)$, the minimum contrast estimate of $(\sigma^{2},\beta)$ is given by

[TABLE]

where $0\leq r_{l}<r_{u}$ and $p>0$ are user-specified tuning parameters, and $\hat{T}$ is an empirical estimate of the summary function. In our case, $\smash{T_{(\sigma^{2},\beta)}}$ is given by $K$ or $g_{0}$, and $\hat{T}$ is given as in Ang et al., (2012). A frequently seen choice is $r_{l}=0$, while general recommendations of $r_{u}$ and $p$ can be found in Guan, (2009) and Diggle, (2014) for point patterns on the Euclidean space.

The minimum contrast procedure can easily be extended to include estimation of $k$ too, but a simulation study indicated that it may be difficult to estimate $\sigma^{2}$ and $k$ simultaneously as an increase in $k$ seemingly can be balanced out by an increase in $\sigma^{2}$. In practice we may therefore simply make a choice of $k$; for simplicity we chose $k=1$ in the following. Note that for both $T=K$ and $T=g_{0}$ the estimate $\hat{T}$ in (4.8) depends on the intensity (Ang et al., , 2012); here we simply plug-in the estimated intensity obtained in the first step of the estimation procedure.

A drawback of using $T=g_{0}$ is the need of choosing a bandwidth for the non-parametric kernel estimate $\hat{g}_{0}$ presented in Ang et al., (2012). However, the simulation study in Appendix B suggests that using $T=g_{0}$ with the default bandwidth and kernel from the spatstat-package generally yield estimates with lower bias and variance than $T=K$ when fitting the proposed Cox process model. This is consistent with results from a simulation study in Lavancier and Møller, (2016) for point processes on a Euclidean space.

In the simulation study found in Appendix B, we also investigated how different choices of $r_{l}$, $r_{u}$, and $p$ affect the estimates of $\sigma^{2}$ and $\beta$ given by (4.8). We observed that the choice of $r_{u}$ often is a matter of trade-off between bias and variance: a large value of $r_{u}$ may entail a large bias, while a small $r_{u}$ often leads to a greater variance of the estimates. The best choice of $r_{u}$ also seems to be quite depending on what the true underlying model parameters are. For example, a larger range of correlation in the retention probabilities, that is, a smaller value of $\beta$, requires a larger $r_{u}$. Naturally we should also take the size of the network into consideration when choosing $r_{u}$. In the simulation study we found that $r_{l}=0$ gives the best estimates in terms of bias and variance, and that $p=1$ or $p=1/2$ behave equally well for $T=g_{0}$, while $p=1/4$ is preferred over $p=1/2$ for $T=K$. For parameter values yielding models close to the Poisson process model, that is, when $\sigma^{2}$ is close to zero or $\beta$ is large, the estimation procedures were not very successful regardless of the tuning parameters. This does not come as a surprise as many combinations of $\sigma^{2}$- and $\beta$-values yield similar Poisson processes. Lastly, the estimation procedure seems quite stable with respect to the choice of start parameter values for the optimisation algorithm (optim in R) used to minimize (4.8).

4.2.3 Model fit and model check

The Cox process model was fitted to each of the spine data sets using the two-step procedure with $k=1$ fixed, cf. Section 4.2.2. For the minimum contrast procedure we let $T=g_{0}$, $p=1$, and $r_{l}=0$ in accordance with the simulation results discussed in Section 4.2.2. Further, we initially let $r_{u}=15$ and obtained a set of initial parameter estimates for each data set. Then we performed a small simulation study based on 500 simulations from the initially fitted models to investigate which of $r_{u}=15$, $30$, $50$ results in the best estimates (in terms of bias and variance) for these specific models. For dendrite 1, 3, and 6, the initially fitted models are close to the Poisson process case, and as a consequence the model parameters are hard to estimate regardless of the choice of $r_{u}$. However, for dendrite 2, 4, and 5, it seems that $r_{u}=50$ is the best choice. Using $r_{u}=50$ for all six data sets, we obtained the parameter estimates in Table 2. The fitted model for dendrite 3 is practically a Poisson process model (in consistency with the conclusions made in Section 4.1) while the remaining fitted models are not. In Figure C.2 in Appendix C, one simulation from each of the fitted random fields $\Pi$ is shown to illustrate the behaviour of the retention probabilities. For example, $\hat{\sigma}^{2}$ is considerably larger for dendrite 5 than dendrite 1, resulting in more fluctuating retention probabilities.

As discussed in Section 3.2, for the statistical analyses of point patterns on linear networks there is only a limited number of options for model checking, especially when the $K$-function or the pair correlation function have already been used to estimate the model parameters. One simple option is to look at simulations from the fitted model as in Figures C.3–C.8 in Appendix C. Comparing these simulations visually to the observed point patterns, it seems that the simulations mimic the behaviour of larger empty areas without spines seen in the data. For a more rigorous model checking, we performed global rank envelope tests with a concatenation of $\hat{F}$, $\hat{G}$, and $\hat{J}$ as test function, where distances less than $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ were disregarded as discussed in Section 4.1; results are shown in Figure 4. Except from dendrite 4, where the data curve for $\hat{J}$ falls below the global rank envelope for some of the smaller $r$-values, the tests do not reveal any evidence against the fitted models. However, for some of the dendrites (especially dendrite 6) the global rank envelopes for the part concerning $\hat{F}$ and $\hat{G}$ are very wide due to a large variance in the number of points.

5 Discussion

Disregarding the possible small-scale repulsion discussed in Section 4.1, the fitted inhomogeneous Poisson models seemed adequate for dendrites 1 and 3, while the Cox model proposed in Section 4.2 is preferred for the remaining dendrites. To incorporate the small-scale repulsion in the Cox model, the inhomogeneous Poisson process $Y$ used to build the Cox process could be replaced by an inhomogeneous and repulsive point process. The simplest case may be to use a dependent thinning as in a Matérn hard core process of type I (Matérn, , 1960, 1986): let $\tilde{Y}$ be an inhomogeneous Poisson process (with constant intensity on the main branch respectively the side branches), and let

[TABLE]

where $h>0$ is a hard core parameter; that is, a point in $\tilde{Y}$ is included in $Y$ if and only if no other point in $\tilde{Y}$ is within distance $h$. However, it is doubtful whether an expression for the $K$-function or the pair correlation function can be found for such a Cox process model, posing new challenges with respect to parameter estimation.

To avoid using the rather ad hoc created summary functions $\hat{F}$, $\hat{G}$, and $\hat{J}$, the summary functions proposed by Cronie et al., (2019) should be used instead. These were first introduced after the preparation of this paper, but global rank envelope tests based on a concatenation hereof can be found in Figure C.9 in Appendix C for the fitted Cox model. In our case, conclusions based on global rank envelopes using a concatenation of the $\hat{F}$, $\hat{G}$, and $\hat{J}$-functions presented in Section 3.2 are consistent with those using the summary functions from Cronie et al., (2019): except for dendrite 4, the observed summary functions stay strictly inside the 95% global rank envelopes.

Rakshit et al., (2017) discussed the importance of how distance is measured when analysing point patterns on a linear network and they generalised the $K$-function to allow the use of any distance metric. In fact, following Anderes et al., (2020) all methods as well as Poisson and Cox process models in this paper immediately apply for more general linear networks, called graphs with Euclidean edges, when the correlation function $c$ is isotropic with respect to the shortest path distance as well as another metric called the resistance metric. For the dendrite networks or any other tree network, the resistance metric is equivalent to the shortest path distance. Anderes et al., (2020) showed that correlation functions that are isotropic with respect to the shortest path distance only are guaranteed to be valid for a small class of linear networks, whereas they are valid for any linear network when considering the resistance metric instead. Thus, depending on the network, it may be preferable to consider the resistance metric over the shortest path distance when specifying a correlation function. Anderes et al., (2020) provided a list of valid isotropic covariance functions for graphs with Euclidean edges.

Appendix A Expressions for $K$

For the Cox point process presented in Section 4.2 with $c$ equal to the exponential correlation function in (4.4), the $K$-function is

[TABLE]

Let $\alpha=\left(1+1/\sigma^{2}\right)^{-2}$, then $K(r)$ is given by

[TABLE]

if $k=1$,

[TABLE]

if $k=2$,

[TABLE]

if $k=3$,

[TABLE]

if $k=4$, and

[TABLE]

if $k=5$.

Appendix B Simulation study concerning estimation procedure

B.1 Second order composite likelihood

In Section 4.2 we fitted the parameters of the Cox process models using a two step procedure involving a minimum contrast procedure for estimating $\sigma^{2}$ and $\beta$. Another option is to consider a second order composite likelihood approach (adapted from Waagepetersen, , 2007, to a point pattern on a linear network). That is, for an observed point pattern $x\subset L$, the maximum composite likelihood estimate is obtained by maximizing the log composite likelihood

[TABLE]

where $w$ is a weight function and $\neq$ over the summation sign means that $u\neq v$. We can for example let

[TABLE]

for some user specified value $r_{0}$. To estimate $(\sigma^{2},\beta)$, we either directly maximise (B.1) or alternatively solve the associated estimating equation obtained by setting the score equal to 0. The score function is in our set-up given by

[TABLE]

To improve the composite likelihood estimation procedure, Lavancier et al., (2018) suggested an adaptive version of (B.2) where the weight function $w$ depends on the model parameters. We may for example let $w$ be the indicator function given by

[TABLE]

where $M(u,v)=\max_{s\in\{u,v\}}\lvert g(s,s)-1\lvert$ and $\epsilon\in(0,1)$ is a small user-specified number, e.g. $\epsilon=0.01$ or $\epsilon=0.05$. Note that for an isotropic correlation function $g(u,v)=g_{0}\{d_{L}(u,v)\}$, we have $M(u,v)=|g_{0}(0)-1|$. Another weight function suggested in Lavancier et al., (2018) is

[TABLE]

where $h(u,v)=\epsilon M(u,v)/\{g(u,v)-1\}$.

For approximating the double integral in (B.2) (or in the adapted version), note that this is of the form $\int_{L}\int_{L}f(u,v)\,\mathrm{d}_{L}u\,\mathrm{d}_{L}v$. We split up the integration area into the line segments constituting $L$, that is,

[TABLE]

Note that $\hat{\rho}(\cdot)$ for the spine data is constant on any line segment $L_{i}$ (as the line segment is either fully contained in $L_{m}$ or $L_{s}$). Further, if $L$ is a tree, $w$ is given by (B.1), (B.3), or (B.4), and $g$ is isotropic, then $f(u,v)=f_{0}\{d_{L}(u,v)\}$ depends only on distance; this will ease the approximation of the integral:

[TABLE]

where $d_{i,j}=\min_{u\in L_{i},v\in L_{j}}d_{L}(u,v)$. Each of these integrals can then be approximated by Monte Carlo integration using uniform variables on $[0,|L_{i}|]$ and $[0,|L_{j}|]$.

B.2 Simulation study

In the following we describe and summarise results from a simulation study investigating how well $\sigma^{2}$ and $\beta$ are estimated using either the minimum contrast procedure with $T=K$ or $T=g_{0}$ (as described in Section 4.2.2) or the adaptive composite likelihood procedure using (B.3) or (B.4). We considered the network for dendrite 4 and simulated from the Cox process described in Section 4.2 with different parameter values given in Table 3 and when $k=1$ was fixed.

For the minimum contrast procedure we investigated different values of the tuning parameters $r_{u}$, $r_{l}$, and $p$ as well as different start parameters for the optimisation algorithm, see Table 3. Here run no. 1 is the reference run from which one (or two) model or tuning parameters are changed at the time. For run no. 1, we chose $\sigma^{2}=5$ and $\beta=0.1$ resulting in a model rather far away from the case of a Poisson process. Note that decreasing $\beta$ will increase the range of correlation in the thinning probability, whilst increasing $\sigma^{2}$ will increase the probability of thinning. Thus a small $\beta$ and a large $\sigma^{2}$ yield a model very different from the Poisson process. For each choice of model parameters we simulated 500 point patterns and estimated $(\sigma^{2},\beta)$ using minimum contrast and for a few selected runs we also estimated $(\sigma^{2},\beta)$ using the adaptive composite likelihood method.

For the adaptive composite likelihood method the integral in (B.2) was approximated using $10^{6}$ simulations. Estimates of $(\sigma^{2},\beta)$ were found by minimising the length of the score over a $100\times 100$ grid centred around the true values of $\sigma^{2}$ and $\beta$. The finer and broader grid, the better, but as a $100\times 100$ grid was already quite time consuming we settled with that.

B.2.1 Results using the minimum contrast procedure

In the following we let MCE-$K$ and MCE-$g$ refer to minimum contrast estimation with $T=K$ and $T=g_{0}$, respectively.

Histograms of the obtained estimates are shown in Figures B.1–B.4. In general it seems that MCE-$g$ performs better than MCE-$K$. For parameter values that result in models close to the Poisson model, neither of the estimation procedures estimate $(\sigma^{2},\beta)$ successfully. For simulations more distinguishable from Poisson (as for example run no. 1), the MCE-$g$ gives more satisfactory estimates.

Neither MCE-$g$ or MCE-$K$ seem to be sensitive to the choice of start parameters for the reference run. Further, for the reference run it seem that $r_{u}=20$ was the best choice, but in general this depend on the model we are trying to fit and naturally on the size of the network. The best choice of $r_{l}$ seem to be $r_{l}=0$, while $p=0.25$ seems preferable over $p=0.5$ for MCE-$K$, and $p=0.5$ and $p=1$ perform equally well for the MCE-$g$. The choice of $r_{u}$ seem to be important with respect to bias and variance: a too high $r_{u}$ may lead to a large bias, while a too smale $r_{u}$ may lead to a large variance in the estimates. It is therefore recommendable to perform a small simulation study for the specific network and proposed model at hand, such that the best choice of $r_{u}$ can be made.

B.2.2 Results using the adaptive composite likelihood

For the CLE procedure we restricted ourselves to a small number of runs as the grid search was very time consuming. Specifically, we simulated 500 point patterns from the Cox process models with the parameters from run no. 1, 10, and 11 in Table 3. For each choice of model parameters we estimated $(\sigma^{2},\beta)$ using both weight functions and $\epsilon=0.05$. Further, we also ran the CLE procedure with $\epsilon=0.01$ for the model parameters from run no. 1. Figure B.5 shows histograms of the resulting estimates. It is clear that for all three choices of model parameters the grid should be broader in order to find the parameter values that give the smallest length of (B.2) and that the estimates are worse than the ones obtained by the minimum contrast procedures. Finally, there is no seemingly advantage of choosing one weight function over the other or of choosing $\epsilon=0.01$ over $\epsilon=0.05$.

Appendix C Analysis of spine locations

This Appendix contains figures related to the analysis of the six spine data sets. Figure C.1 shows $95\%$ global rank envelopes under the fitted inhomogeneous Poisson model using a concatenation of $\hat{F}$, $\hat{G}$, and $\hat{J}$ as test function, where distances less than $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ are disregarded. Figure C.2 show one simulation of the fitted random field $\Pi$ for each network, and each of Figures C.3–C.8 display the data along with five simulated point patterns from the fitted Cox process model. Finally, Figure C.9 shows $95\%$ global rank envelopes under the fitted Cox model using a concatenation of the estimated $F$-, $G$-, and $J$-functions proposed by Cronie et al., (2019) as test function, where distances less than $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ are disregarded.

Acknowledgements

We would like to thank Abdel-Rahman Al-Absi who collected the spine data. We have been supported by The Danish Council for Independent Research | Natural Sciences, grant DFF – 7014-00074 “Statistics for point processes in space and beyond”, and by the “Centre for Stochastic Geometry and Advanced Bioimaging”, funded by grant 8721 from the Villum Foundation.