Sparse Reduced-rank Regression Methods for Spatially Misaligned Data with Application to Spatial Transcriptomics

Abstract

Understanding the spatiotemporal dynamics of disease progression in relation to transcriptomic profiles provides key insights into complex conditions such as Alzheimer’s disease. To enable such investigations, STARmap PLUS technology offers joint profiling of high-resolution spatial transcriptomics and protein detection within the same tissue section. Motivated by data from Zeng et al. (2023), we develop a novel kernel-weighted regression framework that models plaque size as a collective effect of the spatial transcriptomics of neighboring cells, automatically integrating across cell types and tissue samples from different disease states. To further strengthen interpretability and efficiency, we incorporate a sparse low-rank factorization that enables gene selection while borrowing strength across genes, cell types, and time points. The proposed approach is implemented in a fully automated manner with data-driven specification of key model components. Through simulation studies, we demonstrate the robustness of the proposed method and its superiority across a range of specification scenarios. Applied to Alzheimer’s data, the proposed framework uncovers biologically meaningful associations, highlighting its potential for advancing the understanding of disease mechanisms.

Keywords: Cellular microenvironment; Kernel-weighted regression; Sparse Low-rank factorization; STARmap PLUS

1 Introduction

Alzheimer’s disease (AD) comprises interacting pathologies that vary across space, cell types, and time (De Strooper and Karran, 2016; Jack Jr et al., 2018, 2024). Amyloid-$\beta$ plaques and hyperphosphorylated tau accumulate within a dynamic microenvironment shaped by microglia, astrocytes, oligodendrocyte-lineage cells, and neurons (Long and Holtzman, 2019; Hampel et al., 2021). Spatially resolved studies in mouse and human reveal reproducible plaque-adjacent niches and layer-specific molecular patterns, linking microglial/astrocytic remodeling and tau-aligned neuronal changes to local pathology (Chen et al., 2020, 2022; Zeng et al., 2023; Mallach et al., 2024; Miyoshi et al., 2024). Motivated by these spatially organized niches, we adopt a lesion-centric formulation: we treat plaque centers as outcome locations and quantify how nearby cellular gene expression relates to the plaque size (mean radius).

Traditional transcriptomic designs struggle to capture this heterogeneity. Bulk RNA-seq averages across mixed cell types and loses cellular heterogeneity and the spatial context of the cells, while single-cell RNA-seq restores cellular resolution but still discards in-tissue organization. Spatially resolved transcriptomics (SRT) maps gene expression to spatial coordinates via two broad classes: (i) barcoded capture arrays (Spatial Transcriptomics/Visium) that provide transcriptome‐wide profiles at multi-cell spot resolution (Ståhl et al., 2016; Vickovic et al., 2019; Maynard et al., 2021); and (ii) imaging-based assays (e.g., MERFISH(Chen et al., 2015), seqFISH(Shah et al., 2016), STARmap(Wang et al., 2018)) that achieve single-cell or subcellular resolution with targeted gene panels.

Among imaging-based SRT, STARmap PLUS (Zeng et al., 2023) co-registers targeted RNAs with amyloid-$\beta$ and p-tau immunostains within the same tissue sample, yielding single-cell expression aligned to plaque segmentations and centers, and enabling direct linkage between molecular readouts and neuropathology in situ. Each tissue section contains many plaques and thousands of cells with known coordinates; plaque size varies markedly within a section, and this variability is embedded in distinct cellular neighborhoods defined by which cells lie nearby and what genes they express. This provides exactly the data structure needed for a plaque-centric analysis.

Beyond Alzheimer’s disease, co-registering histopathology with SRT enables tumor-region delineation, quantification of tissue architecture and burden, and pathology-informed disease modeling (Ni et al., 2022; Arora et al., 2023; Bassiouni et al., 2023; Xun et al., 2023; Jin et al., 2024). Related statistical questions also arise in environmental applications (Alexeeff et al., 2016). This breadth underscores the generality and translational relevance of our framework for modeling disease-related spatial outcomes from spatial molecular measurements.

The primary statistical challenges are: (1) spatial misalignment between plaque locations and SRT measurement locations, (2) the high dimensionality of SRT data, (3) effect heterogeneity across collection times and cell types, i.e., regression coefficients linking plaque size to neighborhood features vary by time point and cell type, with gene-level coefficients exhibiting shared latent structure, (4) the need for a spatially structured characterization of the collective influence of different neighboring cells on plaque size, and (5) between-sample heterogeneity in plaque-adjacent microenvironments (e.g., shifts in local cell-type composition and activation states).

For regressing spatially varying outcomes on predictors from nearby locations, spatial measurement error models have been proposed, such as two-stage correction approaches (Szpiro and Paciorek, 2013), and spatial simulation extrapolation methods (Alexeeff et al., 2016). Local/kernel regressions and geographically weighted regression (GWR) model local heterogeneity by weighting nearby observations more heavily; multiscale GWR further allows each predictor have its own spatial scale parameter. However, they typically treat constructed predictors as error-free and do not scale or pool information effectively in high-dimensional, structured SRT settings (Nadaraya, 1964; Fotheringham et al., 2017; Fan, 2018; Oshan et al., 2019). Change-of-support and downscaler models propagate aggregation/registration uncertainty, yet classical formulations remain low-dimensional and do not accommodate high-dimensional features or latent commonalities across groups (Gotway and Young, 2002; Berrocal et al., 2010). Together, these approaches address parts of the problem but do not provide a unified, high-dimensional regression that learns predictor-specific spatial scales while quantifying the collective impact of neighboring cells. This gap motivates our framework.

In an asynchronous longitudinal regression setting, Cao et al. (2015); Li et al. (2022, 2023) proposed a kernel-weighted approach to handle mismatched observation times between the response and covariate within a given subject. We follow a similar direction by introducing a kernel function to borrow information from all possible pairs of outcomes and the predictors in the misaligned spatial setting. This particularly helps to characterize the collective effect of different cells from varying cell types on the plaque size. An isotropic proximity kernel with sample-specific bandwidths aggregates local signals, accommodating section-to-section enrichment and sampling differences without overfitting directionality. Gene effects are parameterized by a low-rank coefficient tensor over gene $\times$ cell type $\times$ time, together with gene-mode regularization to capture shared structure and heterogeneity across samples and time points under high dimensionality. To achieve interpretable gene-level discovery, we apply a LASSO penalty only on the gene mode of this factorization: genes with weak signal are shrunk toward zero, while the associated cell-type and time patterns for the retained genes remain flexible. All tuning parameters (e.g. kernel bandwidths, rank, and penalty levels) are selected in data-driven manner; details are provided in Section 3.1. The resulting framework yields interpretable component-level profiles and signed effects, where the sign of the estimated regression coefficient indicates direction (positive = larger plaque size; negative = smaller) across cell types and time points.

Applied to STARmap PLUS AD sections, the model identifies 13-month multigene expression patterns dominated by a Microglia–Oligodendrocyte axis: a microglial Apoe/Tyrobp module and an oligodendrocyte-lineage (OPC-like) pattern are positively associated with plaque size, whereas an astrocytic transport-related signature is negatively associated with plaque size (see Section 5).

We first describe the STARmap PLUS data and exploratory analyses motivating a plaque-centric estimand and isotropic neighborhood kernel. We then develop the kernel-weighted regression with sample-specific bandwidths and low-rank, gene-mode selection), assess performance in spatially realistic simulations, analyze AD sections, and close with limitations and extensions.

2 Data and Exploratory Analysis

2.1 Data

In this study, we use the single-cell spatial transcriptomics data integrated with amyloid-$\beta$ plaque histopathology from a mouse model of Alzheimer’s disease (AD). Specifically, data were obtained from the brains of TauPS2APP transgenic mice at two critical stages of disease progression: 8 months (8 mo) and 13 months (13 mo) of age, providing a temporal snapshot of AD pathology and cellular heterogeneity. Each of the four samples analyzed includes comprehensive gene expression data at single-cell resolution from all identified cells, alongside spatial localization of extracellular amyloid-$\beta$ plaques.

