Spatially Varying Graphical Models for Cell-Cell Interaction Networks in Multiplexed Tissue Imaging

The paper introduces GP-GHS, a scalable Bayesian framework that models cell-cell interactions as spatially varying networks using Gaussian processes and group horseshoe priors, enabling the detection of context-specific immunosuppressive mechanisms in colorectal cancer that existing global or marginal methods miss.

The Gist

The Big Picture: Mapping the City of Cells

Imagine a tumor isn't just a lump of bad cells, but a bustling, chaotic city. In this city, there are different neighborhoods: the "Tumor Core" (the downtown business district), the "Invasive Margin" (the suburbs), and the "Stroma" (the parks and infrastructure).

In this city, there are many types of citizens:

• Tumor cells (the troublemakers).
• T-cells (the police trying to stop the trouble).
• Macrophages (the janitors or meddlers who can either help or hinder).
• Tregs (the peacekeepers who sometimes stop the police from doing their job).

For a long time, scientists could take a photo of this city and count how many of each citizen lived there. But they couldn't tell who was talking to whom. They knew the police and the troublemakers were in the same city, but they didn't know if they were actually interacting, or if they were just both hanging out near the same coffee shop (a third party) by coincidence.

The Problem: The "Group Chat" Mistake

Previous methods for studying these cities were like looking at a group chat log and saying, "Oh, Person A and Person B both sent messages at 2:00 PM, so they must be friends."

This is flawed because:

1. It ignores the third person: Maybe Person A and Person B both sent messages because Person C (the coffee shop owner) told them to. They aren't friends; they just both like the coffee shop.
2. It assumes everyone is the same everywhere: It assumes the rules of friendship in the "Tumor Core" are exactly the same as in the "Suburbs." But in reality, interactions change depending on where you are in the tissue.

The Solution: GP-GHS (The Smart City Planner)

The authors created a new tool called GP-GHS. Think of it as a super-smart city planner who can look at a map of the city and draw a dynamic map of who is actually talking to whom, and how that conversation changes as you walk from the city center to the edge.

Here is how it works, broken down into three simple parts:

1. The "Smooth Map" (Hilbert Space Gaussian Process)

Imagine trying to draw a map of temperature across a city. If you just measured the temperature at 10 specific spots, you'd have a jagged, messy map.

• Old way: Connect the dots with straight lines.
• GP-GHS way: It uses a "smooth brush" (a mathematical tool called a Gaussian Process) to paint a continuous, flowing map. It knows that temperature doesn't jump instantly from hot to cold; it flows. This allows the model to see that a cell interaction might be strong in one corner of the tissue and weak in another, rather than forcing a single "average" answer for the whole tissue.

2. The "All-or-Nothing" Rule (Group Horseshoe Prior)

This is the most important trick.
Imagine you are trying to figure out if two people, Alice and Bob, are friends. You have 20 different clues (like: did they walk together? did they eat lunch together? did they text?).

