Advanced Features

This page documents optional input files and parameters that activate advanced modes of CardamomOT. All optional files are placed in my_project/Data/ (or set directly in adata.obs).

Project directory structure (full)

my_project/
├── Data/
│   ├── data.h5ad                    # required
│   ├── gene_list.txt                # optional — restrict to a gene subset
│   ├── stimulus_schedule.txt        # optional — custom stimulus values per timepoint
│   ├── stimulus_schedule_simul.txt  # optional — overrides schedule for simulation only
│   ├── ref_network.csv              # optional — prior interaction graph (sparsity mask)
│   ├── basal_init.npy / .csv        # optional — warm-start for basal parameters
│   ├── inter_init.npy / .csv        # optional — warm-start for interactions
│   ├── basal_ref.npy / .csv         # optional — regularisation target for basal
│   ├── inter_ref.npy / .csv         # optional — regularisation target for interactions
│   ├── KO_OV_inference.txt          # optional — per-sample KO/OV prior (requires dataset_id)
│   ├── KO_OV_simulate.txt           # optional — in-silico perturbations to simulate
│   └── transition_rates.csv         # optional — cell-type transition cost matrix for OT
└── cardamomOT/                      # generated by the pipeline
    ├── basal.npy                    # (n_samples, G, n_networks)
    ├── inter.npy                    # (G, G, n_networks)
    └── basal_ref_mask.npy           # (n_samples, G) bool — from KO_OV_inference.txt

---

Gene subset (Data/gene_list.txt)

By default select_DEgenes_and_split.py selects differentially expressed genes automatically. To fix the gene set manually, place a plain-text file with one gene name per line:

CHGA
STMN2
PHOX2B
HAND2
...

cardamomot step select_DEgenes_and_split -i my_project -s full -c 0 --mean-forcing 0.5 --force-basins 1.0 --temporal-basins 1

When gene_list.txt is present, those genes are added to the final selection — they are not a replacement. The full behaviour depends on the -c flag:

• -c 1: genes variable across time, conditions, and cell types are selected automatically, then the genes in gene_list.txt are appended (provided they are non-constant in time).

• -c 0: no automatic differential selection is performed; only the genes in gene_list.txt that are non-constant in time are kept.

---

Read-depth correction (infer_rd)

When cells have heterogeneous sequencing depths, run this optional step before gene selection:

cardamomot step infer_rd -i my_project

This identifies Poisson-like genes (Chronocell method, Fang et al. 2024) and stores a per-cell read-depth factor in adata.obs['rd']. Subsequent steps (infer_mixture, infer_network_structure) use it automatically when the column is present.

Key parameters: --var_threshold 1.2 (Poisson variance threshold), --min_mean 0.1.

Requires at least 1 000 genes. Below this threshold the step is silently skipped.

---

Prior interaction network (prepare_reference_network)

To bias network inference toward known biology, first build a prior network from public databases (OmniPath, etc.):

cardamomot step prepare_reference_network -i my_project -d 3

This generates Data/ref_network.csv, a gene × gene binary (or weighted) interaction matrix. The --prior parameter then controls how strongly edges absent from this prior are penalised.

--prior value	Effect
1.0 (default)	No penalisation — all edges equally possible
0.5	Soft penalisation — absent edges discouraged
0.0	Hard mask — only edges present in the prior are allowed

--prior must be passed to both steps — it constrains inference in infer_network_structure and the simulation reference network in infer_network_simul. Always use the same value in both.

# Easiest: let the pipeline handle it (prior network with default depth 3)
cardamomot pipeline -i my_project -s full --ref --prior 0.5

# With a shallower graph (faster):
cardamomot pipeline -i my_project -s full --ref --ref-depth 1 --prior 0.5

# Or step by step — pass --prior to both:
cardamomot step infer_network_structure -i my_project -s full --prior 0.5
cardamomot step infer_network_simul     -i my_project -s full --prior 0.5

---

Custom stimulus schedule (Data/stimulus_schedule.txt)

By default one stimulus is 0 at the first timepoint and 1 at all subsequent ones. To define a custom schedule, create a tab/space-separated file where rows = timepoints and columns = stimulus channels:

# stimulus_schedule.txt
0.0    0.0
1.0    0.0
1.0    1.0

• Values between 0 and 1 are allowed (partial stimulus strength).