The SRT data employed the STARmap PLUS method (Zeng et al., 2023), which simultaneously captures spatially resolved single-cell expression and protein-level histopathological markers within the same tissue section at subcellular resolution. In this study, STARmap PLUS profiled 2,766 genes while simultaneously labeling extracellular amyloid-$\beta$ plaques (X-34 dye) and hyperphosphorylated tau (AT8 immunostaining), enabling precise co-registration of cell bodies, mRNA, and pathology. Although both amyloid-$\beta$ and tau signals are available, the present study focuses exclusively on amyloid-$\beta$ plaques. Voxel size was $95\times 95\times 350$ nm with $56-64$ fields of view per sample. The method preserves morphology and allows joint RNA-protein mapping at subcellular scales, unlike chip-based whole-transcriptome approaches with coarser pixels ($100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ pixel size) (Chen et al., 2020).

In the original dataset, STARmap PLUS yielded a high-resolution spatial atlas across cortex and hippocampus and identified $13$ major cell types over $72$k cells. Our analysis uses the diseased subset to focus specifically on plaque-proximal biology. Two features make this dataset uniquely informative for plaque-centric analyses: (i) exact spatial location of cells, gene expression, and plaque pathology within the same section, and (ii) clear distance-dependent patterns in cellular composition around plaques.

Figure 1: Co-registered spatial maps from two distinct tissue sections collected at two time points: (a) 8 months and (b) 13 months. Cells are shown as gray points, and amyloid plaques are overlaid as red filled circles with radii encoding plaque sizes. Gene expression is measured for each cell at its observed spatial coordinate.

From Zeng et al. (2023), we obtain the plaque size data. In that, they first identified the plaques’ centers and then computed their mean radii. Fig 1 shows co-registered spatial maps of two mouse brain tissue sections collected at (a) 8 months and (b) 13 months. Gray points denote the spatial locations of all cells, and amyloid plaques are shown as red filled circles with circle radius equal to plaque sizes in $\mathrm{\SIUnitSymbolMicro m}$.

2.2 Exploring the Spatial Variability

To characterize plaque-adjacent microenvironments, we re-examined plaque-centric organization using concentric rings around plaques (0–30, 30–70, 70–150 µm) in Figs. 2 and 3. Motivated by the concentric-ring analyses of Zeng et al. (2023), who defined distances with respect to plaque boundaries, we use the Euclidean distances from the cell centroids to the nearest plaque center ${\mathbf{U}}$ in our analysis to define the plaque microenvironment.

We focus on three cell-types: Microglia, Astrocytes, and Oligodendrocyte-lineage cells, because both our exploratory analyses (Figs. 2 and 3) and prior work by Zeng et al. (2023) revealed robust, plaque-centric organization for these types with sufficient counts for stable stratification. Microglia are enriched right at plaques center with “core-shell” (as defined by Zeng et al. (2023) concerning the nearest concentric-ring) architecture closely contacting amyloid-$\beta$ plaques ([range-phrase = –, range-units = single]030 from the plaque center), while Astrocytes and Oligodendrocytes are enriched in surrounding shells ([range-phrase = –, range-units = single]30150 from the plaque center). Across sections and between 8 and 13 months, enrichment magnitudes generally increase, while beyond 150 µm no systematic enrichment is detected.

Figure 2: Plaque-centric distance maps. Shaded annuli show Euclidean distance bands from the nearest plaque center: 0–30, 30–70, 70–150, and $\geq$ 150 µm. Points are cell locations colored by three cell types (legend). Each cell is placed in the ring at which it falls; these maps define the radial stratification used in downstream summaries. The two tissue sections are slightly different in dimension.

Figure 3: Cell density by distance to the nearest plaque. Bars show mean cell density (cells/mm²) within each radial band (0–30, 30–70, 70–150, $\geq$$150\text{\,}\mathrm{\SIUnitSymbolMicro m}$) and overall, for Microglia, Astrocyte, and Oligodendrocyte-lineage cells. Densities are computed per sample and averaged across samples. Enrichment near plaques is most pronounced for Microglia in [range-phrase = –, range-units = single]030 and for Astrocyte/Oligodendrocyte-lineage cells in [range-phrase = –, range-units = single]30150.

Guided by these spatial cell enrichments, we next construct bulk-like plaque-level summaries of near-plaque gene expression stratified by plaque size for descriptive visualization, to assess whether near-plaque expression varies with plaque size in Fig 4. The nearest plaque for each cell is identified by sorting Euclidean distances from the cell centroid to different plaque centers. A cell is considered near-plaque if its distance to the nearest plaque center is $\leq$ 150 µm. For visualization only, plaque size is dichotomized based on the global median radius, namely Small: plaques with radii $\leq$ median; Large: plaques with radii $>$ median. For each gene, heatmap values in Fig 4 are computed via two-stage averaging: first, we compute Plaque-level (bulk-like) mean for each plaque as the average gene expression over all near-plaque cells, yielding one value per (plaque, gene); and subsequently, we compute the average of plaque-level means for each time (8 mo vs 13 mo), plaque-size (Small vs Large). This finally gets us one value per time, size, and gene, which forms the heatmap entries.

Fig 4 tells us: (i) the difference-in-differences for each gene expression, showing how the 8$\rightarrow$13 mo change near Large plaques compares with the change near Small plaques (positive bars = larger increase or smaller decrease near Large plaques), (ii) a heatmap of the per-size change from 8 mo to 13 mo shown separately for Small and Large, and (iii-iv) heatmaps of expression at 8 mo and 13 mo (columns = Small vs Large). Row centering (within time) is applied only to the expression heatmaps (iii–iv); the bar plot and the change heatmap show raw contrasts. Across the 30 most varying genes over time, these summaries reveal size-associated differences in near-plaque expression and motivate our regression problem below, where plaque radius is regressed on the local gene expression levels.

Figure 4: Near-plaque gene expression by plaque size. Left: bar plot of $\Delta_{\text{Large}}-\Delta_{\text{Small}}$ where $\Delta$ denotes the 13 mo minus 8 mo change in near-plaque expression within each plaque-size group (positive indicates a larger temporal increase near Large plaques). Middle-left: heatmap of $\Delta$ shown separately for Small and Large plaques. Middle-right / right: heatmaps of near-plaque expression at 8 months and 13 months, respectively (columns: Small vs Large). Rows are the 30 genes with the largest temporal variation. Row centering (within time) is applied only to the expression heatmaps; the bar plot and change heatmap show raw contrasts.

3 Model

We propose a novel kernel-weighted regression framework to analyze the effect of spatial transcriptomic profiles on amyloid-$\beta$ plaque deposition. For clarity, we first introduce some notations that are used throughout the article.

Let $n$ stands for the total number of tissue samples. For each sample $i=1,\dots,n$, let $M_{i}$ denotes the number of outcome locations (e.g., plaques) and $\{{\mathbf{U}}_{i,j}\}_{j=1}^{M_{i}}$ be their spatial coordinates. At each of these locations, we observe spatially indexed scalar outcomes, such as plaque sizes, denoted by $\{y_{i,j}\coloneqq y({\mathbf{U}}_{i,j}),\,j=1,\dots,M_{i}\}$. Let $N_{i}$ stands for the number of cells in $i$-$th$ sample and $\{{\mathbf{V}}_{i,k}\}_{k=1}^{N_{i}}$ be their spatial coordinates. The spatial transcriptomics data, recording gene expression levels of $p$ genes, are available at individual cell locations, denoted by $\{{\bm{x}}_{i,k}\coloneqq{\bm{x}}({\mathbf{V}}_{i,k})\in\mathbb{R}^{p},\,k=1,\dots,N_{i}\}$. Furthermore, the samples are obtained at different stages of the disease progression, where $t_{i}\coloneqq t(i)$ represents the corresponding collection time of the $i$-$th$ sample. Note that all tissue samples are independent, as they come from different mice. Each cell further belongs to one of the $C$ cell-types. For clarity, let $c_{i,k}\coloneqq c({\mathbf{V}}_{i,k})$ represent the cell type associated with the cell at location ${\mathbf{V}}_{i,k}$ in $i$-$th$ sample.