• Old way: If any single clue looks slightly suspicious, you say "Maybe they are friends." This leads to lots of false alarms.
• GP-GHS way: It treats all 20 clues as a single package. It asks: "Is the entire package of evidence strong enough to say they are friends?"
  • If the evidence is weak, it shrinks the whole package to zero (No, they aren't friends).
  • If the evidence is strong, it keeps the whole package (Yes, they are friends).
  • Why this matters: It prevents the model from getting confused by random noise. It forces a clear "Yes" or "No" decision on the relationship, rather than a fuzzy "maybe" on every little detail.

3. The "Parallel Workers" (Nodewise Regression)

Instead of trying to solve the puzzle of the whole city at once (which is too hard), the model breaks the city into 15 smaller puzzles (one for each cell type). It hires 15 different workers to solve their own small puzzle simultaneously. Then, it combines their answers to build the final map. This makes the math fast enough to run on a normal computer.

What Did They Find? (The Colorectal Cancer Discovery)

The team tested this tool on real data from Colorectal Cancer patients. They looked at two types of "cities" (tumor microenvironments):

1. CLR: A city with organized, distinct neighborhoods (like a planned suburb).
2. DII: A city with a chaotic, scattered mix of people (like a crowded festival).

The Discovery:
Using GP-GHS, they found a massive difference in the "DII" cities.

• In the DII cities, the Tregs (the peacekeepers) were tightly connected to everyone. They were talking to the police (T-cells), the janitors (Macrophages), the troublemakers (Tumor cells), and even the infrastructure (Blood vessels).
• This created a suppression network. The Tregs were essentially putting a "Do Not Disturb" sign on the whole immune system, stopping the police from fighting the cancer.
• In the CLR cities, these connections were weak or non-existent. The immune system was more compartmentalized and less suppressed.

Why This Matters

Before this tool, scientists might have just said, "There are a lot of Tregs in the DII tumors."
With GP-GHS, they can say: "The Tregs are actively forming a secret alliance with the Macrophages and the Tumor cells specifically in the DII tumors, and this alliance is the reason the immune system is failing."

This changes the game. Instead of just counting cells, we can now map the social network of the tumor. This helps doctors understand why some tumors are resistant to treatment and might point to new ways to break up these "secret alliances" to help the body fight the cancer.

Summary Analogy

• Old Methods: Counting how many people are in a room and guessing who is friends based on who is standing near the same wall.
• GP-GHS: Putting a camera on every person, recording their conversations, and drawing a live, moving map of who is actually talking to whom, while ignoring the background noise.

The paper proves that this new method is much better at finding the real connections and ignoring the fake ones, especially when the "city" is messy and complex.

Technical Summary

1. Problem Statement

The tumor microenvironment (TME) is a spatially heterogeneous ecosystem where cell-cell interactions vary across different tissue compartments (e.g., tumor core vs. invasive margin). While multiplexed tissue imaging (e.g., CODEX, Imaging Mass Cytometry) provides single-cell resolution data for dozens of cell types, existing computational methods fail to adequately model the conditional interaction networks governing these systems.

Current approaches suffer from three main limitations:

1. Marginal Statistics: Methods like neighborhood enrichment (e.g., Squidpy, Giotto) rely on pairwise co-occurrence without conditioning on third cell types, leading to spurious correlations.
2. Global Homogeneity: Existing graphical model methods (e.g., SVCA, SpaCeNet) estimate a single global interaction coefficient per cell type pair, ignoring the fact that interactions vary spatially across the tissue.
3. Incoherent Edge Selection: Standard Bayesian shrinkage priors applied to spatial basis functions treat coefficients independently. This fails to enforce a coherent "all-or-nothing" decision for edge existence, resulting in high false discovery rates (FDR) when aggregating multiple basis coefficients to represent a single spatial field.

Goal: Develop a framework to infer sparse, spatially varying conditional independence graphs where edge strengths are functions of location, and edge inclusion is a group-level decision.

2. Methodology: GP-GHS

The authors propose GP-GHS (Gaussian Process - Group Horseshoe), a Bayesian nodewise regression framework. The method consists of three core components:

A. Nodewise Spatial Regression

Instead of estimating a global precision matrix, the problem is decomposed into $p$ parallel nodewise regressions (one for each cell type $s$):
$$Y_s(l) = \sum_{s' \neq s} \beta_{s'}(l) Y_{s'}(l) + \epsilon(l)$$
Here, $\beta_{s'}(l)$ is a spatially varying coefficient function representing the conditional dependency between cell type $s$ and $s'$ at location $l$. An edge exists if $\beta_{s'}(\cdot)$ is not identically zero.

B. Hilbert Space Gaussian Process (HSGP) Approximation

To model $\beta_{s'}(l)$ as a Gaussian Process (GP) while maintaining scalability:

• The authors use the HSGP approximation (Riutort-Mayol et al., 2023).
• The spatial field is expanded using eigenfunctions of the Laplacian operator (spectral basis functions $\phi_{jk}$) on a bounded domain.
• $\beta_{s'}(l) \approx \sum_{j,k} w_{jk}^{(s')} \phi_{jk}(l)$.
• This reduces the computational complexity of exact GP inference from $O(n^3)$ to $O(nm^2)$, where $n$ is the number of spatial locations and $m$ is the number of basis functions per dimension.
• A Matérn 3/2 kernel is used to enforce spatial smoothness, allowing for sharp transitions at tissue boundaries.

C. Group Horseshoe Prior for Edge Selection

This is the critical innovation for handling edge selection in a spatial context.

• The Problem: If independent shrinkage parameters are applied to each spectral basis coefficient ($w_{jk}$), the model cannot make a unified decision about whether the entire interaction field (edge) exists.
• The Solution: A Group Horseshoe Prior is placed on the block of coefficients corresponding to a specific neighbor $s'$.
  • A single local shrinkage parameter $\lambda_j$ is shared across all $m^2$ basis coefficients for a given edge.
  • Mechanism: If $\lambda_j \to 0$, the entire block of coefficients collapses to zero simultaneously, declaring the edge absent globally. If $\lambda_j$ is large, the coefficients are free to vary, governed by the spectral prior.
  • This separates the roles: $\lambda_j$ controls edge existence (group level), while the spectral prior controls spatial smoothness (within-edge level).
• Graph Construction: The $p$ nodewise regressions are fit independently. A symmetric undirected edge $(s, s')$ is included only if both directional regressions ($s \to s'$ and $s' \to s$) declare the edge active (AND rule).

D. Posterior Inference

• Inference is performed via a closed-form block Gibbs sampler.
• The sampler cycles through updates for regression coefficients ($\theta$), group shrinkage parameters ($\lambda$), global shrinkage ($\tau$), and noise variance ($\sigma^2$).
• The $p$ nodewise regressions are parallelized across CPU cores.

3. Key Contributions

1. First Spatially Varying Graphical Model: GP-GHS is the first framework to simultaneously estimate sparse conditional independence graphs and allow edge strengths to vary continuously across tissue space.
2. Group Horseshoe Innovation: The paper demonstrates that group-level shrinkage is mathematically necessary for spatial signal aggregation. Ablation studies show that replacing the group prior with a standard scalar horseshoe prior causes performance to collapse, as independent shrinkage cannot distinguish a weak spatial signal from noise.
3. Scalable HSGP Integration: By combining HSGP with group shrinkage, the method achieves computational tractability for tissue sections with hundreds of spots, avoiding the $O(n^3)$ bottleneck of exact GPs.
4. Differential Testing Framework: The authors propose a Linear Mixed Model (LMM) approach to test for differences in interaction networks between pathological groups, using the posterior shrinkage score ($\kappa$) as the response variable to account for within-patient correlation.

4. Results

Simulation Studies

• Setup: Simulated data with $n=600$ locations and $p=15$ or $25$ cell types, using scale-free graphs (Barabási-Albert) to mimic tumor immune networks.
• Performance: GP-GHS significantly outperformed all competitors (Graphical Lasso, Nodewise Lasso, Standard Horseshoe, Exact GP, Correlation Threshold) in F1 score and Matthews Correlation Coefficient (MCC) across all sparsity levels.
• Key Findings:
  • Exact GP (EGP): While EGP had high sensitivity (TPR), it suffered from catastrophic False Discovery Rates (FDR) in sparse graphs because it lacked global regularization.
  • Standard Horseshoe (SHS): Failed completely (F1 $\approx$ 0) because independent shrinkage on basis coefficients could not enforce coherent edge decisions.
  • GP-GHS: Maintained high sensitivity while controlling FDR, proving that the combination of spatial smoothness (GP) and group sparsity (Horseshoe) is essential.
  • Computational Cost: GP-GHS was significantly faster than Exact GP (seconds vs. tens of seconds) due to the HSGP approximation.

Real Data Application: Colorectal Cancer (CRC)

• Dataset: 140 CODEX images from 35 CRC patients, stratified into two subtypes: Crohn's-like reaction (CLR) and Diffuse Inflammatory Infiltrate (DII).
• Preprocessing: Cell counts were converted to continuous spatial abundance surfaces using 2D Kernel Density Estimation (KDE) on a $30 \times 30$ grid.
• Findings:
  • GP-GHS identified 13 differentially active edges (FDR < 0.05) between the two subtypes.
  • Biological Insight: All significant edges showed stronger spatial interaction in the DII subtype.
  • Treg-Centered Network: The DII subtype exhibited a highly amplified immunosuppressive network centered on Regulatory T cells (Tregs). Tregs showed strong spatial co-localization with CD68+ macrophages, CD8+ T cells, B cells, vasculature, and stroma.
  • CLR Subtype: Showed a globally sparse interaction structure, consistent with the known "compartmentalized" lymphoid aggregates of Crohn's-like reactions.
  • Validation: These findings align with known biological mechanisms where DII tumors are characterized by diffuse, immunosuppressive Treg recruitment, whereas CLR tumors have more focal immune infiltration.

5. Significance

• Methodological Advancement: GP-GHS fills a critical gap in spatial omics analysis by providing a rigorous statistical tool to move beyond static, global interaction maps to dynamic, spatially resolved networks.
• Biological Discovery: The application to CRC data reveals that the spatial architecture of immune suppression (specifically Treg-macrophage coupling) is a defining feature of the aggressive DII subtype, offering new mechanistic hypotheses for therapeutic resistance.
• Practical Utility: The method is computationally efficient enough for whole-slide imaging datasets and provides a framework for differential network analysis that accounts for patient-level correlation, a common challenge in clinical imaging studies.

In conclusion, the paper establishes that modeling spatial heterogeneity and enforcing group-level sparsity are not just optional enhancements but fundamental requirements for accurately recovering cell-cell interaction networks from multiplexed tissue imaging data.