Circadian Workbench analysis explained

For a methods section

Activity recordings were imported from ClockLab AWD or delimited timestamp/elapsed-time tables; timestamps were sorted and deduplicated, negative activity sentinels were treated as missing, and the positive median sampling interval was inferred when it was not declared. The selected channel was restricted to the inclusive analysis interval and inclusive exclusions, then count activity was summed into 15-minute bins while continuous channels were averaged. Unless changed in the saved normalized configuration, analyses used a 24-hour reference cycle, a period search from 18 to 30 hours, a 0.05-hour period step, Lomb–Scargle as the primary period estimator, and a folded 24-hour activity profile. Daily activity onset and offset used the 20th percentile of positive binned activity with six-hour inactive and active templates; accepted marker times were circularly unwrapped and regressed on day to estimate tau. A fixed-period cosinor estimated the midline estimating statistic of rhythm, amplitude, acrophase, coefficient of determination, F statistic, and p value; non-parametric analysis estimated L5, M10, relative amplitude, interdaily stability, and intradaily variability. Question-specific analyses used their declared recording, point, trace, stimulus, schedule, temperature, effect-size, or fitted-model inputs and the normalized settings stored with the result. Missing, under-covered, non-identifiable, or scientifically undefined quantities were withheld rather than replaced by zero. Result JSON and exports recorded the software version, source SHA-256 fingerprint, normalized settings, inputs, and generation time.

Overview

Think of the workbench as a railway junction: import and selection create one checked track, then independent analyses branch from that same track and reunite only in the result and provenance record.

file or declared points/traces/model
        |
        v
validate input -> normalize settings -> select valid time/channel
        |                                  |
        |                                  +-> daily markers and onset fits
        |                                  +-> profile and fixed cosinor
        |                                  +-> spectral period estimates
        |                                  +-> non-parametric rhythm measures
        |                                  +-> rhythmicity or Hilbert phase
        |
        +-> direct-input science ---------> derived analysis or prediction
                                             |
                                             v
                              result + warnings + provenance + replay script

The branches are independent where the scientific question is independent. For example, temperature compensation consumes already estimated periods, while phase angle consumes daily markers and a dated zeitgeber schedule.

The plot gallery is the visual companion to this execution map. These representative panels are produced by the same renderers used by the package's publication export:

Compare period estimatesSee what unlike estimators can and cannot report. Track period driftSee a time-local result instead of one global period. Separate damping from desynchronisationSee population coherence and amplitude together. Fit a phase-response curveSee animal points and the periodic fitted response.

The analysis in order

Step 1 — Import, select, and bin

ClockLab declares its interval in the header; timestamped and elapsed tables use the median positive time difference. The source is sorted and duplicate timestamps keep the final row. Missing activity stays missing. Selection and exclusions are inclusive. For a bin (j), counts are summed and continuous levels are averaged over the (n_j) observed samples.

(1) count_bin[j] = SUM(x[i]); level_bin[j] = SUM(x[i]) / n[j]

$$ \tag{1} C_j=\sum_{i\in j}x_i,\qquad L_j=\frac{1}{n_j}\sum_{i\in j}x_i $$

Engine: io.py importers and analysis.py:1147-1299 selection/resampling.

Step 2 — Fold the activity profile

Each selected bin is assigned a phase modulo reference period (P). Mean, sample standard deviation, observation count, and standard error are calculated separately at every phase bin.

(2) phase[i] = time[i] MOD P; SEM[j] = SD[j] / SQRT(n[j])

$$ \tag{2} \phi_i=t_i\bmod P,\qquad SEM_j=\frac{SD_j}{\sqrt{n_j}} $$

Engine: analysis.py:1719-1780.

Step 3 — Detect daily markers and estimate onset tau

The onset template rewards activity after a candidate and inactivity before it. A flat score curve is withheld. Accepted circular marker times are unwrapped and fitted against day (d); one hour of daily drift changes tau by one hour.

(3) score(c) = (SUM(after_active) - SUM(before_inactive)) / (n_on + n_off); onset[d] = a + b*d; tau = 24 + b

$$ \tag{3} S(c)=\frac{\sum x_{after}-\sum x_{before}}{n_{on}+n_{off}},\qquad u_d=a+bd,\qquad \tau=24+b $$

Engine: analysis.py:1453-1716.

Step 4 — Fit the fixed-period cosinor

Ordinary least squares fits an intercept and sine/cosine pair at angular frequency ω. Amplitude is vector length and acrophase is its circular direction. The model-versus-residual variance ratio gives the F statistic.