Our goal is to estimate the effect of gene expression levels on plaque sizes, measured at different sample collection times. Specifically, we aim to estimate $\mbox{$\beta$}_{c,t}\in\mathbb{R}^{p}$, whose $j$-$th$ value quantifies the effect of the $j$-$th$ gene’s expression (among the $p$ genes) in the $c$-$th$ cell type on plaque sizes collected at the $t$-$th$ time point. However, to run this analysis, we must take into account the following three characteristics in our proposed model:

1) The spatial locations of the cells and plaques do not match in general.

2) Each plaque is surrounded by cells with different cell types.

3) The gene-specific regression effects for different cell types may share some commonalities.

4) Samples collected at different time points are independent, but gene-specific effects may share commonalities through underlying latent factors that capture shared biological or regulatory dependencies.

In order to resolve these issues, we utilize the inherent spatial dependence in the data. Specifically, aligning with the realistic scenarios arising from integrating single-cell spatial transcriptomics data with amyloid-beta plaque data in Alzheimer’s disease, spatial proximity informs similarities in gene expression levels and biological response patterns, using sample-specific, density-adaptive neighborhoods.

Thus, to characterize the relationship between spatially mismatched $\bm{x}_{i,k}$ and $y_{i,j}$, we first propose the following kernel-weighted objective function,

$$\sum_{i=1}^{n}\sum_{k=1}^{N_{i}}\sum_{j=1}^{M_{i}}K_{i,j,k}\cdot\left(y_{i,j}-{\bm{x}}_{i,k}^{\top}\mbox{$\beta$}_{c_{i,k},t_{i}}\right)^{2},$$ (1)

where $c_{i,k}\in\{1,\dots,C\}$ standing for the cell-type of the cell present at location ${\mathbf{V}}_{i,k}$, $t_{i}\in\{1,\dots,T\}$ is the collection time for sample $i$, and $K_{i,j,k}\coloneqq K_{h_{i}}\left(\|{\mathbf{U}}_{i,j}-{\mathbf{V}}_{i,k}\|_{2}\right)$ are nonnegative kernel weights with possibly sample-specific bandwidth $h_{i}$. For the kernel $K_{h}(\cdot)$, we consider the Epanechnikov kernel to assign weights based on the Euclidean distance between ${\mathbf{U}}_{i,j}$ and ${\mathbf{V}}_{i,k}$. This compactly supported kernel determines the contribution of each observation to the local fit. The Epanechnikov kernel is a widely used choice due to its optimality in minimizing mean squared error. The unscaled kernel is defined as $K(d)=0.75(1-d^{2})_{+}$, where $(\cdot)_{+}$ denotes the positive part. Given a bandwidth parameter $h>0$, the kernel is scaled as $K_{h}(d)=K(d/h)/h$. This formulation ensures that observations farther than $h$ units apart receive zero weight, while closer observations have smoothly decreasing influence as their distance approaches $h$. Now, there are $pGT$ unknown parameters to be estimated in $\{\mbox{$\beta$}_{c,t}\in\mathbb{R}^{p}\}_{c=1}^{C}{}_{t=1}^{T}$. Direct estimation of these parameters is computationally infeasible and fails to exploit the inherent dependencies among genes, cell types, and time. Moreover, sparsity-inducing structures are required to identify important genes. Thus, we propose a sparse reduced-rank factorization of the regression coefficient to identify important genes while borrowing information across different cell types and time.

Low-rank factorization of the regression coefficient:

In order to judiciously estimate the regression coefficient $\bm{\beta}_{c,t}$, we utilize a reduced-rank factorization that leverages shared structure across cell types, time, and genes. In that, we first rewrite it as a tensor $\mathcal{B}\in\mathbb{R}^{p\times C\times T}$ with three dimensions: Genes $\times$ Cell-types $\times$ Times, where the $(l,c,t)$-$th$ entry is $\beta(l,c,t)$, for $l=1,\dots,p$, $c=1,\dots,C$, and $t=1\dots,T$. One may then apply the tensor-factorization methods, such as CANDECOMP/PARAFAC (CP) Carroll and Chang (1970); Harshman (1970); Kiers (2000), or more generally the Tucker decomposition Tucker (1966). In this work, we employ a CP decomposition on $\mathcal{B}$ as:

$$\beta(l,c,t) =\sum_{r=1}^{R}w_{r}q_{l,r}^{(1)}q_{c,r}^{(2)}q_{t,r}^{(3)}$$

$$\mathcal{B} =\sum_{r=1}^{R}w_{r}{\bm{q}}^{(1)}_{r}\circ{\bm{q}}^{(2)}_{r}\circ{\bm{q}}^{(3)}_{r}$$

$${\bm{q}}^{(1)}_{r} \coloneqq\left(q_{l,r}^{(1)}\right)_{l=1}^{p},{\bm{q}}^{(2)}_{r}\coloneqq\left(q_{c,r}^{(2)}\right)_{c=1}^{C},{\bm{q}}^{(3)}_{r}\coloneqq\left(q_{t,r}^{(3)}\right)_{t=1}^{T},$$

$$\text{such that }w_{r}>0,\|{\bm{q}}^{(m)}_{r}\|_{2}=1,\text{for all }r=1,\ldots,R;m=1,2,3,$$

where $\|\cdot\|_{2}$ stands for the $\ell_{2}$ norm which is defined as $\|\mathbf{a}\|=\sqrt{\sum_{i}a^{2}_{i}}$. By restricting ${\bm{q}}^{(1)}_{r}$, ${\bm{q}}^{(2)}_{r}$, and ${\bm{q}}^{(3)}_{r}$, to be normalized individually, it ensures that $w_{r}$ captures the scale for each rank. This further helps in rank selection and introducing sparsity, as discussed later. The CP decomposition reduces parameter dimensionality while enabling information sharing across genes, cell types, and time. Specifically, we represent $\{\mbox{$\beta$}_{c,t}\}$ as a sum of $R$ components; $r=1,\ldots,R$ indexes the component (“rank”). Note that under the CP decomposition, the same rank $R$ appears across all modes by construction, since $R$ denotes the number of rank-one components. For greater flexibility, Tucker factorization with mode-specific ranks may be considered. However, for sparse tensors, CP decompositions are usually adequate.

In the gene direction, we place a LASSO penalty on $\bm{q}^{(1)}_{r}=\{q^{(1)}_{l,r}\}_{l=1}^{p}$ for each rank $r$, inducing sparsity in $\mathcal{B}$. The resulting objective function is:

$$F(\mathcal{B})=\frac{1}{2}\sum_{i=1}^{n}\sum_{k=1}^{N_{i}}\sum_{j=1}^{M_{i}}K_{i,j,k}\cdot\left(y_{i,j}-{\bm{x}}_{i,k}^{\top}\mbox{$\beta$}_{c_{i,k},t_{i}}\right)^{2}+\frac{\lambda}{Rp}\sum_{r=1}^{R}\|{\bm{q}}_{r}^{(1)}\|_{1},$$ (2)

where $\|\cdot\|_{1}$ stands for the $\ell_{1}$ norm which is defined as $\|{\bm{a}}\|_{1}=\sum_{i}|a_{i}|$. It allows us to obtain a sparse solution and select the important genes.

3.1 Estimation

