---
title: Neural Particle Automata
url: https://www.emergentmind.com/topics/neural-particle-automata-npa
type: topic
---

# Neural Particle Automata

Neural Particle Automata (NPA) are a Lagrangian extension of Neural Cellular Automata (NCA), generalizing the concept of neural automata from static lattice-based cells to dynamic systems of particles with continuous positions and internal states. NPA replaces the Eulerian paradigm—where state is tied to fixed pixels or voxels—with a particle-based approach wherein each particle is individually simulated and interacts locally via neural rules. This framework enables explicit individuation, heterogeneous dynamics across particles, and focused computation in regions of activity. A differentiable Smoothed Particle Hydrodynamics (SPH) formulation, with CUDA-accelerated kernels, allows NPA to scale efficiently to large particle ensembles and supports end-to-end learning of self-organizing behaviors across morphogenesis, texture synthesis, and classification tasks [2601.16096].

## 1. Formal Model Specification

NPA operates on a set of $N$ particles indexed by $i=1,\dots,N$ in $D$-dimensional space. Each particle is defined by a continuous position $x_i \in \mathbb{R}^D$ and an internal state $S_i \in \mathbb{R}^C$ representing hidden channels. The evolution of the system is governed by a shared, learnable local update rule $f_\theta$, typically implemented as a two-layer MLP. At each discrete timestep $t$:

1. For each particle $i$, a local perception vector $Z_i^t$ is constructed from its $\epsilon$-neighborhood $N(i)$:
   - $S_i^t$: current state
   - $\bar{S}_i^t$: density-normalized, locally smoothed state
   - $\nabla S_i^t$: state gradient
   - $\nabla \rho_i^t$: density gradient
   - Expressed as $Z_i^t = [ S_i^t ;  \bar S_i^t ;  \nabla S_i^t ;  \nabla\rho_i^t ]$

2. The neural rule computes $\Delta y_i^t = f_\theta(Z_i^t)$, yielding updates for state (static) or both position and state (dynamic tasks).

3. A stochastic update is applied per particle using $\delta_i^t \sim \mathrm{Bernoulli}(p)$, typically $p = 0.5$ for asynchrony:
   - $S_i^{t+1} = S_i^t + \delta_i^t \cdot \Delta S_i^t$
   - $x_i^{t+1} = x_i^t + \delta_i^t \cdot \Delta x_i^t$

This model permits both heterogeneity and locality in updates, supporting individuation and dynamic adaptation even in unstructured particle ensembles.

## 2. Differentiable SPH Perception Operators

To enable fully differentiable, locality-aware perception and interaction, NPA utilizes SPH for estimating densities, smoothing states, and calculating gradients:

- **Density Estimation:**
  $$
  \rho_i = \sum_{j \in N(i)} m_j W(\|r_{ji}\|; \epsilon)
  $$
  where $r_{ji} = x_j - x_i$, kernel $W$ is Poly6 in 2D, and $m_i$ is usually normalized $(m_i = 1/N)$.

- **State Smoothing:**
  $$
  \bar{S}_i = \sum_{j \in N(i)} \frac{m_j}{\rho_j} S_j W(\|r_{ji}\|; \epsilon)
  $$

- **Gradient Estimation (Difference Form):**
  $$
  \nabla_0 S_i = \sum_{j \in N(i)} \frac{m_j}{\rho_j} (S_j - S_i) \otimes W^\nabla(r_{ji}; \epsilon)
  $$

- **Gradient Correction (Moment Matrix):**
  $$
  M_i = \sum_j \frac{m_j}{\rho_j} r_{ji} W^\nabla(r_{ji}; \epsilon)^\top
  $$
  $$
  \nabla_1 S_i = M_i^{-1} \nabla_0 S_i
  $$
  In practice, this is implemented as: $\nabla S_i = \nabla_0 S_i + \mathrm{detach}(\nabla_1 S_i - \nabla_0 S_i)$ to ensure numerical stability in backpropagation.

- **Density Gradient:**
  $$
  \nabla \rho_i = \sum_{j \in N(i)} m_j W^\nabla(r_{ji}; \epsilon)
  $$

### SPH Kernel Functions

| Kernel          | Formula                                                | Notes                |
|-----------------|--------------------------------------------------------|----------------------|
| Poly6 $W(r;\epsilon)$ | $\frac{4}{\pi \epsilon^8}(\epsilon^2 - r^2)^3$  | Support: $0 \leq \|r\| \leq \epsilon$ |
| Spiky Gradient $W^\nabla(r;\epsilon)$ | $\frac{10}{\pi \epsilon^5}(\epsilon - \|r\|)^2 (r/\|r\|)$ | Support: $0 < \|r\| \leq \epsilon$      |

A uniform hash-grid bins particles by position for efficient neighbor querying; two CUDA kernel strategies are implemented:
- **Particle-centric**: One thread per particle, scanning adjacent grid cells.
- **Grid-centric**: One block per cell, staging per-cell particles in shared memory for coalesced access.

Neighborhoods are discovered on-the-fly; memory grows $\mathcal{O}(N)$ and computational cost is $\mathcal{O}(\text{avg neighbors})$ per particle.

## 3. Training Regime and Loss Construction

End-to-end NPA training employs backpropagation through SPH kernels and adapts several practices from NCA:

