par2-circadian

AR(2) eigenvalue analysis for gene expression time series

Fits second-order autoregressive models to gene expression data and computes the eigenvalue modulus |λ|, a single number that quantifies how strongly a gene's past determines its future (temporal persistence). Discovers the three-layer hierarchy: Clock > Target > Background.

Installation

Requires Python 3.9 or later.

pip install par2-circadian

Or install from source:

git clone https://github.com/mickwh2764/par2discovery.git
cd par2discovery
pip install .

Quick Start

Python API

import par2

# Single gene
result = par2.fit_ar2([1.2, 3.4, 2.1, 4.5, 3.2, 5.1, 2.8, 4.9, 3.5, 5.2, 2.9, 4.7])
print(f"|λ| = {result['eigenvalue']:.3f}, type = {result['root_type']}")

# Whole matrix from CSV
matrix, genes = par2.load_expression_matrix("my_data.csv")
results = par2.fit_ar2_batch(matrix, genes)

# Discover the three-layer hierarchy
hierarchy = par2.discover_hierarchy(results)
print(f"Clock median:  {hierarchy['clock_median']:.3f}")
print(f"Target median: {hierarchy['target_median']:.3f}")
print(f"Gearbox gap:   {hierarchy['gearbox_gap']:.3f}")
print(f"Health grade:  {hierarchy['health_grade']}")

# Save results
par2.save_results(results, "ar2_results.csv")

Confidence intervals

A point estimate of |λ| from a short series is imprecise, so fit_ar2 can return a residual-bootstrap 95% interval alongside it:

result = par2.fit_ar2(expression, n_bootstrap=2000, seed=1)
print(result["eigenvalue"], result["eigenvalue_ci"])

The residuals are resampled, the series regenerated from the fitted recursion and the model refitted, so only the one observed series is needed — no biological replicates, which most circadian designs do not provide. Intervals are opt-in and absent unless n_bootstrap > 0.

Command Line

# Analyse a CSV file (genes as rows, timepoints as columns)
par2 my_data.csv -o results.csv

# Show top 20 genes by eigenvalue
par2 my_data.csv --top 20

Example Dataset

An example dataset is included at data/example_circadian.csv (30 genes x 12 timepoints) with known clock, target, and background genes showing realistic circadian dynamics:

import par2

matrix, genes = par2.load_expression_matrix("data/example_circadian.csv")
results = par2.fit_ar2_batch(matrix, genes)
h = par2.discover_hierarchy(results)
print(f"Hierarchy preserved: {h['hierarchy_preserved']}")  # True
print(f"Health grade: {h['health_grade']}")  # A

Input Format

CSV file with:

• First row: header (timepoint labels)
• First column: gene names
• Remaining columns: expression values (minimum 6 timepoints)

Example:

Gene,ZT0,ZT2,ZT4,ZT6,ZT8,ZT10,ZT12,ZT14,ZT16,ZT18,ZT20,ZT22
Bmal1,12.3,14.1,15.8,16.2,14.5,11.2,9.8,8.1,7.5,8.9,10.2,11.5
Per2,8.1,7.2,6.5,7.8,10.2,13.5,15.1,14.8,13.2,10.5,9.1,8.5

Output

Per-Gene Results

Each gene gets:

• eigenvalue: |λ|, the eigenvalue modulus (0 to ~1). Higher = more persistent.
• phi1, phi2: AR(2) coefficients
• r2: goodness of fit
• root_type: 'Complex' (oscillatory) or 'Real' (monotone decay)
• half_life: persistence half-life in sampling intervals
• eigenperiod: intrinsic oscillation period (complex roots only)

Hierarchy Discovery

discover_hierarchy() returns:

• clock_median / target_median / background_median: layer-wise eigenvalue medians
• gearbox_gap: clock_median − target_median (the circadian health metric)
• hierarchy_preserved: True if clock > target > background
• health_grade: A (gap ≥ 0.15) through F (gap < 0.02)
• clock_genes / target_genes: per-gene eigenvalue lists

Interpreting |λ|

These bands describe groups of genes, not individual ones. At 24 evenly sampled timepoints the 95% interval on a single gene's |λ| is typically ~0.4 wide, and the estimator is biased upward below ~24 points, so a lone value can easily fall in the wrong row. Report eigenvalue_ci with any per-gene number, compare groups rather than genes, and do not interpret differences smaller than the interval.