• If fewer rows than timepoints are provided, the last row is repeated.

• To use a different schedule for simulation (e.g. a novel protocol), place Data/stimulus_schedule_simul.txt; it takes priority over stimulus_schedule.txt during simulate_network and simulate_network_KOV.

The --stimulus parameter controls how strongly the stimulus regulates genes in the inferred network:

--stimulus value	Effect
1.0 (default)	Full stimulus influence
0.0	Stimulus has no regulatory influence

---

Multiple experimental samples (dataset_id)

If your experiment contains several biological conditions that share a gene regulatory network but have different basal transcription rates (e.g. different cell lines, donors, or perturbation backgrounds):

adata.obs['dataset_id'] = ...   # string or integer label per cell

CardamomOT then:

• Builds one set of per-sample basal parameters θ_basal(s) for each sample.

• Runs a single joint optimisation over shared interaction weights and all per-sample basals simultaneously.

• Saves basal.npy with shape (n_samples, G, n_networks).

Keeping per-sample basals close (constrain_basal_uniform): set model.constrain_basal_uniform = λ (e.g. λ = 100–1000) to add an L2 penalty pushing each sample’s basals toward their common mean. Genes with a non-zero basal_ref entry for a given sample (e.g. a KO prior) are excluded from the penalty for that sample.

---

Regularisation targets (basal_ref, inter_ref)

Place reference arrays in Data/ to anchor the optimiser:

File	Shape	Role
basal_ref.npy / .csv	(G,) or (n_samples, G, n_networks)	Regularisation target for basal parameters
inter_ref.npy / .csv	(G, G) or (G, G, n_networks)	Regularisation target for interactions
basal_init.npy / .csv	same as basal_ref	Warm-start values (not regularised)
inter_init.npy / .csv	same as inter_ref	Warm-start values for interactions

CSV files must have gene names matching adata.var_names as row and column indices.

---

KO/OV inference with per-sample priors (Data/KO_OV_inference.txt)

When dataset_id is present, you can encode prior knowledge about which genes are knocked out or overexpressed in each sample:

# KO_OV_inference.txt  (tab-separated, header required)
sample_id           KO      OV
wt                  0       0
ko_CHGA             CHGA    0
ov_STMN2            0       STMN2
ko_CHGA_ov_STMN2    CHGA    STMN2

• sample_id must match adata.obs['dataset_id'].

• KO genes are forced to basal = −100 (silent during inference).

• OV genes are forced to basal = +100 (always active).

• Multiple genes per cell: CHGA,POSTN.

During inference, this replaces basal_ref for the affected genes/samples and saves cardamomOT/basal_ref_mask.npy for use by simulate_network_KOV.

cardamomot step infer_network_structure -i my_project -s full
# KO_OV_inference.txt is read automatically when present

---

In-silico perturbation simulation (Data/KO_OV_simulate.txt)

After training, simulate arbitrary KO/OV combinations by defining them in:

# KO_OV_simulate.txt  (tab-separated, header required)
KO          OV
none        none
CHGA        STMN2
POSTN       S100B,STMN2
CHGA-80     STMN2-60

Each row produces one independent simulation saved as cardamomOT/adata_sim_KO_*_OV_*_stim*.h5ad.

Partial KO/OV (GENE-X syntax): append -X (0 < X < 100) for a partial perturbation of strength X%:

Mode	Creation rate factor	Effect
KO X%	× (1 − X/100)	Reduced production → lower steady-state
OV X%	× 1/(1 − X/100)	Increased production → higher steady-state

Output files use the pctX suffix to distinguish partial from complete perturbations (KO_CHGApct80_OV_STMN2pct60).

cardamomot step simulate_network_KOV -i my_project -s full

Visualise results with:

from CardamomOT import plot_results_sim_kov, compare_cell_types

p, stim, prior = "my_project/", 1.0, 1.0
compare_cell_types(p, "KO_CHGA_OV_none", split="full", stim=stim, prior=prior)
plot_results_sim_kov(p, "KO_CHGA_OV_none", stim=stim, prior=prior)

---

Population dynamics: proliferation, death, and transition rates

Per-cell proliferation/death (adata.obs)

adata.obs['prolif_rate'] = ...   # net proliferation rate per cell
adata.obs['death_rate']  = ...   # net death rate per cell