(4) y = M + beta*COS(omega*t) + gamma*SIN(omega*t); A = SQRT(beta^2 + gamma^2); phase = ATAN2(gamma,beta)/omega MOD P

$$ \tag{4} y_i=M+\beta\cos(\omega t_i)+\gamma\sin(\omega t_i),\quad A=\sqrt{\beta^2+\gamma^2},\quad \phi=\frac{\operatorname{atan2}(\gamma,\beta)}{\omega}\bmod P $$

Engine: analysis.py:1980-2011. The source contains no scholarly citation for this implementation.

Step 5 — Calculate non-parametric circadian rhythm analysis

Intradaily variability (IV) uses successive differences, interdaily stability (IS) compares hourly means with total variance, and relative amplitude (RA) contrasts the lowest five-hour and highest ten-hour circular windows.

(5) IV = N*SUM((x[i]-x[i-1])^2)/(A*SUM((x[i]-mean)^2)); IS = N*SUM((hour_mean[h]-mean)^2)/(24*SUM((x[i]-mean)^2)); RA = (M10-L5)/(M10+L5)

$$ \tag{5} IV=\frac{N\sum(x_i-x_{i-1})^2}{A\sum(x_i-\bar x)^2},\quad IS=\frac{N\sum_{h=0}^{23}(\bar x_h-\bar x)^2}{24\sum(x_i-\bar x)^2},\quad RA=\frac{M10-L5}{M10+L5} $$

Engine: analysis.py:2496-2551. The source contains no scholarly citation for this implementation.

Step 6 — Estimate and compare periods

Every estimator searches the same declared period range after its stated detrending and resampling. The reported candidate maximizes that estimator's score; comparison preserves failures rather than substituting the primary method.

(6) period_hat[m] = ARGMAX over P in search_grid of score[m](P)

$$ \tag{6} \widehat P_m=\underset{P\in\mathcal P}{\operatorname{argmax}}\;Q_m(P) $$

Engine: analysis.py:2014-2298, period_methods.py:394-1983. Implemented estimators include Lomb–Scargle, Sokolove–Bushell chi-square, Dörrscheidt–Beck F, FFT-NLLS, Burg maximum entropy, mFourFit, and spectrum resampling.

Step 7 — Test rank-based rhythmicity

JTK_CYCLE compares observed and template orderings with Kendall tau-b. Empirical JTK_CYCLE permutes the observed ordering; the plus-one correction prevents a zero p value.

(7) tau_b = (concordant-discordant)/SQRT((pairs-ties_x)*(pairs-ties_y)); p_emp = (1 + COUNT(abs(tau_null)>=abs(tau_obs))) / (1 + permutations)