The overall implementation is fully automated with data-driven efficient tuning. We update all the parameters using a blocked coordinate descent algorithm. To do that, we first partition the full parameter set into $R$ blocks with $r$-$th$ block: $\{{\bm{q}}^{(1)}_{r},{\bm{q}}^{(2)}_{r},{\bm{q}}^{(3)}_{r},w_{r}\}$. In each iteration, we cycle through $r=1,\dots,R$, and update the components sequentially in each block. The blocked coordinate descent updates are given below. After each update, we normalize ${\bm{q}}^{(k)}_{r}$ by its corresponding $\ell_{2}$ norm and adjust $w_{r}$ accordingly.

We minimize Eq.(2), where ${\bm{x}}_{i,k}=(x_{1,ik},\dots,x_{p,ik})^{\top}$. Write the CP form

$$\mbox{$\beta$}_{c,t}=\sum_{\begin{subarray}{c}s=1\end{subarray}}^{R}w_{s}\,{\bm{q}}^{(1)}_{\cdot s}\,q^{(2)}_{c,s}\,q^{(3)}_{t,s},\qquad\mbox{$\beta$}_{c,t}^{(-r)}=\mbox{$\beta$}_{c,t}-w_{r}\,{\bm{q}}^{(1)}_{\cdot r}\,q^{(2)}_{c,r}\,q^{(3)}_{t,r}.$$

For a fixed component $r$ define

$$z_{ik,r}^{(1)}\coloneqq w_{r}q^{(2)}_{c_{i,k},r}q^{(3)}_{t_{i},r},\qquad s_{ik,r}\coloneqq{\bm{x}}_{i,k}^{\top}{\bm{q}}_{\cdot r}^{(1)},\qquad R_{i,j,k}^{(-r)}\coloneqq y_{i,j}-{\bm{x}}_{i,k}^{\top}\mbox{$\beta$}_{c_{i,k},t_{i}}^{(-r)}.$$

We update $q^{(1)}_{l,r}$ by holding all others fixed, and use the leave-one-coordinate-out residual

$$\tilde{R}_{i,j,k}^{(l,r)}\coloneqq R_{i,j,k}^{(-r)}-z_{ik,r}^{(1)}\left(s_{ik,r}-x_{l,ik}\cdot q^{(1)}_{l,r}\right).$$

Define the sufficient statistics

$$a_{l,r}\coloneqq\sum_{i,k,j}K_{i,j,k}\cdot z_{ik,r}^{(1)}\cdot x_{l,ik}\cdot\tilde{R}_{i,j,k}^{(l,r)},\qquad b_{l,r}\coloneqq\sum_{i,k,j}K_{i,j,k}\left(z_{ik,r}^{(1)}\cdot x_{l,ik}\right)^{2}.$$

Let $S_{\tau}(u)=\mathrm{sign}(u)(|u|-\tau)_{+}$ be the soft-$th$resholding operator where $\tau=\lambda/(Rp)$. Then the closed-from update is

$$q_{l,r}^{(1)}\leftarrow\frac{S_{\lambda/(Rp)}(a_{l,r})}{b_{l,r}},$$ (3)

with the convention $q_{l,r}^{(1)}\leftarrow 0$ if $b_{l,r}=0$. For fixed $(l,r)$, the objective function reduces to $\frac{1}{2}b_{l,r}q^{2}-a_{l,r}q+\frac{\lambda}{Rp}|q|$, whose minimizer is the soft-$th$resholding operator above.

For updating the entry $q_{c,r}^{(2)}$ in the cell-type loadings, we set $z_{ik,r}^{(2)}\coloneqq w_{r}q_{t_{i},r}^{(3)}s_{ik,r}$. This block is unpenalized weighted least squares:

$$q_{c,r}^{(2)}\leftarrow\dfrac{\displaystyle\sum_{i,k:\,c_{i,k}=c}\sum_{j=1}^{M_{i}}K_{i,j,k}\cdot z_{ik,r}^{(2)}\cdot R_{i,j,k}^{(-r)}}{\displaystyle\sum_{i,k:\,c_{i,k}=c}\sum_{j=1}^{M_{i}}K_{i,j,k}\cdot\left(z_{ik,r}^{(2)}\right)^{2}}.$$ (4)

Similarly for the time loadings, we set $z_{ik,r}^{(3)}\coloneqq w_{r}q_{c_{i,k},r}^{(2)}s_{ik,r}$. Then

$$q_{t,r}^{(3)}\leftarrow\dfrac{\displaystyle\sum_{i:\,t_{i}=t}\sum_{k=1}^{N_{i}}\sum_{j=1}^{M_{i}}K_{i,j,k}\cdot z_{ik,r}^{(3)}\cdot R_{i,j,k}^{(-r)}}{\displaystyle\sum_{i:\,t_{i}=t}\sum_{k=1}^{N_{i}}\sum_{j=1}^{M_{i}}K_{i,j,k}\cdot\left(z_{ik,r}^{(3)}\right)^{2}}.$$ (5)

To update entry $w_{r}\geq 0$, we set $z_{ik,r}^{(2)}\coloneqq q_{c_{i,k},r}^{(2)}q_{t_{i},r}^{(3)}s_{ik,r}$. Then

$$w_{r}\leftarrow\max\left\{0,\cdot\dfrac{\displaystyle\sum_{i,k,j}K_{i,j,k}\cdot z_{ik,r}^{(w)}\cdot R_{i,j,k}^{(-r)}}{\displaystyle\sum_{i,k,j}K_{i,j,k}\cdot\left(z_{ik,r}^{(w)}\right)^{2}}\right\}.$$ (6)

Computational efficiency and stability:

In our implementation, we keep track of the residuals $R_{i,j,k}=y_{i,j}-{\bm{x}}_{i,k}^{\top}\mbox{$\beta$}_{c_{i,k},t_{i}}$ and update them incrementally. When a coordinate $q_{l,r}^{(1)}$ changes by $\Delta$, set $R_{i,j,k}\leftarrow R_{i,j,k}-z_{ik,r}^{(1)}x_{l,ik}\Delta$; analogous updates are used for $q_{c,r}^{(2)}$, $q_{t,r}^{(3)}$ and $w_{r}$. This reduces the per-coordinate cost from $O(p\sum_{i}M_{i}N_{i})$ for recomputation to $O(\sum_{i}M_{i}N_{i})$, yielding substantial speed-ups without any approximation.

After updating the factors for each rank, we rescale them to unit $\ell_{2}$ norm and absorb the entire scale into the scaler $w_{r}$. Let $\alpha_{m}=\big\|{\bm{q}}^{(m)}_{\cdot r}\big\|_{2}$ for $m\in\{1,2,3\}$. We set

$${\bm{q}}^{(m)}_{\cdot r}\leftarrow\frac{{\bm{q}}^{(m)}_{\cdot r}}{\alpha_{m}}\quad\text{for }m=1,2,3,\qquad w_{r}\leftarrow w_{r}\,\alpha_{1}\alpha_{2}\alpha_{3},$$

so that the contribution from the $r$-$th$ rank, $w_{r}\,{\bm{q}}^{(1)}_{\cdot r}\,q^{(2)}_{c,r}\,q^{(3)}_{t,r}$, (and thus fitted values) is unchanged. This prevents factor columns from arbitrarily exploding/vanishing with compensating shrink/expansion in $w_{r}$, keeps the coordinate-descent curvature well conditioned, and makes the $\ell_{1}$ penalty on ${\bm{q}}^{(1)}$ act consistently rather than depending on arbitrary scaling of ${\bm{q}}^{(1)}$.

Stopping rules:

One may continue the above iterative process up to a fixed number of iterations. Instead, we continue until the following predefined convergence criterion is met:

$$1) \max\limits_{c,t}\dfrac{\|\mbox{$\beta$}^{(new)}_{c,t}-\mbox{$\beta$}^{(old)}_{c,t}\|_{\infty}^{2}}{\|\mbox{$\beta$}_{c,t}^{(old)}\|_{\infty}^{2}}<10^{-6}.$$

$$2) \max_{r,m}\left\{\|{\bm{q}}^{(m)}_{r,new}-{\bm{q}}^{(m)}_{r,old}\|\right\}\leq 10^{-4}.$$ (7)