When both columns are present, the OT marginals are corrected so that faster-growing cells carry proportionally more weight as trajectory sources. This implements demographically corrected OT (as in Waddington OT, Schiebinger et al. 2019).

Cell-type transition rates (Data/transition_rates.csv)

To bias the OT cost toward biologically plausible transitions:

# transition_rates.csv — rows = source type at t1, cols = target type at t2
         ,TypeA,TypeB,TypeC
TypeA    ,  0.3,  0.1,  0.01
TypeB    ,  0.05, 0.2,  0.05
TypeC    ,  0.01, 0.05, 0.3

Row/column names must match adata.obs['cell_type']. At each consecutive timepoint pair, transition probabilities are computed as exp(rate × Δt) and rescaled so that the mean weight per row equals 1 — transitions with weight > 1 become cheaper (preferred), weight < 1 become more expensive (penalised).

---

Proliferation-aware simulation (--proliferation)

CardamomOT can learn and simulate net proliferation rates (R = birth − death) directly from the inferred optimal-transport couplings, without requiring external annotations.

How it works

Learning R from the OT coupling.
During trajectory inference (infer_network_structure), at each consecutive timepoint pair (t → t+1), the row marginals of the optimal coupling encode the effective mass transported out of each source cell. The optimal net growth rate per cell is recovered via the self-consistent zero-mean inversion of the WOT marginal formula:

m_n  =  (1/N) · exp(R_n · Δt/2) / C          (WOT: source marginal, normalised)

R_n  =  (2/Δt) · (log(m_n) − mean_k(log(m_k)))

The factor 2/Δt comes from the WOT convention where both the source and the target marginals are corrected by exp(R · Δt/2) each — the formula recovers the true net rate R from the source marginal alone. With balanced OT (default) all m_n = 1/N so R_n = 0 exactly; with unbalanced OT (model.unbalanced_reg > 0) they vary and encode genuine differential growth. The resulting R_opt array (one value per cell per timepoint transition) is always saved to cardamomOT/data_R_opt.npy.

Training the ProliferationMLP.
infer_network_simul --proliferation reads data_R_opt.npy and trains a lightweight two-hidden-layer MLP (64 units, Tanh activation, MSE loss):

R(P) : protein levels → net proliferation rate

The network is trained on (prot[:, ns:], R_opt) pairs after any recomputation of protein trajectories, so it captures the optimal R consistent with the final inferred network. The weights are saved to cardamomOT/prolif_network.pt.

Branching simulation.
simulate_network --proliferation (and simulate_network_KOV --proliferation) loads the MLP and applies a trapezoidal forward+backward resampling after each simulated interval — mirroring the structure of the WOT correction used during inference (which corrects both source and target marginals):

log_weight_n  =  (R(P_start_n) + R(P_end_n)) / 2 · Δt

P_start and P_end are the protein states of cell n at the beginning and end of the interval. The N cells are then resampled multinomially (with replacement, weighted by exp(log_weight_n)). This avoids the catastrophic variance of Poisson sampling (which would kill ~37 % of cells per step even at R = 0) and is consistent with the inference-time coupling.

Pipeline usage

# Step 4 — adapt parameters and learn R
cardamomot step infer_network_simul -i my_project -s full --proliferation

# Step 5 — simulate with branching
cardamomot step simulate_network -i my_project -s full --proliferation

# KO/OV perturbations with branching
cardamomot step simulate_network_KOV -i my_project -s full --proliferation

Or directly:

python infer_network_simul.py -i my_project -s full --proliferation
python simulate_network.py    -i my_project -s full --proliferation

Notes

• --proliferation is off by default. infer_network_structure always computes and saves R_opt as a by-product regardless — enabling the flag later has no cost.

• The flag must be passed to both infer_network_simul (trains the MLP) and simulate_network / simulate_network_KOV (uses it). Passing it only to the simulation scripts while prolif_network.pt is absent prints a warning and falls back to standard simulation.

• Resampling is currently applied to the PDMP stochastic mode only; ODE mode ignores it silently.

• For balanced OT the estimated R values are exactly zero for all cells; meaningful proliferation signals require unbalanced OT (model.unbalanced_reg > 0) or prior annotations via adata.obs['prolif_rate'] / adata.obs['death_rate'] (see section above).

• The trapezoidal rule (R_start + R_end)/2 naturally handles cells whose protein levels — and thus growth rate — change significantly within the interval.