Method

The AR(2) model fits:

x(t) = φ₁·x(t-1) + φ₂·x(t-2) + ε

The characteristic equation r² − φ₁r − φ₂ = 0 yields eigenvalues whose modulus |λ| quantifies temporal persistence. Expression values are mean-centred before fitting.

For complex roots: |λ| = √(−φ₂) For real roots: |λ| = max(|r₁|, |r₂|)

The three-layer hierarchy emerges because clock genes (strong autonomous oscillation) have higher |λ| than clock-controlled target genes (driven oscillation), which in turn have higher |λ| than background genes (no circadian regulation).

See: Whiteside M (2026). "AR(2) eigenvalue modulus as a measure of temporal persistence in gene expression: circadian hierarchy emerges from two coefficients." Research Square [Preprint]. doi:10.21203/rs.3.rs-9283100/v1

Benchmark and scope

benchmarks/ holds a ground-truth simulation suite comparing |λ| with JTK_CYCLE, ARSER, Lomb–Scargle (via MetaCycle) and RAIN on series whose answer is known. In summary:

• |λ| is the most sensitive of the methods tested on noise-driven AR(2) oscillators (92% vs 40–62% at true |λ| = 0.618, n = 24, matched 5% false-positive rate) and the only one that detects non-24 h rhythms without being told the period.
• The p-value methods beat it on clean deterministic waveforms (100% vs 83% on a sustained cosine).
• |λ| is a persistence measure, not a rhythm test. Calibrated on white noise it calls plain AR(1) red noise at 85%, because red noise genuinely is persistent. Pair |λ| with a rhythm test before calling a gene an oscillator.
• Point estimates are biased upward at low persistence (+0.21 at true 0.30, n = 24), which compresses rather than manufactures a clock-versus-target gap.

See benchmarks/README.md for the full tables and reproduction steps.

Researcher Profile

Michael Whiteside — independent computational systems researcher
Researcher Profile · ORCID 0009-0000-0643-5791 · par2discovery.com

AI Usage Statement

An AI coding agent (Replit Agent, Anthropic/Replit) was used to assist with implementation of platform modules and translation of analytical specifications into code. Each module was specified by the author in terms of mathematical requirements and validation criteria, implemented by the agent, and verified against known analytical results before incorporation.

Large language models (Claude, Anthropic; GPT-4, OpenAI) were used for drafting, structural editing, and refinement of documentation and manuscript text. All scientific content — the analytical framework, hypothesis design, dataset selection, validation architecture, and scientific interpretation — originated with and was decided by the author.

The AR(2) framework, the eigenvalue hierarchy hypothesis, and all scientific judgements are the author's own. The responsibility for the correctness and adequacy of all results remains exclusively with the human author.

License

Apache License 2.0 — free for any use, including commercial, see LICENSE. Releases up to and including 1.1.8 were published under PolyForm Noncommercial 1.0.0; from 1.2.0 the package is open source, so that it can be reviewed and depended on by the venues and pipelines that require an OSI-approved licence.

The PAR(2) methodology is the subject of pending UK patent application GB2518973.9. Section 3 of Apache-2.0 grants you a patent licence covering use of this software; it does not grant rights to practise the methodology by other means (see NOTICE). For that, contact mickwh@msn.com.

Citation

Machine-readable metadata is in CITATION.cff — GitHub's "Cite this repository" button reads it, as do Zenodo, cffconvert and most reference managers. Please cite the software and the method preprint together:

Whiteside M (2026). par2-circadian: AR(2) eigenvalue analysis for gene expression time series. Version 1.2.0. Zenodo. doi:10.5281/zenodo.21963192

Whiteside M (2026). "AR(2) eigenvalue modulus as a measure of temporal persistence in gene expression: circadian hierarchy emerges from two coefficients." Research Square [Preprint]. doi:10.21203/rs.3.rs-9283100/v1

Every release is archived on Zenodo. 10.5281/zenodo.21963192 is the concept DOI and always resolves to the latest version; cite a version DOI only when you need to pin exactly what you ran.

In a Methods section, identify the software as par2-circadian (RRID:SCR_028837). The package is also registered in bio.tools as biotools:par2-circadian.

The full publication list lives at https://par2discovery.com/publications.bib — one entry per work, every DOI verified against Crossref or DataCite.