$$ \tag{7} \tau_b=\frac{C-D}{\sqrt{(N_p-T_x)(N_p-T_y)}},\qquad p_{emp}=\frac{1+#(|\tau_0|\ge|\tau_{obs}|)}{1+N_{perm}} $$

Engine: rhythm_detection.py:54-532. Dense data are reduced to at most 96 points by default.

Step 8 — Track instantaneous phase

A zero-phase band-pass filter isolates the configured circadian band. The Hilbert transform forms analytic signal (z(t)); unwrapped phase slope gives local period and magnitude gives local amplitude. Edge cycles are discarded.

(8) z(t) = filtered_x(t) + i*HILBERT(filtered_x(t)); period(t) = 2*pi / DERIVATIVE(unwrapped_phase(t)); amplitude(t) = ABS(z(t))

$$ \tag{8} z(t)=x_f(t)+i\mathcal H[x_f(t)],\qquad P(t)=\frac{2\pi}{d\theta(t)/dt},\qquad A(t)=|z(t)| $$

Engine: instantaneous.py:134-532. This project-specific interpretation guard has no citation in source.

Step 9 — Fit a phase-response curve

Pre-scored shifts or record-derived onset shifts are placed on circadian time. One or two periodic harmonics are fitted by least squares. Advance is positive; insufficient valid points returns no fit.

(9) CT = (CT_onset + (stimulus-projected_onset)*24/tau) MOD 24; shift(CT) = a0 + SUM(a[k]*COS(2*pi*k*CT/24) + b[k]*SIN(2*pi*k*CT/24))

$$ \tag{9} CT=\left[CT_{on}+\frac{24(t_s-\widehat t_{on})}{\tau}\right]_{24},\qquad \Delta(CT)=a_0+\sum_k\left[a_k\cos\frac{2\pi kCT}{24}+b_k\sin\frac{2\pi kCT}{24}\right] $$

Engine: phase_response.py:65-468, orchestrated by actions.py:1951-2000. The source names the Aschoff sign convention and Winfree curve definition without a full citation.

Step 10 — Measure phase angle

Each daily marker is differenced from the selected zeitgeber transition and wrapped to ±12 hours. Unit vectors give circular mean and resultant length (R); Rayleigh concentration plus absence of material drift supports the entrainment verdict.

(10) psi[i] = WRAP24(marker[i]-zeitgeber[i]); R = ABS(SUM(EXP(i*2*pi*psi[i]/24))/n); circular_SD = SQRT(-2*LN(R))*24/(2*pi)

$$ \tag{10} \Psi_i=[m_i-z_i]_{24},\qquad R=\left|\frac1n\sum_i e^{i2\pi\Psi_i/24}\right|,\qquad SD_c=\sqrt{-2\ln R}\frac{24}{2\pi} $$

Engine: entrainment.py:99-409, using circular.py:86-173.

Step 11 — Measure re-entrainment

Tolerance is the larger of a multiple of baseline circular spread and one analysis bin. The first sustained run inside tolerance is reported only if the final day remains settled. When supported, an exponential describes approach.

(11) tolerance = MAX(k*baseline_SD, bin_minutes/60); deviation(day) = A*EXP(-day/tau_days); half_life = tau_days*LN(2)

$$ \tag{11} \epsilon=\max(kSD_{base},B/60),\qquad d(t)=Ae^{-t/\tau_r},\qquad t_{1/2}=\tau_r\ln2 $$

Engine: entrainment.py:411-703. The exponential decision rule is original to this project; no citation is present in source.

Step 12 — Measure acute masking

Activity during the pulse is compared with the same clock window on prior control days. The prediction-interval denominator treats the pulse as a new observation, not another member of the control mean.

(12) t = (pulse-control_mean)/(control_SD*SQRT(1+1/n_control)); difference = pulse-control_mean

$$ \tag{12} t=\frac{x_p-\bar x_c}{s_c\sqrt{1+1/n_c}},\qquad \Delta x=x_p-\bar x_c $$

Engine: entrainment.py:705-900. This test choice is original to this project; no citation is present in source.

Step 13 — Score sleep, food anticipation, and ultradian activity

Sustained immobility uses a run-length threshold. Food anticipation contrasts declared pre- and post-meal windows. Ultradian estimation subtracts a fitted slow component before searching the faster band and checks simple harmonics.

(13) min_samples = CEIL(immobility_minutes/sample_interval_minutes); FAA_index = premeal/(premeal+postmeal); ultradian_resolution = fast_period^2/record_span

$$ \tag{13} n_{min}=\left\lceil\frac{T_{immobile}}{\Delta t}\right\rceil,\qquad I_{FAA}=\frac{x_{pre}}{x_{pre}+x_{post}},\qquad \delta P=\frac{P_f^2}{T_{span}} $$

Engine: sleep.py:66-289, feeding.py:42-356, and ultradian.py:73-285. Pack et al. (2007) is cited for the video-immobility convention; food and ultradian decision rules have no citation in source.

Step 14 — Test temperature compensation

Rate is the reciprocal of period. With at least three temperatures, the log-rate slope gives Q10; repeated unit labels are retained for clustered uncertainty.

(14) rate = 1/period; LN(rate) = a + b*temperature; Q10 = EXP(10*b)

$$ \tag{14} r_i=P_i^{-1},\qquad \ln r_i=a+bT_i+\epsilon_i,\qquad Q_{10}=e^{10b} $$

Engine: compensation.py:50-292. The 0.8–1.2 acceptance band is a configurable convention and has no citation in source.

Step 15 — Measure population phase coherence

Each component trace supplies an analytic phase. The Kuramoto order parameter averages unit phase vectors; ensemble amplitude is decomposed into mean component amplitude and phase coherence.

(15) R(t) = ABS(SUM(EXP(i*phase[cell,t]))/N); ensemble_amplitude(t) ~= mean_cell_amplitude(t)*R(t)

$$ \tag{15} R(t)=\left|\frac1N\sum_j e^{i\theta_j(t)}\right|,\qquad A_{ens}(t)\approx\overline{A_j(t)}R(t) $$

Engine: synchrony.py:1-745. The damping decomposition follows Yoo et al. (2004).

Step 16 — Map space, lag, locking, and period dispersion

Spatial phase fits the plane that maximizes residual circular concentration. Pairwise lag and phase-locking value (PLV) use circular phase difference. Period dispersion winds starting phases forward at each component's period to predict loss of coherence.

(16) spatial_Q(b,c) = ABS(MEAN(EXP(i*(phase-b*x-c*y)))); lag = P*ARG(MEAN(EXP(i*(phase_B-phase_A))))/(2*pi); predicted_phase[i,t] = phase0[i] + 2*pi*t/tau[i]

$$ \tag{16} Q(b,c)=\left|\frac1N\sum_j e^{i(\theta_j-bx_j-cy_j)}\right|,\quad \ell=\frac{P}{2\pi}\arg\overline{e^{i(\theta_B-\theta_A)}},\quad \theta_i^{pred}(t)=\theta_{i0}+\frac{2\pi t}{\tau_i} $$

Engine: spatial.py:1-479, coupling.py:118-578, and synchrony.py:749-1195. Whole cycles or component traces, not timepoints, form the uncertainty units.

Step 17 — Plan prospective power

The simulator repeatedly generates data under the declared smallest meaningful effect and runs the real target test. The first recommended size is based on the Wilson lower confidence bound, not the noisier point estimate.

(17) power_hat = rejections/simulations; wilson_centre = (power_hat + z^2/(2*S))/(1+z^2/S)

$$ \tag{17} \widehat\pi=R/S,\qquad c_W=\frac{\widehat\pi+z^2/(2S)}{1+z^2/S} $$

Engine: power.py:1-605. Post-hoc observed power is refused following Hoenig and Heisey (2001).

Step 18 — Predict entrainment, re-entrainment, and sleep pressure

Phase-response predictions iterate a once-per-cycle phase map; a fixed point is stable only when its response-curve derivative lies between -2 and 0. The two-process model evolves homeostatic pressure with separate wake and sleep time constants and a sinusoidal circadian threshold.

(18) phase[n+1] = (phase[n] + 24*T/tau + shift(phase[n])) MOD 24; S_wake_next = 1-(1-S)*EXP(-dt/tau_rise); S_sleep_next = S*EXP(-dt/tau_fall)

$$ \tag{18} \phi_{n+1}=\left[\phi_n+\frac{24T}{\tau}+\Delta(\phi_n)\right]{24},\quad S_w(t+\Delta t)=1-[1-S(t)]e^{-\Delta t/\tau $$}},\quad S_s(t+\Delta t)=S(t)e^{-\Delta t/\tau_{fall}

Engine: prediction.py:1-658 and two_process.py:1-987. The phase-map implementation and two-process fitting choices have no complete citation in source and are labelled original to this project where they extend the named models.

Where each step can mislead

• Import cannot infer the scientific meaning of an unfamiliar numeric column; inspect retained labels and units.
• A folded profile hides drift and transients; use instantaneous phase when stationarity is doubtful.
• A period at a search boundary, or at half/twice another estimate, can be a range artefact or harmonic.
• Low analytic amplitude makes Hilbert phase unstable even when an array can be drawn.
• Calendar analyses are withheld for elapsed-hours files because their dates were synthesized.
• Phase-response zones outside observed phase coverage are extrapolations.
• Masking is immediate activity change, not clock resetting.
• Immobility is not electroencephalographic sleep.
• Correlation, phase locking, and spatial gradients do not prove causal coupling.
• Simulation and prediction inherit every assumption of the fitted input; they are not measurements.
• null means withheld or unidentified, never zero and never evidence of no effect.

Worked examples

Analytic phase-response curve

Eight points sampled from -cos(2πCT/24) every three circadian hours recover coefficients [0, -1, 0], R² 1, and extrema ±1 hour under the one-harmonic model. The case is committed as validation/cases/phase-response-sinusoid.json.

Temperature identity

Periods of 24 hours at 20, 25, and 30 °C have zero log-rate slope and Q10 1. A confidence interval is correctly withheld because one unit per temperature contains no within-temperature spread.

Expected model refusal

A two-process prediction supplied with an empty fit returns input_error. The absence of parameters is not silently replaced by defaults.

Interpretation guide

Result state	Meaning
Value with diagnostics	The declared calculation was identified on the supplied input.
Warning	A value exists but a named assumption, coverage, or interpretation boundary needs review.
null field	The quantity was not identified or was deliberately withheld.
Input refusal	The requested number has no defined meaning for those inputs.
Independent PASS	A predeclared committed reference agreed within its predeclared tolerance.
Not externally validated	Implementation tests pass, but no independent distributable result supports an external agreement claim.

Symbols

Symbol	Definition
(x_i)	observed activity or measurement at sample (i)
(t_i)	elapsed time of sample (i)
(B)	analysis-bin width in minutes
(P), τ	reference or free-running period in hours
(M)	cosinor midline estimating statistic of rhythm
(A)	amplitude; in Equation 5 only, finite adjacent-pair count
φ, θ	circular phase in hours or radians as stated
Ψ	signed marker-to-zeitgeber phase angle
(R)	circular resultant length or rejection count where explicitly subscripted
(S)	template score, homeostatic pressure, or simulation count as locally labelled
(Q_{10})	rate change for a 10 °C increase
PLV	phase-locking value, the resultant length of pairwise phase differences

Technical reference

Files and configuration

The registry and parameter contract are in actions.py; shared defaults and normalization are in analysis.py; import is in io.py. Each specialised analysis lives in its named scientific module. Numerical modules do not import FastAPI or browser code. The browser catalogue is a reviewed allow-list in ui_catalogue.py.

Normalized settings clamp declared numeric ranges and fill defaults before calculation. Every browser run uses the active recording and the same current normalized configuration as the core dashboard. Direct point/model actions skip recording import but still record their inputs and version.

Naming and outputs

Action results are JSON-clean dictionaries. Recording results include source name and fingerprint; workbench downloads wrap action, title, inputs, configuration, result, version, and generation time. Publication bundles preserve exact plotted data and source SVGs. Projects preserve replay state but not source bytes.

Validation evidence

The manifest-driven runner dispatches the real registry action. Current independent hand-calculated passes cover a timestamped import, an analytic phase-response curve, and the constant-period Q10 identity. All other user-facing methods are explicitly not externally validated in the distributable evidence matrix, even when their unit and known-answer tests pass.

Verification commands

python scripts/run_validation.py --all
python scripts/check_docs.py
python -m pytest tests/test_validation_cases.py tests/test_analysis_explanation.py
mkdocs build --strict

References

• Barrodale I, Erickson RE. Algorithms for least-squares linear prediction and maximum entropy spectral analysis. 1980. Citation as recorded in source.
• Burg JP. Maximum entropy spectral analysis. 1967. Citation as recorded in source.
• Costa MJ et al. Spectrum resampling for period estimation. Biostatistics. 2013. Citation as recorded in source.
• Dörrscheidt GJ, Beck L. F-periodogram method. 1975. Citation as recorded in source.
• Edwards KD et al. Quantitative analysis of regulatory flexibility under changing environmental conditions. Molecular Systems Biology. 2010;6:424.
• Enright JT. The search for rhythmicity in biological time-series. 1965. Citation as recorded in source.
• Fisher NI. Statistical Analysis of Circular Data. 1993.
• Hoenig JM, Heisey DM. The abuse of power: the pervasive fallacy of power calculations for data analysis. The American Statistician. 2001.
• Hughes ME, Hogenesch JB, Kornacker K. JTK_CYCLE: an efficient nonparametric algorithm for detecting rhythmic components. Journal of Biological Rhythms. 2010;25:372.
• Hutchison AL et al. Improved statistical methods enable greater sensitivity in rhythm detection for genome-wide data. PLoS Computational Biology. 2015;11:e1004094.
• Lomb NR. Least-squares frequency analysis of unequally spaced data. Astrophysics and Space Science. 1976.
• Pack AI et al. Novel method for high-throughput phenotyping of sleep in mice. Physiological Genomics. 2007.
• Plautz JD et al. Quantitative analysis of Drosophila period gene transcription in living animals. Journal of Biological Rhythms. 1997;12:204.
• Scargle JD. Studies in astronomical time series II: statistical aspects of spectral analysis of unevenly spaced data. The Astrophysical Journal. 1982.
• Sokolove PG, Bushell WN. The chi-square periodogram: its utility for analysis of circadian rhythms. Journal of Theoretical Biology. 1978.
• Ulrych TJ, Bishop TN. Maximum entropy spectral analysis and autoregressive decomposition. 1975. Citation as recorded in source.
• Yoo SH et al. PERIOD2::LUCIFERASE real-time reporting of circadian dynamics reveals persistent circadian oscillations in mouse peripheral tissues. Proceedings of the National Academy of Sciences. 2004;101:5339.
• Zar JH. Biostatistical Analysis. 4th ed. 1999.
• Zielinski T et al. Strengths and limitations of period estimation methods for circadian data. PLoS ONE. 2014;9:e96462.

Where the engine names an algorithm without bibliographic detail, the relevant step above says that the citation is absent rather than inventing one.