The second criterion in Eq.(7) was also used in Anandkumar et al. (2015) and Sun et al. (2017) in the context of tensor-valued response regression.

Automatic CP-rank update:

In our approach, the procedure of selecting CP rank is seamlessly integrated into the blocked Coordinate Descent algorithm. We do not need to specify or tune the rank during model fitting manually. Instead, the method begins by initializing the model with a sufficiently large CP rank, which we recommend setting to the product of the number of cell type ($C$) and time points ($T$), i.e., $C\times T$. The general upper bound on the maximum CP rank for any three-way tensor $\bm{\mathcal{X}}\in\mathbb{R}^{I\times J\times K}$ satisfies (Kruskal, 1989; Kolda and Bader, 2009)

$$\mathrm{rank}(\bm{\mathcal{X}})\leq\min\{IJ,IK,JK\}.$$ (8)

A rank-dropping strategy is then applied during the first 500 iterations, and the algorithm actively prunes unnecessary components in these steps. This strategy ensures that the subsequent updates focus on stabilizing the estimation across iterations based on the selected rank after the rank-dropping. Specifically, the CP rank of $\mathcal{B}$ is reduced from $R$ to $R-1$ by discarding $r$-$th$ component if any of the following conditions are met:

$$1)\frac{w_{r}}{\sum_{i}w_{i}}<10^{-5},\qquad 2)w_{r}=0,\qquad 3)\|{\bm{q}}_{r}^{(1)}\|_{\infty}>0.99.$$ (9)

By automatically discarding redundant or negligibly small components, the rank-dropping strategy allows the model to achieve both good estimation accuracy and computational efficiency by selecting a final rank that balances model complexity and goodness-of-fit automatically.

Lambda selection:

We first fit the method for different choices of $\lambda$ and then select the optimal $\lambda$ based on a BIC-like criteria. To ensure computational efficiency, we begin with the penalty parameter $\lambda$ that yields the all-zero solution and progressively decrease it, as done in the path solution algorithms (Efron et al., 2004; Friedman et al., 2010). The successive penalty parameters are then set recursively as $\lambda_{k+1}=0.9\lambda_{k}$. Finally, for each penalty, we evaluate the following BIC-like criterion to do the final selection:

$$N^{*}\log\left[\frac{1}{N^{*}}\sum_{i=1}^{n}\sum_{k=1}^{N_{i}}\sum_{j=1}^{M_{i}}K_{i,j,k}\cdot\left(y_{i,j}-{\bm{x}}_{i,k}^{\top}\mbox{$\beta$}_{c_{i,k},t_{i}}\right)^{2}\right]+\nu\log(N^{*}),$$ (10)

where $N^{*}=\sum_{i=1}^{n}\sum_{k=1}^{N_{i}}\sum_{j=1}^{M_{i}}\mathbf{1}\{K_{i,j,k}>0\}$ and $\nu$ is the number of nonzero entries in ${\bm{q}}^{(1)}_{r}$, ${\bm{q}}^{(2)}_{r}$, and ${\bm{q}}^{(3)}_{r}$. Although one can choose a smaller decay factor than 0.9 for a finer exploration, in our experiments, this choice works reasonably well.

Bandwidth selection:

The kernel bandwidth $h_{i}$ is set based on an elbow rule, as discussed below. First, we set the candidate bandwidths as the medians of the distances from each plaque to its $L$-$th$ closest neighboring cell for different choices of $L$’s. Let the median distance based on $L$-$th$ closest neighboring cell in $i$-$th$ sample be $H_{i}(L)=\text{Median}(\{e^{(L)}_{i,j}\}_{j=1}^{N_{i}})$, where $e^{(L)}_{i,j}=D_{(L)}({\mathbf{U}}_{i,j})$ where $D_{(1)}({\mathbf{U}}_{i,j})\leq D_{(2)}({\mathbf{U}}_{i,j})\leq\cdots\leq D_{(M_{i})}({\mathbf{U}}_{i,j})$ stands for the sorted sequence of the distances of different cells from $j$-$th$ plaque i.e. $\{\|{\mathbf{U}}_{i,j}-{\mathbf{V}}_{i,k}\|_{2}\}_{k=1}^{M_{i}}$.

Finally, we plot the normalized loss

$$\dfrac{1}{N^{*}}\sum_{i=1}^{n}\sum_{k=1}^{N_{i}}\sum_{j=1}^{M_{i}}K_{H_{i}(L)}(\|{\mathbf{U}}_{i,j}-{\mathbf{V}}_{i,k}\|_{2})\left(y_{i,j}-{\bm{x}}_{i,k}^{\top}\mbox{$\beta$}_{c_{i,k},t_{i}}\right)^{2}$$

as a function of the neighborhood size $L$ and select $L$ using the elbow (point-of-diminishing-returns) criterion, defined as the value of $L$ maximizing the perpendicular distance to the line connecting the endpoints of the loss curve (Satopaa et al., 2011). Even with the same $L$ for all samples, the sample-specific bandwidths $(H_{i}(L))$ may be different.

3.2 Initialization

To initialize the model parameters $w_{r}$, ${\bm{q}}^{(1)}_{r}$, ${\bm{q}}^{(2)}_{r}$, and ${\bm{q}}^{(3)}_{r}$, we apply a weighted Ridge regression without the CP decomposition structure. For each pair $(c,t)$, define the objective

$$\mathcal{L}_{c,t}(\mbox{$\beta$})=\sum_{i:t_{i}=t}\sum_{k:c_{i,k}=c}\sum_{j=1}^{M_{i}}K_{i,j,k}\cdot\left(y_{i,j}-{\bm{x}}_{i,k}^{\top}\mbox{$\beta$}\right)^{2}+\lambda_{c,t}\|\mbox{$\beta$}\|_{2}^{2},$$

and the estimator

$$\hat{}\mbox{$\beta$}_{c,t}^{(\text{Ridge})}\in\arg\min_{\mbox{$\beta$}\in\mathbb{R}^{p}}\mathcal{L}_{c,t}(\mbox{$\beta$}).$$

Here, we put the $l_{2}$ penalty on the coefficient vector. A CP decomposition with the initial rank described above is then applied to the resulting $\hat{\mathcal{B}}^{(\text{Ridge})}$ to obtain the initial model parameters, where $\left(\hat{\mathcal{B}}^{(\text{Ridge})}\right)_{l,c,t}=\left(\hat{}\mbox{$\beta$}_{c,t}^{(\text{Ridge})}\right)_{l}$, for $l=1,\dots,p$.

Algorithm 1 Comprehensive Implementation Algorithm for the Proposed Method

Step 1: Select a sequence of $\lambda$ and $L$. For each choice of $(\lambda,L)$ run the following:

1:$\{y({\mathbf{U}}_{i,j})\}_{j=1,\ldots,N_{i};i=1,\ldots,n}$,$\{{\bm{x}}({\mathbf{V}}_{i,k})\}_{k=1,\ldots,M_{i};i=1,\ldots,n}$, $R_{\max}$, $\lambda$, $h_{i}(L)$

2:(Initialization) Run a weighted Ridge regression applying glmnet and get $\hat{\mathcal{B}}^{(\text{Ridge})}$ as described in Sec 3.2, where $\lambda_{\text{Ridge}}$ is selected via cross-validation (Friedman et al., 2010).

3:while $R_{\max}\geq 1$ do

4:  $R\leftarrow R_{\max}$

5:  Initialize $w_{r}$, ${\bm{q}}^{(1)}_{r}$, ${\bm{q}}^{(2)}_{r}$, and ${\bm{q}}^{(3)}_{r}$ by applying CP on $\hat{\mathcal{B}}^{(\text{Ridge})}$ with initial rank $R$.

6:  repeat

7:   Update parameters as equations from Eq.(3) to Eq.(6).

