Metadata-Version: 2.4
Name: DynaVista
Version: 0.5.0rc1
Summary: Stochastic discrete dynamical systems simulator for Boolean and multi-state biological networks
Author: Daniel Plaugher
License: MIT
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Requires-Dist: sympy>=1.12
Requires-Dist: matplotlib>=3.7
Requires-Dist: pandas>=2.0
Requires-Dist: joblib>=1.3
Requires-Dist: openpyxl>=3.1
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: jupyter>=1.0; extra == "dev"
Requires-Dist: ipykernel; extra == "dev"
Provides-Extra: docs
Requires-Dist: mkdocs>=1.5; extra == "docs"
Requires-Dist: mkdocs-material>=9.5; extra == "docs"
Requires-Dist: mkdocstrings[python]>=0.24; extra == "docs"
Requires-Dist: black>=24.0; extra == "docs"
Dynamic: license-file

# DynaVista

**Stochastic Discrete Dynamical Systems (SDDS) Simulator for Biological Networks**

A Python framework for Boolean and multi-state network simulation, mutation/control
analysis, attractor identification, and Markov chain analysis.

**Authors:** Daniel Plaugher  
**Based on:** SMATA pipeline (Plaugher 2022) and BNPBN toolbox (Shmulevich, Aguilar, Murrugarra)  
**Version:** 0.5.0