- Persistent **state pool** for sampling simulation rollouts.
- **Variable rollout lengths** with $T \sim [T_{\min}, T_{\max}]$.
- **Overflow regularization** with $L_2$ penalties enforcing $S_i \in [-1,1]$.
- **Vector-input scaling** for stability:
  $$
  v \leftarrow \log(1 + |v|) \cdot \frac{v}{|v| + \eta}
  $$
- **Stop gradients** through positions $x$ for perception computations.
- **Displacement regularization**: penalizes $\sum_t \|\Delta x_i^t\|$ for smooth motion.

### Representative Tasks and Losses

| Task                         | Particle Count ($N$) | State Channels ($C$) | Loss Function                                      |
|------------------------------|---------------------|----------------------|---------------------------------------------------|
| 2D Morphogenesis             | 4096                | 16                   | $L_{1+2}(D, D^*) + \lambda \cdot L_{1+2}(C, C^*) \cdot \text{detach}(e^{-|D-D^*|})$ |
| 3D Morphogenesis             | 16384               | 24                   | Multi-scale SSIM + isotropy regularizer ($\|A_i\|_1$) |
| Texture Synthesis (particles)| 4096                | 16                   | VGG OT-loss (on $C$ and $D$ across channels), power law transform $\overline{D} = (D+\epsilon)^\alpha, \alpha = 0.1$ |
| Self-classifying Point Clouds| 512                 | 16                   | $L_2$ on one-hot labels, result: $98.42\%$ test accuracy |

Each task utilizes customized rendering—Gaussian splatting for morphogenesis and textures, GSplat-based multi-view for 3D—plus curriculum strategies (progressive difficulty via blurred/sharpened targets) and per-task regularization.

## 4. Observed Dynamics and Empirical Behavior

NPA demonstrates several salient behaviors across evaluated tasks:

- **Robustness to Discretization**: Dynamics remain stable across varied $\epsilon \in [0.05, 0.2]$ and particle number $N$, degrading only under extreme downsampling.
- **Stochastic Updates**: Varying Bernoulli $p$ during inference (from $0.1$ to $1.0$) preserves morphogenetic convergence.
- **Self-Regeneration**: Local perturbations (state zeroing, slit cuts, clumping) are autonomously repaired by continued application of local update rules.
- **Heterogeneous Multi-Species Interaction**: Independently trained NPA rules interact within the same simulation, exhibiting cooperation, mixing, and disruption, despite no joint training.
- **Hidden State Structure and Persistent Flows**: Visualization of RGB-projected hidden states reveals internally consistent, directional channel patterns ("rainbow" flows) and persistent vortex-like movement post-target attainment, which suggests emergent memory/maintenance mechanisms.
- **PointMNIST Dynamics**: Hidden state clusters, tracked by UMAP, form per digit class over time, indicating emergence of global semantic structure from local exchanges.

## 5. Essential Algorithms and Mathematical Formulations

Key pseudocode for NPA update:

```python
# Given particles {x_i, S_i}, for t = 1...T:
1. Build hash-grid index from {x_i}
2. For each particle i (in parallel):
    # Gather neighborhood
    N(i): all j in 3^D grid cells around x_i, with ||x_j - x_i|| < ε
    # Compute local perception via SPH
    ρ_i, \bar S_i, ∇S_i, ∇ρ_i via SPH sums
    # Scale vectors
    ∇S_i ← scale(∇S_i), ∇ρ_i ← scale(∇ρ_i)
    # Combine perception
    Z_i = [S_i; \bar S_i; ∇S_i; ∇ρ_i]
    # Execute neural rule
    Δy_i = f_θ(Z_i)
3. For each i:
    δ_i ∼ Bernoulli(p)
    S_i ← S_i + δ_i·ΔS_i ; x_i ← x_i + δ_i·Δx_i
```

Summary mathematical constructs central to NPA:

- SPH Poly6 kernel:
  $$
  W(r;\epsilon) = \frac{4}{\pi \epsilon^8} (\epsilon^2 - r^2)^3
  $$
- Spiky gradient kernel:
  $$
  W^\nabla(r;\epsilon) = \frac{10}{\pi \epsilon^5} (\epsilon - r)^2 \frac{r}{\|r\|}
  $$
- Logarithmic scaling:
  $$
  \text{scale}(v) = \log(1+|v|) \cdot \frac{v}{|v|+\eta}
  $$

CUDA constants:
- Each SPH kernel (forward/backward) scans $3^D$ cells per particle.
- Particle-centric (thread/particle), grid-centric (block/cell) strategies—no explicit adjacency lists, neighborhoods discovered per step, memory growth $\mathcal{O}(N)$.

## 6. Theoretical and Practical Significance

NPA synthesizes compactness and robustness characteristic of NCA with the flexibility and resolution of particle-based, Lagrangian models. By leveraging differentiable SPH operators and efficient CUDA implementations, NPA accommodates large-scale, self-organizing particle dynamics with local computation and persistent individuation. Its empirical properties—stability under discretization, resilience to stochastic updates, self-regenerative capacity, and emergent structure—suggest potential utility in a range of simulation, synthesis, and classification domains. The particle framework inherently avoids global synchronization, handles unstructured domains, and empowers the learning of both robust and heterogeneous behaviors from local rules [2601.16096].

Source: https://www.emergentmind.com/topics/neural-particle-automata-npa