8:   $w_{r}\leftarrow w_{r}\|{\bm{q}}^{(1)}_{r}\|\|{\bm{q}}^{(2)}_{r}\|\|{\bm{q}}^{(3)}_{r}\|$; ${\bm{q}}^{(m)}_{r}\leftarrow{\bm{q}}^{(m)}_{r}/\|{\bm{q}}^{(m)}_{r}\|$ for $m=1,2,3$.

9:   Update CP-rank $R$ automatically by Eq.(9).

10:  until Stopping rule Eq.(7) satisfied, or $R=0$.

11:  if $R=0$ then

12:   $R_{\max}\leftarrow R_{\max}-1$; Repeat Steps 2-9

13:  else

14:   $R_{\mathrm{final}}\leftarrow R$; break $\triangleright$ converged at current rank $R$   

15:Obtain final estimate as $\mathcal{B}_{\lambda,L}^{(\mathrm{final})}=\sum_{r=1}^{R_{\mathrm{final}}}w_{r}{\bm{q}}^{(1)}_{r}\circ{\bm{q}}^{(2)}_{r}\circ{\bm{q}}^{(3)}_{r}$ for a given $\lambda$ and $L$.

Step 2: Select the optimal $\hat{\lambda}$ from the pre-specified sequence for each choice of $L$ using the criteria in Eq.(10).

Step 3: Select the optimal $\hat{L}$ applying the elbow rule, and get the final estimate $\mathcal{B}^{(\mathrm{final})}_{\hat{\lambda},\hat{L}}$.

4 Simulations

To evaluate our method’s ability to recover the important genes and quantify gene-plaque relationships, we design a simulation setting relying on simulated spatial transcriptomic data using the STRsim package (Zhu et al., 2023). Following the real data, the plaque sites are set well-separated in the simulation too.

We compare our estimates with those from a LASSO-regularized linear regression model. Since fitting LASSO requires the data in a paired form as (outcome, predictors), we first organize the data as $(y_{i},{\bm{x}}_{i})$, where ${\bm{x}}_{i}$ is the gene expression of the cell nearest to $i$-$th$ outcome $y_{i}$. Then, we fit the LASSO regression model.

4.1 Data Generation

While generating the simulated data, we mimic our real data and consider two time points. Thus, for each replicated dataset, we generate two square-shaped SRT samples using SRTsim package, one for each time point. For each sample $i\in\{1,2\}$, let $\mathcal{S}_{i}$ denote the set of spatial spots. For each sample $i$, we simulate a square section with an expected $|\mathcal{S}_{i}|\approx 5000$ spots and $C=3$ cell-type groups (SRTsim “number of groups”=3). SRTsim returns a gene-by-spot expression matrix ${\bm{x}}_{i}\in\mathbb{R}^{p\times|\mathcal{S}_{i}|}$, 2D coordinates $\{{\mathbf{s}}_{i,j}\}_{j\in\mathcal{S}_{i}}$, and group labels $c_{ij}\in\{1,2,3\}$. Each spot has a unique location within the section.

For each sample $i$, we designate a plaque set $\mathcal{P}_{i}\subset\mathcal{S}_{i}$ of size $M_{i}\in\{50,100,200\}$, chosen to be well separated in space and balanced across cell types (i.e., $|\mathcal{P}_{i,c}|\;\;\text{is the same for }c=1,\dots,C$, where $\mathcal{P}_{i,c}=\{j\in\mathcal{P}_{i}:c_{ij}=c\}$). The remaining spots $\mathcal{C}_{i}=\mathcal{S}_{i}\setminus\mathcal{P}_{i}$ serve as the pool of neighboring predictor cells, reflecting the spatial misalignment in real data where outcomes are observed at plaques while predictors are measured at nearby SRT locations.

For each plaque $j\in\mathcal{P}_{i}$, we generate the outcom using the linear model:

$$y_{i,j} ={\bm{x}}(s_{i,j})^{\top}\mbox{$\beta$}_{c,t}+e_{i,j},$$

where $\mbox{$\beta$}_{c,t}\in\mathbb{R}^{p}$ is the regression coefficient for $c$-$th$ cell-type at time $t$. Let $\mathbf{e}_{i}=(e_{i,j})_{j\in\mathcal{P}_{i}}$. Rather than independent noise, we add spatially correlated errors,

$$\mathbf{e}_{i}\sim\mathcal{N}(0,\sigma^{2}_{e}\mathbf{K}_{i}),\qquad\sigma^{2}_{e}\in\{1,5,10,100,200\},$$

with $K_{i}(j,j^{\prime})=\exp\{-d_{i}(j,j^{\prime})/\hat{\phi}\}$. We estimate $\hat{\phi}$ once from the real plaque by fitting an exponential variogram using geoR (variog/variofit; Ribeiro Jr and Diggle (2025)) and then hold $\hat{\phi}$ fixed across samples and replicates. Although $y_{i,j}$ is generated from the expression at the same plaque location ${\bm{x}}(s_{i,j})$, in the downstream analysis we treat plaque-site expression as unobserved and remove $\{{\bm{x}}(s_{i,j}):j\in\mathcal{P}_{i}\}$ from the predictor set, so $\{y_{i,j}:j\in\mathcal{P}_{i}\}$ are modeled using only neighboring non-plaque spots $\{{\bm{x}}(s_{i,j}):j\in\mathcal{C}_{i}\}$.

We set the total number of predictors at 50, with 5 out of 50 having non-zero effects, and assume a CP-based low-rank form with rank 4. In order to achieve this, we set a randomly selected 5 rows in ${\bm{q}}_{1}$ to be non-zero and the rest all zero, with non-zero entries generated from $\mathcal{N}(0,2^{2})$ The entries in ${\bm{q}}_{2}$ and ${\bm{q}}_{3}$ are generated from $\mathcal{N}(5,2^{2})$ and $\mathcal{N}(0,0.5^{2})$ respectively.

4.2 Estimation

We set the $R_{\text{max}}$ at $6$ following maximum rank characterization in Eq.(8), and vary bandwidth over choosing $L\in\{5,10,12,15,20,25,30,35,40,45,50,70\}$. We specify $\lambda_{\text{max}}=10^{10}\times 0.9^{-30}$. We then apply the procedure described in Algorithm 1, which performs data-driven selection of these parameters and yields the final estimate.

To assess variable selection performance, we compute the true positive rate (TPR) and false positive rate (FPR) at threshold $\tau\geq 0$ in Table 1 as

$$\text{TPR}(\tau) =\frac{\sum\mathbb{I}\left\{|\hat{\beta}_{lgt}|\geq\tau,\ \beta_{lgt}\neq 0\right\}}{\sum\mathbb{I}\left\{\beta_{lgt}\neq 0\right\}};\quad\text{FPR}(\tau)=\frac{\sum\mathbb{I}\left\{|\hat{\beta}_{lgt}|\geq\tau,\ \beta_{lgt}=0\right\}}{\sum\mathbb{I}\left\{\beta_{lgt}=0\right\}}.$$

By thresholding the estimated coefficients $\hat{\beta}(i,c,t)$ at a sequence of cutoff values, we trace out an ROC curve in the TPR–FPR plane. We augment it with the extreme points (0,0) and (1,1) to ensure a full range of rate values, and compute the area under the curve (AUC) via the trapezoidal rule.

To quantify overall estimation accuracy, we calculate the mean squared error (MSE) of the coefficient estimates:

$$MSE=\underset{l,c,t}{\mathrm{mean}}\left[\widehat{\beta}(l,c,t)-\beta(l,c,t)\right]^{2}.$$

Across 30 replicates over plaque numbers $M_{i}\in\{50,100,200\}$ and error variances $\sigma_{e}^{2}\in\{1,5,10,100,200\}$, the proposed method consistently delivers better discrimination than the paired LASSO: AUC = 0.54–0.68 (mean 0.606) versus 0.50–0.54 (mean 0.521), and typically lower MSE. The gains persist under high noise and grow with larger $M_{i}$, indicating accurate recovery of the underlying cell-type-specific gene effects. As shown in Table 1, the estimates from our proposed method are better both in estimation accuracy and in selecting the genes with non-zero effects, whereas the paired LASSO often collapses to all-zero solutions.

