A framework for simulating interactive particle systems with no underlying geometrical structure, using gillespie SSA with global-bucket optimization, Fisher Information Geometry, and explicit graph topology tracking to derive a statistical distance metric by which to calculate interaction potentials.
Useful for modeling biomolecular dynamics in cellular simulations.
- User Guide (recommended): USER_GUIDE.md
- API Reference: See docstrings in each module
- Mathematical Details: This README (sections below)
Information-Geometric Stochastic Particle Assembly (IGSPA) simulates the self-assembly of particles with colored binding sites without relying on spatial coordinates. The system maintains a dual geometric representation:
| Representation | Type | Purpose |
|---|---|---|
| Implicit Information Manifold | Continuous | Tracks availability fractions θᵢᶜ with Fisher Information Metric |
| Explicit Graph Topology | Discrete | Tracks exact bond structure as a multi-graph |
Standard Gillespie SSA requires O(N²|C|) pairwise evaluation. NS-IGAA achieves O(|C| + |E|) per step by:
- Global Buckets: B_c = {i | particle i has available sites of color c}
- Aggregate Propensities: α_bind^{c,-c} = k_on × (Σᵢ∈B_c availableᵢᶜ) × (Σⱼ∈B₋c availableⱼ₋c)
- Two-Stage Sampling: Channel → Weighted particle selection
Each particle i has availability fractions θᵢᶜ = available/capacity. The manifold uses the Shahshahani form of the Fisher Information Metric:
ds² = Σ_c (dθᶜ)² / θᶜ
As θ → 0, distances diverge → ∞, geometrically encoding scarcity.
Bonds form a multi-graph G(t) = (V, E(t)) with colored edges. Conservation law:
θᵢᶜ(t) = 1 - (Σⱼ Aᵢⱼᶜ(t)) / Kᵢᶜ
| Metric | Standard SSA | IGSPA (Global Bucket) |
|---|---|---|
| Time/step | O(N² | C |
| Space | O(N²) | **O(N |
| Scaling | Quadratic | Linear/Constant |
Dependencies: numpy>=1.24, networkx>=3.0, scipy>=1.10, lxml>=6.0 (for GraphML export)
from igspa import create_system_from_strings
# Create system from string color lists
system = create_system_from_strings([
['white', 'red'], # Particle 0
['green', 'blue'], # Particle 1
['black', 'green', 'yellow'], # Particle 2
['red', 'blue'], # Particle 3
], k_on=1.5, k_off=0.3)
# Run simulation
stats = system.run(max_steps=100, verbose=True)
# Inspect complexes
system.print_complexes()Output:
=== COMPLEXES (Time: 7.878s, Step: 12) ===
Complex 1 (Size 1): P0
Complex 2 (Size 3): P1 -- P2 -- P3
Each particle has binding sites of specific colors. Colors come in complementary pairs:
white↔blackred↔greenblue↔yellow
A particle can bind to another if it has an available site of a color complementary to the other particle's available site.
NS-IGAA maintains two synchronized representations:
-
Implicit Information Manifold (Continuous)
- Tracks availability fractions θᵢᶜ ∈ [0, 1] for each particle and color
- Equipped with Fisher Information Metric (Shahshahani geometry)
- As θ → 0, informational distance → ∞ (scarcity geometry)
-
Explicit Graph Topology (Discrete)
- Multi-graph G(t) = (V, E(t)) with colored edges
- Tracks exact bond structure
- Enables complex detection via connected components
- TwoStageSampler (default): O(|C| + |E|) two-stage Gillespie
- CompositionRejectionSampler: O(1) dyadic binning for large systems (N > 50)
- SteadyStateAnalyzer: Equilibrium constants, gel fraction, cluster sizes
- Manifold Geometry: Geodesic distances, saturation distances, metric condition numbers
- Thermodynamics: Entropy, free energy, chemical potentials
- Ensemble Statistics: Multi-trajectory statistical analysis
- Logging: JSONL streaming, NumPy .npz, Pandas DataFrames
- Export: JSON, Pickle, NPZ, GraphML (Gephi/Cytoscape), CSV
- Checkpointing: Full state snapshots for resumable simulations
- Standard 6-color palette (3 complementary pairs)
- Custom palettes (e.g., DNA: A↔T, C↔G)
- String-based color specification
src/igspa/
├── __init__.py # Public API exports
├── core/
│ ├── colors.py # Color enum, palette, complementarity
│ ├── particle.py # ParticleBlueprint, Particle, factories
│ └── system.py # ParticleSystem, SimulationConfig, factory
├── geometry/
│ └── manifold.py # InformationManifold, FisherMetric
├── topology/
│ └── graph.py # BondGraph, Bond, ComplexAnalyzer
├── algorithms/
│ └── gillespie.py # GlobalBucketSampler, TwoStageSampler, CR Sampler
├── io/
│ └── serialization.py # SimulationLogger, StateExporter
└── utils/
└── analysis.py # Validation, SteadyStateAnalyzer, ensemble tools
Let a particle system be defined by the tuple
-
$\mathcal{V} = {1, \dots, N}$ is the bounded set of particle vertices. -
$\mathcal{C} = {c_1, c_2, \dots, c_m}$ is the alphabet of binding site colors. A bijective mapping function$\text{comp}: \mathcal{C} \to \mathcal{C}$ satisfies comp(c) = -c and comp(-c) = c. -
$\mathbf{K}_i = [K_i^{c_1}, \dots, K_i^{c_m}]^T \in \mathbb{N}^m$ specifies the static total site capacities. -
$\boldsymbol{\theta}_i(t) = [\theta_i^{c_1}(t), \dots, \theta_i^{c_m}(t)]^T \in [0, 1]^m$ details the continuous tracking coordinate on the information manifold$\mathcal{M}_i$ , representing the fraction of active, unbonded sites:
The distance element on this manifold is governed by the localized Shahshahani form of the Fisher Information Metric Tensor (
The actual physical structure of the system is a combinatorial multi-graph
Standard stochastic simulation loops over all pairs
The global binding propensity (
The global breaking propensity scales linearly with the cardinal volume of the active graph edges:
The sum of all active pathways yields the system parameter
Each simulation iteration proceeds via the following sequence:
-
Time Delta Evacuation: Draw next event occurrence interval
$\tau = \frac{1}{a_0} \ln\left(\frac{1}{r_1}\right)$ , where r₁ ~ Uniform(0,1). -
Channel Selection: Draw a random channel from the combined propensity array
$\mathbf{\alpha}$ proportional to its total weight. -
Two-Stage Node Resolution (If Binding Event Chosen):
- Stage 3a: Extract an index i from bucket
$\mathcal{B}_c$ with probability proportional to its localized residual capacity:$P(i) \propto K_i^c \theta_i^c$ . - Stage 3b: Extract an index j from bucket
$\mathcal{B}_{-c}$ with probability proportional to its localized residual capacity:$P(j) \propto K_j^{-c} \theta_j^{-c}$ .
- Stage 3a: Extract an index i from bucket
-
State Mutation Execution: Append an edge to
$\mathcal{E}$ connecting (i, j). Update internal parameters$\theta_i^c$ and$\theta_j^{-c}$ . If a parameter reaches zero, purge that node index from its corresponding bucket$\mathcal{B}$ in$\mathcal{O}(1)$ time.
To introduce dependencies where a bond formation at one site alters the affinity of another site on the same particle, define an Allosteric Configuration Tensor (
This dynamically updates the global buckets
While the system is non-spatial, checking graph properties can introduce virtual geometry constraints. For instance, to penalize or favor the formation of closed molecular rings (cycles), the algorithm can check the shortest path distance between candidate nodes i and j in the explicit graph before confirming a bond. If a path already exists, the binding rate can be scaled by a cyclization factor γ:
MIT License - See LICENSE file for details.
Information-Geometric Stochastic Particle Assembly (IGSPA)
Carson Scott, 2024