[![PyPI](https://img.shields.io/pypi/v/DynaVista?label=PyPI)](https://pypi.org/project/DynaVista/)
[![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](LICENSE)
[![Docs](https://img.shields.io/badge/docs-mkdocs-green)](https://drplaugher.github.io/DynaVista)

---

## Table of Contents

1. [Installation](#1-installation)
2. [Quick Start](#2-quick-start)
3. [Project Structure](#3-project-structure)
4. [Core Concepts](#4-core-concepts)
5. [Module Reference](#5-module-reference)
6. [Config File Guide](#6-config-file-guide)
7. [Network Format Guide](#7-network-format-guide)
8. [Mutation and Control Guide](#8-mutation-and-control-guide)
9. [Mutation Sweep Guide](#9-mutation-sweep-guide)
10. [Markov Chain Guide](#10-markov-chain-guide)
11. [Scalability Guide](#11-scalability-guide)
12. [MATLAB Equivalence Table](#12-matlab-equivalence-table)
13. [Citation](#13-citation)

---

## 1. Installation

### Requirements

- Python 3.10 or higher
- pip (included with Python)

### Step-by-step install

**Step 1.** Download the `dynavista` folder (from GitHub or as a zip).

**Step 2.** Open a terminal and navigate to the folder containing `dynavista/` and `pyproject.toml`.

**Step 3.** Install:

```bash
pip install -e .
```

The `-e` flag installs in *editable* mode — changes to source files take effect immediately without reinstalling.

All required packages are installed automatically:

| Package | Purpose |
|---------|---------|
| `numpy` | Array math, state vectors |
| `scipy` | Sparse matrices (Markov chain) |
| `sympy` | Symbolic polynomial evaluation in sdds_build |
| `matplotlib` | All plotting |
| `pandas` | Sweep results, attractor tables |
| `joblib` | Parallel simulation runs |
| `openpyxl` | Excel workbook export (config pipeline) |

### Verify

```python
import dynavista
print(dynavista.__version__)  # should print 0.5.0
```

### Documentation

Full documentation — including the auto-generated API reference — can be read
two ways:

**Option A — Read it online (no setup).**
The published documentation is hosted at
**[https://drplaugher.github.io/DynaVista](https://drplaugher.github.io/DynaVista)**.
This is the easiest option: just open the link in a browser. (If the link does
not load yet, the online docs have not been published for this release — use
Option B in the meantime.)

**Option B — Build and read it locally.**
The documentation is built with [MkDocs](https://www.mkdocs.org/) from the
`docs/` folder and the source docstrings. To read it on your own machine:

```bash
pip install -e ".[docs]"     # installs mkdocs + the docs toolchain
mkdocs serve                 # serves the docs locally
```

Then open **http://127.0.0.1:8000** in your browser. Press `Ctrl+C` in the
terminal to stop the server. To produce a static copy instead of serving it,
run `mkdocs build` — the finished site is written to a `site/` folder you can
open directly (`site/index.html`).

The API reference is generated automatically from the docstrings in the source
code, so it always matches the installed version.

### IDE setup

DynaVista works in any Python environment with no special configuration:
- **Jupyter Notebook / JupyterLab** — recommended for interactive analysis
- **VS Code** — install the Python extension, select your Python interpreter
- **PyCharm** — open the project folder, mark `dynavista/` as a source root
- **Spyder** — install via Anaconda, then `pip install -e .` in the Spyder terminal
- **Plain Python scripts** — `python my_analysis.py`

---

## 2. Quick Start

### Minimal example

```python
import numpy as np
from dynavista import sdds_build, simulate, plot_simulation

# Define network
node_names = ["A", "B", "C"]
rules = ["C", "A", "A OR B"]         # Boolean or polynomial form

# Build SDDS structures
varF, nv, F = sdds_build(rules, node_names, p=2)

# Propensity matrix (probability of state change per step)
n = len(node_names)
c = 0.9 * np.ones((2, n))

# Simulate — noise=0.0 selects the clean path automatically
Y, My = simulate(F, varF, nv, p=2, c=c, nsteps=20, nins=1000, noise=0.0)

# Plot
plot_simulation(Y, node_names)
```

### One-file pipeline

```python
from dynavista.config import run_from_config

run_from_config("dynavista_config.py")
```

---

## Customising figures

`PlotConfig` controls palette, font size, DPI, and save format for every plot
function in dynavista:

```python
from dynavista import PlotConfig, plot_simulation

cfg = PlotConfig(
    palette="coolwarm",        # matplotlib colormap or list of hex colours
    fontsize=14,
    dpi=300,
    save_path="trajectory.pdf",
    save_format="pdf",
)

plot_simulation(Y, node_names, plot_cfg=cfg)
```

---

## 3. Project Structure

```
dynavista/
├── __init__.py               Clean public API — import everything from here
├── sdds_build.py             Build truth tables from network rules
├── sdds_next_state.py        Single SDDS update step (stochastic + deterministic)
├── sdds_run.py               Single trajectory runner (with/without noise)
├── sdds_sim.py               Ensemble simulation over many initializations
├── truth_table_mutations.py  Node and edge perturbations (mutations, controls)
├── attr_search.py            Attractor identification from simulation endpoints
├── markov.py                 Markov chain analysis (transition matrix, basins)
├── sweep.py                  Automated single-node and combinatorial sweeps
├── parsers.py                Multi-format network file parsers
├── plotting.py               Publication-ready figures
├── config.py                 Config-file pipeline runner (run_from_config)
└── utils.py                  Base conversion utilities (dec2multistate, etc.)

examples/
├── TLGL6/                    6-node T-LGL leukemia model (config + notebooks)
├── PCC_small/                22-node network
├── KLK/                      38-node lung-cancer signaling network
├── TLGL_large/               60-node T-LGL model
└── PCC_large/                69-node network (cyclic attractors)
    # each folder has dynavista_config_<NAME>.py, attractors_<NAME>.csv,
    # (TLGL6 and KLK also ship a worked-example notebook)
```

---

## 4. Core Concepts

### Network rules

Each node's next state is defined by an update rule. The two most common input formats are shown below; see §5 for the full set (polynomial, Boolean word, symbol, and mixed):

**Polynomial form** (GF(p) arithmetic, matches MATLAB):
```python
"x1*x2 + x1 + x2"   # x1 OR x2
"x1 + 1"             # NOT x1
"x1*x2"              # x1 AND x2
"1"                  # always ON
"0"                  # always OFF
```

**Boolean form** (human-readable, auto-converted):
```python
"x1 OR x2"
"NOT x1"
"x1 AND x2"
"True"
"False"
```

You can mix both formats in the same rule list. `sdds_build` auto-detects per rule.

### SDDS structures

`sdds_build()` produces three arrays used by all downstream functions:

| Array | Shape | Content |
|-------|-------|---------|
| `varF` | `(MAXnv, n)` | Input variable indices per node (1-indexed, -1 = unused) |
| `nv` | `(n,)` | Number of inputs per node |
| `F` | `(p^MAXnv, n)` | Truth table values |

### Stochastic update (SDDS)

At each step, each node computes its deterministic target `z[i]`, then:
- **Activation** (`x[i] < z[i]`): node updates with probability `c[0, i]`
- **Degradation** (`x[i] > z[i]`): node updates with probability `c[1, i]`
- **No change** (`x[i] == z[i]`): stays as-is

`c = 0.9 * ones(2, n)` means nodes update 90% of the time when a change is warranted.
`c = 1.0` gives deterministic (synchronous) updates.

### Global noise

`noise > 0` replaces the entire state randomly with probability `noise` per step.
This models stochastic perturbations and prevents convergence to attractors.
**Do not run attractor search on noisy simulations.**

### Simulation output

`Y` has shape `(n, nsteps+1)`. `Y[i, t]` is the fraction of simulation runs in
which node `i` was ON at time step `t`. `Y[:, -1]` gives endpoint expression frequencies.

---

## 5. Module Reference

### `sdds_build.py`

```python
varF, nv, F = sdds_build(rules, node_names, p=2)
```

Converts network rules to SDDS truth table structures. This is the entry
point for every analysis — all other functions consume its output. The rule
format (polynomial, Boolean word, symbol, or mixed) is detected automatically,
so you never need to declare it.

| Parameter | Type | Description |
|-----------|------|-------------|
| `rules` | `list[str]` | One rule per node, polynomial or Boolean (auto-detected) |
| `node_names` | `list[str]` | Ordered node names (index 0 = x1) |
| `p` | `int` | States per node (default 2) |

---

### `sdds_sim.py`

```python
Y, My = simulate(F, varF, nv, p, c, nsteps=20, nins=1000,
                 noise=0.0, n_jobs=1, seed=None)
```

Main simulation entry point. Automatically selects clean or noisy mode.

| Parameter | Description |
|-----------|-------------|
| `F, varF, nv, p` | SDDS structures from `sdds_build` |
| `c` | Propensity matrix, shape `(2, n)` |
| `nsteps` | Time steps per trajectory |
| `nins` | Number of random initializations |
| `noise` | Global noise per step (0 = clean) |
| `n_jobs` | Parallel workers (`-1` = all cores) |
| `seed` | Integer for reproducibility |

**Returns:**
- `Y`: `(n, nsteps+1)` — averaged expression frequencies
- `My`: `list` of per-run state matrices, each `(nsteps+1, n)`

Lower-level functions also available:
```python
sdds_sim(F, varF, nv, p, c, n, nsteps, nins)       # no noise
sdds_sim_noise(g, F, varF, nv, p, c, n, nsteps, nins)  # with noise
sdds_run(x0, F, varF, nv, p, c, nsteps)             # single trajectory
sdds_run_noise(g, x0, F, varF, nv, p, c, nsteps)    # single noisy trajectory
sdds_next_state(x, F, varF, nv, p, c)               # one step
```

---

### `truth_table_mutations.py`

#### `fix_node` — recommended for node mutations and controls

```python
FF = fix_node(F, nv, p, node, v)
```

Locks node `node` (1-indexed) to constant value `v` by overwriting its
F column. Equivalent to `f(i) = v` before `SDDS_Build` in MATLAB.

```python
FF = fix_node(F, nv, p=2, node=3, v=1)   # KRAS = 1 (gain-of-function)
FF = fix_node(F, nv, p=2, node=9, v=0)   # RAF inhibition
FF = fix_node(F, nv, p=2, node=2, v=0)   # LKB1 loss of function
```

#### `fix_nodes` — batch node mutation

```python
FF = fix_nodes(F, nv, p, node_value_pairs=[
    (3, 1),   # KRAS = 1
    (2, 0),   # LKB1 = 0
    (1, 0),   # KEAP1 = 0
])
```

#### `delete_edge` — edge-level structural control

```python
FF = delete_edge(F, nv, varF, p, tail=5, head=7, v=0)
```

Removes the structural influence of node `tail` on `head` by fixing
the tail input to `v` in head's truth table. Use for edge-specific
perturbations where the node itself should remain functional.

#### `delete_node` — structural node removal

```python
FF = delete_node(F, nv, varF, p, node=3, v=0)
```

Removes all outgoing edges from `node` and locks its column to `v`.
More thorough than `fix_node`. Use when modeling complete gene deletion.

#### `apply_mutations` — mixed batch

```python
FF = apply_mutations(F, nv, varF, p, mutations=[
    {"type": "node", "node": 3, "value": 1},
    {"type": "edge", "tail": 5, "head": 7, "value": 0},
])
```

---

### `attr_search.py`

```python
FndAttr, AttrFreq = attr_search(nins, n, Natt, My, Attrs)
```

Matches simulation endpoints to known attractors via Hamming distance on
non-oscillating nodes. Only valid for **noise-free** simulations.

| Parameter | Description |
|-----------|-------------|
| `My` | Output from `simulate()` (noise=0) |
| `Attrs` | `(n, Natt)` matrix. Column = attractor, value -1 = oscillating node |

```python
# Load attractor CSV
Attrs   = load_attractors_csv("attractors/WT_attractors.csv")
Natt    = Attrs.shape[1]
FndAttr, AttrFreq = attr_search(nins, n, Natt, My, Attrs)

# Summarize
summary = summarize_attractors(AttrFreq, nins)
print(summary)
```

**Attractor CSV format** — columns = attractors, rows = nodes in order:
```
0,1,0
1,0,-1
0,-1,1
```
`-1` marks a node that oscillates in that attractor (excluded from Hamming distance).

> **Want to *compute* attractors rather than test against known ones?**
> `attr_search` checks whether simulation endpoints match a set of attractors
> *you already supply* — it does not discover them from the rules. To compute
> attractors/fixed points directly from a network's update rules, see
> [`pystablemotifs`](https://github.com/jcrozum/pystablemotifs) (Rozum, Deritei,
> Park, Gómez Tejeda Zañudo, Albert, *pystablemotifs*, *Bioinformatics*, 2021;
> algorithm: Rozum et al., *Sci. Adv.* 7, eabf8124, 2021), which uses trap-space
> percolation to compute exact attractors efficiently on real biological
> networks.

---

### `sweep.py`

#### Single-node sweep

```python
df = sweep_single_nodes(
    F, varF, nv, p, c,
    node_names=NODE_NAMES,
    phenotype_nodes=[33, 34, 35, 36],
    phenotype_labels=["OxyStress", "PRC2", "CC_Survival", "Cell_Death"],
    nsteps=20, nins=1000, noise=0.0,
    values=[0, 1],            # test knockout AND overexpression
    nodes_to_sweep=None,      # None = all nodes
    n_jobs=-1,                # parallel
    seed=42,
)
```

Returns a DataFrame with one row per condition (wildtype first), columns
for each phenotype endpoint.

#### Delta analysis

```python
df_delta = compute_delta(df, phenotype_labels=["OxyStress","PRC2","CC_Survival"])
```

Adds `delta_<phenotype>` columns showing change from wildtype.

#### Combinatorial sweep

```python
df_combo = sweep_combinations(
    F, varF, nv, p, c,
    node_names=NODE_NAMES,
    phenotype_nodes=[35, 36, 37, 38],
    mutation_sets=[
        [(3,1), (2,0)],   # KRAS + LKB1 double
        [(3,1), (1,0)],   # KRAS + KEAP1 double
        [(3,1), (2,0), (1,0)],  # triple
    ],
    mutation_set_names=["KRAS+LKB1", "KRAS+KEAP1", "Triple"],
    nsteps=20, nins=1000,
)
```

#### Pairwise double-mutant screen

```python
df_pairs = sweep_pairwise(
    F, varF, nv, p, c,
    node_names=NODE_NAMES,
    phenotype_nodes=[35, 36],
    value_a=0, value_b=0,       # all double knockouts
    nsteps=20, nins=300,        # use lower nins for speed
)
```

#### Sweep visualizations

```python
# Heatmap of delta values
plot_sweep_heatmap(df_delta, phenotype_labels, delta=True,
                   save_path="sweep_heatmap.png")

# Ranked bar chart for one phenotype
plot_sweep_ranked(df, "CC_Survival", top_n=15,
                  save_path="ranked_CCsurvival.png")
```

---

### `parsers.py`

Convert published network files into DynaVista format.

```python
# BoolNet (.bnet) — most common published format
rules, names = parse_boolnet("network.bnet")

# BooleanNet / pystablemotifs text format
rules, names = parse_booleannet("network.txt")

# Plain text rules (name = rule, one per line)
rules, names = parse_rules_txt("rules.txt")

# Edge list CSV (source, target, sign columns)
rules, names = parse_adjacency_csv("edges.csv")

# Auto-detect format
rules, names = parse_network("my_file.bnet")

# All parsers feed directly into sdds_build
varF, nv, F = sdds_build(rules, names, p=2)
```

**BoolNet format** (`network.bnet`):
```
targets, factors
NodeA, NodeB & NodeC
NodeB, NodeA | !NodeC
NodeC, 1
```

**BooleanNet format** (`network.txt`):
```
NodeA* = NodeB AND NodeC
NodeB* = NodeA OR NOT NodeC
NodeC* = True
```

**Plain text rules** (`rules.txt`):
```
KEAP1 = 1
LKB1 = 1
KRAS = 0
Methionine = Lat1 OR Homocysteine
```

**Edge list CSV** (`edges.csv`):
```
source,target,sign
A,B,1
C,B,-1
A,C,1
```
Sign values: `1, +1, +, activates` = activation; `-1, -1, -, inhibits` = inhibition.

---

### `markov.py`

#### Exact attractor analysis (networks with n ≤ 20)

```python
A, Avec = bn_asparse(F, varF, nv, p)         # sparse transition matrix
ab, d   = bn_attractor(Avec)                  # attractor labels and distances
attrs   = get_attractor_states(ab, Avec, n, p) # attractor states as binary vectors
bs      = basin_sizes(ab)                      # basin size per attractor
```

`ab` encodes: negative values = on attractor (−1 = attractor 1, −2 = attractor 2, ...),
positive values = in basin of that attractor. `d` = steps to attractor from each state.

#### Full stochastic transition matrix (n ≤ 16 recommended)

```python
A   = multistate_A(F, varF, nv, c, p)   # p^n × p^n transition matrix
G   = google_matrix(A, g=0.9)           # add noise for ergodicity
pi  = stationary_distribution(G)        # stationary distribution
```

#### Scalable stationary distribution (n > 20)

For networks too large to enumerate exactly, estimate the dominant attractors
by sampling. It runs random trajectories, builds a sparse transition matrix
over the visited states only, applies Google damping, and solves for the
stationary distribution of that reduced chain — ranking the dominant attractors
without enumerating 2ⁿ states.

> **Note:** this is an attractor-ranking heuristic, not the exact SDDS
> stationary distribution. It targets the deterministic-plus-damped chain over
> sampled states; its accuracy plateaus with more samples and the mid-rank
> probabilities are not expected to match exact values. Use the exact path
> (`multistate_A` → `stationary_distribution`) when the network is small enough.

```python
result = stationary_distribution_sampled(
    F, varF, nv, p,
    node_names=node_names,
    n_samples=10000,
    nsteps=300,
    g=0.85,
    top_k=10,
)
# result["top_states"]: ranked states with probability and type
#   (fixed point / cycle member); result["n_states_visited"]: distinct states seen
```

> **Computing attractors from the rules.** dynavista's `attr_search` only tests
> whether trajectories end in attractors *you supply* — it does not discover
> them. To compute attractors/fixed points directly from a network's update
> rules, use [`pystablemotifs`](https://github.com/jcrozum/pystablemotifs)
> (Rozum et al., *Bioinformatics* 2021; algorithm: Rozum et al., *Sci. Adv.* 7,
> eabf8124, 2021), which uses trap-space percolation to compute exact
> attractors efficiently on real biological networks.

---

### `plotting.py`

Trajectory plots default to 14pt font and the other plot types to 13pt (all configurable), at 300 DPI with publication-clean axes.

```python
# Apply publication style globally (do this once at the top of your notebook)
set_publication_style()

# Trajectory plot — shows average node expression over time
plot_simulation(Y, node_names,
    nodes_to_plot=[33, 34, 35, 36, 37, 38],
    title="KLK Network — Wildtype",
    save_path="trajectory.png")

# Endpoint bar chart — shows final time-step frequencies
plot_endpoint_bar(Y[:, -1], node_names,
    nodes_to_plot=[33, 34, 35, 36],
    title="Phenotype Readouts")

# Attractor basin pie + bar
plot_attractor_basins(AttrFreq, nins,
    attr_labels=["Attractor 1", "Attractor 2"])

# Stationary distribution top states
plot_stationary_distribution(pi, n, p, top_k=20)

# Sweep heatmap
plot_sweep_heatmap(df_delta, phenotype_labels, delta=True)

# Sweep ranked bar
plot_sweep_ranked(df, "CC_Survival", top_n=15)
```

---

## 6. Config File Guide

A config file defines an entire analysis in one place. It is an ordinary
Python file (`dynavista_config_*.py`) whose module-level constants are read by
`run_from_config()`. Every validated network ships one — see
`examples/TLGL6/dynavista_config_TLGL6.py` for a fully annotated 6-node example,
or `examples/KLK/dynavista_config_KLK.py` for a 38-node one. Copy one, edit the
constants, and run it; the file is heavily commented and is the only thing you
need to edit.

### Minimal config

```python
# my_config.py
NETWORK_NAME = "My Network"
P            = 2                    # 2 = Boolean
OUTPUT_DIR   = "output"

NODE_NAMES = ["A", "B", "C"]
RULES      = ["C", "A", "A OR B"]   # one rule per node, same order as NODE_NAMES

NINS       = 1000                   # random initial conditions
NSTEPS     = 20                     # SDDS time steps
NOISE      = 0.0                    # per-step noise probability
PROPENSITY = 0.9                    # update probability (c value)
SEED       = 42                     # int for reproducibility, or None
```

### Full config with mutations and phenotypes

The example configs group settings into clearly commented blocks. The most
commonly edited ones:

```python
NETWORK_NAME = "KLK Network"
P            = 2
OUTPUT_DIR   = "klk_output"

NODE_NAMES = ["KEAP1", "LKB1", "KRAS", ...]   # 38 total
RULES      = ["1", "1", "0", ...]             # one per node

# Simulation
NINS       = 3000
NSTEPS     = 20
NOISE      = 0.01
PROPENSITY = 0.9
N_JOBS     = -1          # -1 = all CPU cores
SEED       = None

# Mutations: 'exhaustive', 'manual', or 'off'
# In 'manual' mode, each entry is a {node_id: value} dict (node IDs are
# 1-indexed). An empty dict {} is the wildtype (no mutations).
MUTATION_MODE   = 'manual'
MUTATION_MANUAL = [
    {},                  # wildtype (no mutations)
    {3: 1},              # node 3 -> ON
    {3: 1, 2: 0},        # node 3 ON, node 2 OFF
]

# Phenotype reads: (1-indexed node, label)
PHENOTYPE_NODES = [
    (33, "Oxidative Stress"),
    (35, "Cell Survival"),
    (36, "Cell Death"),
    (37, "PDL1"),
]
PLOT_NODES = [33, 34, 35, 36, 37, 38]

# Attractor search against a supplied CSV (None to skip)
ATTRACTOR_CSV = None
```

### Running

```python
from dynavista import run_from_config

results = run_from_config("my_config.py", verbose=True,
                          save_figures=True, show_figures=False)

# Access results (keyed by condition name). The wildtype entry ({}) is keyed
# "wildtype"; other conditions use names generated from their mutations.
# Inspect the keys to see the exact names for your run:
print([k for k in results if not k.startswith("_")])

Y_wt      = results["wildtype"]["Y"]
phenotype = results["wildtype"]["phenotype"]   # dict of endpoint values
```

`run_from_config()` reads every setting from the file, builds the SDDS
structures, simulates each mutation and control condition, optionally runs the
sampled stationary-distribution analysis and the perturbation sweep, saves
figures to `OUTPUT_DIR`, and returns a structured results dict.

---

## 7. Network Format Guide

| Format | Extension | Use case |
|--------|-----------|----------|
| BoolNet | `.bnet` | Published models (GINsim, BioModels, Cell Collective) |
| BooleanNet | `.txt` | pystablemotifs, BooleanNet, BoolSim |
| Plain rules | `.txt` | DynaVista-native, easiest for new networks |
| Adjacency list | `.csv` | STRING, SIGNOR, literature curation |
| Truth table | `.csv` | Precomputed F columns (auto-detected by `parse_network`) |
| Polynomial | inline `list[str]` | Direct from MATLAB SMATA pipeline |

```python
# From any format → sdds_build → simulate
from dynavista import parse_network, sdds_build, simulate

rules, names = parse_network("my_model.bnet")
varF, nv, F = sdds_build(rules, names, p=2)
Y, My = simulate(F, varF, nv, p=2, c=0.9*np.ones((2,len(names))),
                 nsteps=20, nins=1000)
```

---

## 8. Mutation and Control Guide

### When to use each function

| Situation | Function |
|-----------|----------|
| Lock node to 0 or 1 (knock out / overexpress) | `fix_node` |
| Multiple node changes at once | `fix_nodes` |
| Remove one specific edge between two nodes | `delete_edge` |
| Complete structural removal of a node | `delete_node` |
| Mixed node + edge changes | `apply_mutations` |
| In a config file | `MUTATION_MANUAL` / `MUTATION_MODE` |

### Method equivalence

All three methods produce **identical F**:

```python
# Method A — fix_node (recommended, most concise)
FF = fix_node(F, nv, p=2, node=3, v=1)

# Method B — pre-build (most transparent, matches original MATLAB)
rules_mut = rules.copy()
rules_mut[2] = "1"   # index 2 = x3 (0-based)
_, _, FF = sdds_build(rules_mut, names, p=2)

# Method C — in a dynavista_config_*.py file
# MUTATION_MODE   = 'manual'
# MUTATION_MANUAL = [
#     {3: 1},   # node 3 (1-indexed) -> ON
# ]
```

### Common mutation patterns

```python
# NSCLC mutation panel
F_WT        = F.copy()
F_KRAS_mut  = fix_node(F, nv, p=2, node=3, v=1)
F_LKB1_ko   = fix_node(F, nv, p=2, node=2, v=0)
F_KEAP1_ko  = fix_node(F, nv, p=2, node=1, v=0)
F_double    = fix_nodes(F, nv, p=2, node_value_pairs=[(3,1),(2,0)])
F_triple    = fix_nodes(F, nv, p=2, node_value_pairs=[(3,1),(2,0),(1,0)])

# Therapeutic controls
F_RAFi  = fix_node(F, nv, p=2, node=9,  v=0)
F_MEKi  = fix_node(F, nv, p=2, node=15, v=0)
F_combo = fix_nodes(F, nv, p=2, node_value_pairs=[(3,1),(9,0)])  # KRAS mut + RAFi
```

---

## 9. Mutation Sweep Guide

### Screen all single-node knockouts

```python
from dynavista import sweep_single_nodes, compute_delta, plot_sweep_heatmap

df = sweep_single_nodes(
    F, varF, nv, p=2, c=c,
    node_names=NODE_NAMES,
    phenotype_nodes=[35, 36, 37, 38],
    phenotype_labels=["CC_Survival","Cell_Death","PDL1","IFNg"],
    nsteps=20, nins=1000,
    values=[0],           # knockout only
    n_jobs=-1, seed=42,
)

df_delta = compute_delta(df, ["CC_Survival","Cell_Death","PDL1","IFNg"])

# Heatmap of all knockouts vs phenotype
plot_sweep_heatmap(df_delta, ["CC_Survival","Cell_Death","PDL1","IFNg"],
                   title="Single-Node Knockout Screen",
                   save_path="knockout_screen.png")

# Which knockout most increases CC_Survival?
plot_sweep_ranked(df, "CC_Survival", target_value=0, top_n=15,
                  save_path="ranked_survival.png")
```

### Save sweep results

```python
df.to_csv("sweep_results.csv", index=False)
df_delta.to_csv("sweep_delta.csv", index=False)
```

---

## 10. Markov Chain Guide

### Small networks (n ≤ 20) — exact analysis

```python
from dynavista import bn_asparse, bn_attractor, get_attractor_states, basin_sizes

A, Avec = bn_asparse(F, varF, nv, p=2)
ab, d   = bn_attractor(Avec)

attrs = get_attractor_states(ab, Avec, n, p=2)
for k, states in attrs.items():
    print(f"Attractor {k}: {[list(s) for s in states]}")

bs = basin_sizes(ab)
for k, size in bs.items():
    print(f"Attractor {k} basin: {size} states ({100*size/2**n:.1f}%)")
```

### Small networks — stationary distribution

```python
from dynavista import multistate_A, google_matrix, stationary_distribution

A   = multistate_A(F, varF, nv, c, p=2)   # exact transition matrix
G   = google_matrix(A, g=0.9)             # regularize for Perron-Frobenius
pi  = stationary_distribution(G)          # stationary probabilities

# Top 10 states by stationary weight
top_idx = pi.argsort()[::-1][:10]
for i in top_idx:
    state = dec2multistate(i, p=2, n=n)
    print(f"State {i}: {list(state)}  pi={pi[i]:.4f}")
```

### Large networks — sampling

```python
from dynavista import stationary_distribution_sampled

result = stationary_distribution_sampled(
    F, varF, nv, p=2, node_names=node_names,
    n_samples=10000, nsteps=300, g=0.85, top_k=5,
)
for s in result["top_states"]:
    on = [k for k, v in s["state_labels"].items() if v == 1]
    print(f"rank {s['rank']}: p={s['probability']:.3f} ({s['type']}) ON: {on}")
```

---

## 11. Scalability Guide

| Network size | Recommended approach | Notes |
|---|---|---|
| n ≤ 12 | Full Markov chain (`multistate_A`) | Dense p^n × p^n matrix; fits comfortably in memory |
| n ≤ 38 | Simulation (`simulate`) | e.g. KLK (n=38): ~5s for 3000 runs × 20 steps |
| n ≤ 69 | Simulation (`simulate`) | e.g. PCC_large (n=69): ~8s for 3000 runs × 20 steps; ~35s for 1000 runs × 300 steps |
| n > 69 | `simulate` with `n_jobs=-1` | Parallelize across all CPU cores |
| Large n, long-run behaviour | `stationary_distribution_sampled` | Sampling-based attractor ranking (see §10) |

### Speeding up large simulations

```python
# Use all available CPU cores
Y, My = simulate(F, varF, nv, p=2, c=c,
                 nsteps=20, nins=5000,
                 n_jobs=-1)

# Reduce nins for sweeps (500 is often sufficient for ranking)
df = sweep_single_nodes(..., nins=500, n_jobs=-1)
```

---

## 12. MATLAB Equivalence Table

| MATLAB | Python (DynaVista) |
|--------|----------------|
| `SDDS_Build(syms, f, p)` | `sdds_build(rules, names, p)` |
| `SDDS_sim(F, varF, nv, p, c, n, nsteps, nins)` | `simulate(F, varF, nv, p, c, nsteps, nins)` |
| `SDDS_simNoise(g, F, varF, nv, p, c, n, nsteps, nins)` | `simulate(F, varF, nv, p, c, nsteps, nins, noise=g)` |
| `SDDSRun(x, F, varF, nv, p, c, nsteps)` | `sdds_run(x, F, varF, nv, p, c, nsteps)` |
| `SDDSRunNoise(g, x, F, ...)` | `sdds_run_noise(g, x, F, ...)` |
| `SDDSNextState(x, F, varF, nv, p, c)` | `sdds_next_state(x, F, varF, nv, p, c)` |
| `f(i) = v` before Build | `fix_node(F, nv, p, i, v)` |
| Multiple `f(i) = v` | `fix_nodes(F, nv, p, [(i,v), ...])` |
| `TruthTable_del_n_temp(...)` | `delete_node(F, nv, varF, p, node, v)` |
| `TruthTable_del_a_temp(...)` | `delete_edge(F, nv, varF, p, tail, head, v)` |
| `AttrSearch(nins, n, Natt, My, Attrs)` | `attr_search(nins, n, Natt, My, Attrs)` |
| `bnAsparse(F, varF, nv)` | `bn_asparse(F, varF, nv, p)` |
| `bnAttractor(Avec)` | `bn_attractor(Avec)` |
| `multistateA(F, varF, nv, c, p)` | `multistate_A(F, varF, nv, c, p)` |
| `dec2multistate(y, p, n)` | `dec2multistate(y, p, n)` |
| `multistate2dec(x, p, n)` | `multistate2dec(x, p, n)` |
| `dtmc(G); asymptotics(mc)` | `stationary_distribution(G)` |
| `figure; plot(X, Y(i,:), ...)` | `plot_simulation(Y, names, nodes_to_plot=[i])` |
| (no equivalent) | `run_from_config("dynavista_config.py")` |
| (no equivalent) | `sweep_single_nodes(...)` |
| (no equivalent) | `parse_boolnet(...)` |

---

## 13. Citation

If you use this software in your research, please cite the accompanying
publication. Citation details (BibTeX and full reference) will be added here
once the manuscript is published.

<!-- TODO: add published citation (BibTeX + reference) at publication time. -->

---

## License

MIT License. See LICENSE file.