$M_{i}$	$\sigma_{e}^{2}$	Proposed MSE	Paired MSE	Proposed AUC	Paired AUC
50	1	6.333	7.086	0.574	0.485
50	5	6.367	7.014	0.563	0.485
50	10	6.392	6.985	0.544	0.484
50	100	6.891	8.029	0.495	0.501
50	200	7.359	9.169	0.510	0.507
100	1	6.256	7.455	0.647	0.514
100	5	6.253	7.422	0.629	0.524
100	10	6.275	7.380	0.650	0.525
100	100	6.403	8.330	0.565	0.539
100	200	6.586	9.256	0.568	0.546
200	1	6.278	6.706	0.626	0.484
200	5	6.235	6.717	0.607	0.489
200	10	6.273	6.705	0.606	0.494
200	100	6.319	6.815	0.565	0.501
200	200	6.486	8.356	0.568	0.498

Table 1: Performance over $30$ replicates across choices of plaque number $M_{i}$ and error variance $\sigma_{e}^{2}$, when the data is generated following Section 4.1. For each condition, penalty parameter $\lambda$ and bandwidth $h$ were selected as described in Section 3. Note that paired LASSO often yields all-zero estimates.

5 A$\beta$ Plaque Size Analysis

In this section, we perform the integrative spatial transcriptomic analysis to examine the relationship between gene expression profiles and Amyloid-$\beta$ (A$\beta$) plaque pathology in Alzheimer’s disease (AD), focusing on two different stages of disease progression and three specific cell types. We analyze four SRT tissue sections ($i=1,2,3,4$) from STARmap PLUS, with two sections at each time point (8 and 13 months). Specifically, sections $i=1,2$ correspond to the two 8-month replicates, and sections $i=3,4$ correspond to the two 13-month replicates. The cell counts are $(N_{1},N_{2},N_{3},N_{4})=(8{,}186,\ 8{,}202,\ 10{,}372,\ 9{,}634)$, and plaque counts are $(M_{1},M_{2},M_{3},M_{4})=(88,\ 99,\ 192,\ 136)$. We use the normalized expression provided by Zeng et al. (2023) to place all $p=2{,}766$ genes on a common scale across sections and time points so that estimated coefficients are comparable across genes and over time.

We set the $R_{\text{max}}$ at $6$, and $\lambda_{\text{max}}=10^{5}$. The final selected values are always found to be smaller than these pre-specified upper bounds. We consider the discrete grid ${1,\dots,20,25,30}$ for the number of neighbors $L$, while setting the bandwidth $H_{i}(L)$ for the $i$-$th$ sample; the bandwidth selection procedure is described in Algorithm 1.

Prior to fitting the plaque-size regressions, we apply an expression filter independent of plaque morphology to exclude transcripts with insufficient detection for stable estimation. Specifically, within each of the three target cell types and for each tissue sample collected at 8 and 13 months, we retain genes detected (nonzero) in $>20\%$ of near-plaque cells (within $150\mu$m of the plaque center). We then augment the retained set with (i) 64 marker genes from Zeng et al. (2023) and (ii) plaque-induced genes (PIGs) from Chen et al. (2020) (32 of the 57 PIGs are present in our dataset). This procedure yields a final gene panel of 182 transcripts for modeling plaque morphology.

From our proposed model, we estimate a 3-dimensional coefficient tensor $\mathcal{B}\in\mathbb{R}^{p\times C\times T}$ (gene $\times$ cell type $\times$ time point), with a CP structure to relate plaque size to local, cell-type-resolved expression across time. The final rank is $R=4$ with penalty $\lambda=10^{5}(0.9)^{21}\approx 1.094\times 10^{4}$. For the bandwidth selection, elbow rule give us the choice of $L=5$ for all four samples. Accordingly, we set the sample-specific bandwidth for the kernel function as $h_{i}=H_{i}(5),i=1,2,3,4$, where $H_{i}(L)$ denotes the median distance to the $L$-$th$ nearest neighboring cell in sample $i$ (defined in Sec 3.1). This yields $h=(41.77,45.42,36.89,37.76)\mu$m. Here, each coefficient $\beta(l,c,t)$ is associated with the normalized expression of gene $l$, while setting other variables in the model fixed for cell type $c$ and time $t$.

For each cell type $c$ and time $t$, we summarize the overall strength of association between gene expression in cell type $c$ and the plaque size at time $t$ by $A_{c,t}=\|\beta_{\cdot ct}\|_{2},$ the $\ell_{2}$ norm of all $p$ gene-effects. We summarize the averaged direction of association by $\bar{\beta}_{c,t}=p^{-1}\sum_{l=1}^{p}\beta(l,c,t),$ i.e., the mean gene effect across all $p$ genes for cell type $c$ at time $t$. At 8 months, $A_{c,t}$ is modest (Astrocyte: 2.30; Oligodendrocyte: 6.20; Microglia: 3.04), but it increases sharply by 13 months (Astrocyte: 19.40; Oligodendrocyte: 39.74; Microglia: 21.81), indicating substantially stronger plaque–expression associations over time. The signed means $\bar{\beta}_{c,t}$ are near zero at 8 months (Astrocyte: $-0.0049$, Oligodendrocyte $-0.0018$), while Microglia shows a small positive mean ($0.0285$), meaning that higher microglial gene expression tends to be associated with larger plaques on average at 8 months. By 13 months, the average effects are predominantly positive for Oligodendrocyte (0.0379) and especially Microglia (0.1536), while Astrocyte remains near zero (0.0017). Overall, these results point to strong, increasing plaque-proximal gene-expression effects in oligodendrocyte-lineage cells and microglia over time, with microglia showing the clearest positive directional signal at 13 months.

Following the CP structure, we pursue a component-level interpretation of the results here. CP Decomposes the coefficient tensor into four rank-1 latent factors that summarize the gene–cell-type–time patterns in the association with plaque size (Table 2). The first three components have substantially larger weights than the fourth. For component $r$ with factors ${\bm{q}}_{r}^{(1)}$ (genes), ${\bm{q}}_{r}^{(2)}$ (cell types), ${\bm{q}}_{r}^{(3)}$ (times) and weight $w_{r}$, we orient it so that its top-loading genes have positive entries in the gene factor, and then summarize its overall direction using the signed-average product

$$\text{NetDir}_{r}=\left(\frac{\sum_{l}q_{l,r}^{(1)}}{\sum_{l}|q_{l,r}^{(1)}|}\right)\left(\frac{\sum_{c}q_{c,r}^{(2)}}{\sum_{c}|q_{c,r}^{(2)}|}\right)\left(\frac{\sum_{t}q_{t,r}^{(3)}}{\sum_{t}|q_{t,r}^{(3)}|}\right)\in[-1,1].$$

We refer to 13 months as the late disease stage. As shown in Table 2, each component is dominated by the late time point, with the 13-month time loading having magnitude close to 1, indicating that the estimated associations concentrate their mass at 13 months (with direction determined by the joint sign across modes). Among them, Component 3 shows a clearly positive direction (Net direction $\approx+0.47$) with largest loadings on Oligodendrocytes and Microglia at 13 months; top genes include Trem2, Tmsb4x, Tyrobp, and Plp1, indicating a late microglia–oligodendrocyte component-derived signature in which higher expression is associated with larger plaques. In contrast, Components 1 and 2 carry negative temporal loadings at 13 months and split cell-type emphasis–one tilting toward Oligodendrocytes (Oligodendrocyte ↑, Microglia ↓) and the other toward Microglia (Microglia ↑, Astrocyte ↓). Despite substantial overlap in top genes (e.g., Trem2, C1qa, Aplp1), their opposite cell-type signs imply opposing directions of association across components at the same time point, consistent with heterogeneity between microglia-enriched expression patterns and oligodendrocyte/myelin related pathways. Finally, Component 4 is uniformly negative (Net direction $\approx-1$), dominated by Oligodendrocyte and Microglia loadings at 13 months with genes such as Tmsb4x, Plp1, C1qa, and Cst3, suggesting a late-stage component pattern in which higher expression accompanies smaller plaques. Thus, in summary, the estimated associations concentrate at 13 months and along a Microglia–Oligodendrocyte axis, but with mixed directions across components: Component 3 aligns with larger plaques, whereas Components 1, 2, and 4 align with smaller plaques, suggesting distinct plaque-proximal expression patterns.

These multigene expression patterns arising from the CP-component-based results align well with the known biological knowledge: microglial activation modules (e.g., Trem2, Tyrobp, C1q complex) rise near plaques (Hong et al., 2016; Keren-Shaul et al., 2017; Krasemann et al., 2017), while oligodendrocyte/myelin genes (e.g., Plp1, Mbp) mreflect white-matter processes implicated in AD progression (McKenzie et al., 2017; Nasrabady et al., 2018). CP combines these signals into a small number of interpretable, signed components that are localized by cell type and disease stage, enabling component-level summaries rather than gene-by-gene effects. For instance, a component dominated by microglia at 13 months shows a positive association with plaque size, whereas a component dominated by oligodendrocytes at 13 months shows a negative association.

Component,$r$	Weight, $w_{r}$	NetDir	Top cells	Top times	Top genes
$1$	$98.27$	$-0.1319$	Oligodendrocyte ($+0.75$), Microglia ($-0.51$)	13 mo ($-0.99$), 8 mo ($-0.15$)	Trem2($+0.91$), Tyrobp($-0.24$), C1qa($-0.19$), Aplp1($+0.13$), Mbp($+0.08$)
$2$	$86.34$	$-0.0287$	Microglia ($+0.83$), Astrocyte ($-0.55$)	13 mo ($-0.99$), 8 mo ($-0.13$)	Trem2($+0.92$), C1qa($-0.20$), Aspa($-0.15$), Aplp1($+0.13$), Tyrobp($-0.11$)
$3$	$66.13$	$+0.4662$	Oligodendrocyte ($+0.93$), Microglia ($+0.35$)	13 mo ($+0.99$), 8 mo ($+0.12$)	Trem2($+0.85$), Tmsb4x($+0.31$), Tyrobp($+0.22$), Plp1($+0.17$), Aplp1($+0.15$)
$4$	$19.49$	$-1.0000$	Oligodendrocyte ($-0.85$), Microglia ($-0.52$)	13 mo ($+1.00$), 8 mo ($+0.09$)	Tmsb4x($+0.85$), Plp1($+0.40$), C1qa($+0.29$), Cst3($+0.19$)

Table 2: CP component summary for the estimated coefficient tensor. For each component $r$, “Top cells/times/genes” are those with the largest-magnitude loadings in the corresponding mode. Component signs are oriented per component by top genes. The direction of association is determined by the reconstructed coefficient contribution (positive = larger plaques with higher expression; negative = smaller plaques with higher expression).

Gene-level examples consistent with prior literature (associational). To summarize temporal changes in estimated regression coefficients (associational gene effects), we write $a\!\to\!b$ for gene $l$ in cell type $c$, where $\hat{\beta}_{l,c,8}=a$ (8 months) and $\hat{\beta}_{l,c,13}=b$ (13 months). Astrocytes: transport-related markers show increasingly negative associations with plaque size (Aqp4 $-0.27\!\to\!-1.85$, Slc13a3 $-0.26\!\to\!-1.67$), while alarmin/stress signals become more positively associated (Il33 $+0.27\!\to\!+1.99$, S100a6 $+0.28\!\to\!+1.82$), consistent with reactive astrocyte frameworks and IL-33–microglia crosstalk (Escartin et al., 2021; Vainchtein et al., 2018). Oligodendrocyte lineage: larger plaques are associated with stronger precursor-like signal (Gpr17 $+0.29\!\to\!+2.39$) and reduced mature oligodendrocyte identity (Sox10 $-0.16\!\to\!-0.50$, Aspa $\approx\!-0.20$). These patterns are consistent with reported links between oligodendrocyte function and iron/transferrin biology (Todorich et al., 2009), and OPC–vasculature coupling provides spatial context (Tsai et al., 2016). Microglia: an early Trem2-associated effect attenuates and reverses (Trem2 $+0.64\!\to\!-0.37$), while a later Apoe/Tyrobp-centered signature becomes strongly positively associated with plaque size (Apoe $+0.21\!\to\!+1.78$; Tyrobp $-0.11\!\to\!+1.20$), mirroring reported DAM-like shifts and the TREM2–APOE pathway (Keren-Shaul et al., 2017; Krasemann et al., 2017; Deczkowska et al., 2018).

Novel or less-established associations requiring validation (hypothesis-generating). Oligodendrocyte “logistics”: strong, late positive coefficients for resource intake and protein production (Trf $+0.53\!\to\!+3.53$, Pabpc1 $+0.50\!\to\!+3.55$, Caskin1 $+0.52\!\to\!+3.95$) suggest a resource-mobilizing oligodendrocyte niche as plaques enlarge; in contrast, Flt1 is strongly negative ($-0.57\!\to\!-4.02$), pointing to a vascular/adhesion signature aligned with smaller plaques. Microglial cytoskeleton and lysosome: Rhoc shows a robust negative association ($-0.52\!\to\!-3.95$), consistent with more contractile morphology near smaller plaques, while cathepsins diverge (Ctsl $+0.20\!\to\!+1.04$ vs Ctss $-0.11\!\to\!-0.93$), indicating targeted lysosomal remodeling rather than uniform activation. Astrocytes: persistence and amplification of transport↓/alarmin↑ patterns across time underscore a simple, testable hypothesis: supporting basic astrocytic transport may align with smaller plaques.

6 Conclusion

In this paper, we introduce a novel kernel-weighted method to model how plaque-proximal single-cell spatial transcriptomic expression patterns relate to plaque size. The proposed framework is designed to capture both cell-type-specific and time-varying effects, thereby accommodating heterogeneity across genes, cell populations, and sample collection times. In addition, the method is fully automated through data-driven specification of the tuning parameters, such as rank, penalty, and bandwidth. We further adopt a low-rank model for the regression coefficients, which naturally captures the shared structure and inherent dependencies among genes, cell types, and time points.

Our analysis reveals several plaque-size-associated expression patterns that are strongest at the 13-month data, mainly in microglia (brain immune cells) and oligodendrocytes (myelin-related cells), indicating that plaque-size associations reflect multiple transcriptional signatures acting in parallel rather than a single dominant signature. At a late disease stage, the microglial pattern is characterized by lipid-handling and complement-related genes, consistent with a plaque-proximal activation signature (Hong et al., 2016; Keren-Shaul et al., 2017; Krasemann et al., 2017); oligodendrocyte-lineage cells show a shift toward less mature (precursor-like) signatures, together with increased iron-homeostasis and protein-synthesis transcripts; and astrocyte signals are comparatively weak in net plaque-size tracking, with transport-related genes showing attenuated associations. These findings suggest testable hypotheses for future validation: reduce the late-stage microglial-dominant lipid/complement activity, support astrocyte water/solute transport, and steer oligodendrocytes toward a mature maintenance state to evaluate whether plaque growth can be limited.

For future work, we will analyze how gene expression changes in comparison to the plaque size. This will let us map which genes and multigene expression patterns rise or fall as plaques enlarge, compare patterns across different disease stages and brain regions. These insights will have important translational implications for understanding disease mechanisms and guiding hypothesis-driven experimental studies. Methodologically, there are opportunities to consider more flexible Tucker factorization for $\beta$, and also other types of effect characterizations beyond linearity.

Data availability

We use publicly available data from Zeng et al. (2023) for our analysis.