Improved estimation and interpretation of correlations in neural circuits, as a DataJoint pipeline in Python.
This repository reproduces the analysis of:
Dimitri Yatsenko, Krešimir Josić, Alexander S. Ecker, Emmanouil Froudarakis, R. James Cotton, and Andreas S. Tolias (2015). Improved estimation and interpretation of correlations in neural circuits. PLOS Computational Biology 11(3): e1004083. doi:10.1371/journal.pcbi.1004083
It is a translation of the study's original code and data,
atlab/cov-est (MIT License, © 2014 Tolias Lab), from DataJoint
for MATLAB to DataJoint 2 for Python. The schema follows the original covest schema table for
table, and the estimators are ported from its +cove package to numpy and scipy.
- The study is by the authors above. The article is published under the Creative Commons Attribution 4.0 License.
- The data are the two-photon recordings published with the original code in atlab/cov-est/data. They are not copied into this repository: the pipeline downloads them from a pinned commit and verifies their SHA-256.
- The original code, atlab/cov-est, is © 2014 Tolias Lab under the MIT License. This repository is a derivative work under the same license; see LICENSE.
- To cite this work, cite the article. CITATION.cff lists the article and the original code.
The pipeline is the covest schema, from Session on. A separate ingest schema handles the
published data files, one import job per file: each job reads its file and fills the pipeline's
first tables, which have no reference back to it (fan-out ingestion), so the diagram draws the
two apart.
| Schema | Table | Tier | One row is | Original |
|---|---|---|---|---|
ingest |
SourceFile |
Lookup | a data file published with the original code (atlab/cov-est at a pinned commit), with its SHA-256 and the number of its first trace set | (new) |
ingest |
SourceImport |
Ingest | the import of one file's sessions, sites, and trial trace sets into the pipeline | (new) |
covest |
Session |
Entry | an imaging session: one mouse on the setup | (new) |
covest |
Site |
Entry | an imaging site: one 3D scan of a session | (new) |
covest |
TrialTraceSet |
Entry | the binned activity of one imaging site in every trial of one stimulus set | covest.Traces |
covest |
CellSelection |
Compute | the cells of a correlation set that stay active throughout the recording | covest.ActiveCells |
covest |
OriTuning |
Compute | the von Mises orientation tuning of every cell of a tuning set, with its significance by 10,000 shuffles | covest.OriTuning |
covest |
Method |
Lookup | a covariance estimator and the grids of its hyperparameters | covest.Method |
covest |
Fold |
Lookup | a fold of 10-fold cross-validation, or all trials (nfolds = 1) |
covest.Fold |
covest |
CovMatrix |
Compute | a correlation matrix estimated on a fold's training trials, its hyperparameters, and its validation loss on the fold's test trials | covest.CovMatrix |
Each table is named for what one row is, so two tables of the original are renamed:
covest.Traces is TrialTraceSet (the name of the lab table it was copied from,
pop.TrialTraceSet), and covest.ActiveCells is CellSelection.
Keys. The original table was keyed by its identifiers in the lab's database: subject,
setup, session and scan start times, preprocessing method, and stimulus set. Here a session and
a site take the numbers the original code and figure scripts label sites by, the last four
digits of their start times: session is mod(session_start_time, 10000), and a site is keyed
by its session and site, mod(aod_scan_start_time, 10000), so its key says which animal it
comes from and matches the original code. A trial trace set has a simple number (trace_set)
and refers to its site; trace sets are numbered in the order of the published files, each file
from the first number recorded in SourceFile (1, 17, 42, 57). Every number is fixed by its
record or its file, so each import job reads only its own file and the jobs can run in any
order; a session with records in two files is inserted by whichever job comes first. The
original identifiers are kept as unique attributes (session_start_time,
aod_scan_start_time), and the analysis tables are keyed by trace_set. The example site of
Fig 3 (292 cells) is session 8828, site 2031, with tuning set 27 and correlation set 28. Random
draws (the hyperparameter search, the tuning shuffles) are seeded from the scan identifier.
Subjects. The published records carry a subject_id, but it is 16 in all 60 of them,
although the study imaged 14 mice, so it does not identify the animals and is not kept. The 27
analyzed sites come from 14 sessions, which agrees with one session per mouse.
Each site is a population of 106–328 neurons in layers 2/3 of mouse primary visual cortex,
recorded with 3D random-access two-photon microscopy. Activity is inferred firing rate, binned in
150 ms bins aligned on stimulus onset. Each site has two stimulus sets (high_repeats):
- Correlation sets (31 sites): drifting gratings in 2 or 5 directions, about 100–300 trials of each, with 1 s of grating and 1 s of blank per trial. Noise correlations are estimated from these.
- Tuning sets (29 sites): 16 directions, 12–34 trials of each, for orientation tuning.
trace_segments is a bins × directions × trials × cells array, padded with NaN where a direction
has fewer trials than others.
Every array (traces, cell positions, correlation matrices, means, variances, search traces) is a
<npy@> attribute: a standard .npy file in the default object store, at a path that mirrors the
schema, fetched lazily (np.asarray loads it). datajoint.json points the default store at
data/store in the repository; any file, S3, GCS, or Azure store works (DJ_STORES or the
stores section of datajoint.json).
method |
Estimator | Hyperparameters (grid) | Python implementation |
|---|---|---|---|
| 0 | sample correlation, all bins | none | numpy |
| 1 | sample correlation, evoked bins only | none | numpy |
| 2, 3 | sample correlation, evoked bins of one direction (two-direction sites) | none | numpy |
| 10 | shrinkage toward a diagonal target (Cdiag) | variance shrinkage toward the median, e−6…e−1; correlation shrinkage toward zero, e−4…1 | numpy |
| 30 | factor analysis: low rank plus diagonal (Cfactor) | variance shrinkage, e−6…e−1; factors, 0…70 | scipy: maximum likelihood by L-BFGS-B, loadings in closed form |
| 80 | sparse precision, graphical lasso (Csparse) | α, e−5…e−1 | gglasso ADMM |
| 90 | sparse plus latent (Csparse+latent) | α, e−7…e−1; β, e−5…e−1 | gglasso ADMM |
| 91 | sparse plus latent with the number of latent units capped (Fig 7) | α, e−7…e−1; at most 0…100 latent units; β = 0.001 | numpy port of the original ADMM |
| 100 | sparse plus latent, narrower grid (Figs 5, 6, S1) | α, e−7…e−1.5; β, e−5…e−1.5 | gglasso ADMM |
All grids step by 0.05 in natural log. As in the original:
- Factor analysis shrinks the variances toward their median before the fit. The fit is the same maximum-likelihood factor model as the original's EM; it converges where the EM's 200-cycle cap stops short.
- Sparse and sparse plus latent fit the covariance scaled to unit mean variance and penalize the whole of S, diagonal included: minimize tr(Σ(S − L)) − log det(S − L) + α Σ|Sij| + β tr(L) with L positive semidefinite (Eq. 18 with the penalties in the original's units, without the 1/p). gglasso's solver reaches the same optimum as regain and cvxpy; at the paper's S1 optimum for site 2031 (α ≈ 0.017, β ≈ 0.15) it gives 93.45% sparsity and 73 latent units (paper: 92.8%, 72).
- Method 91 caps the rank of L instead of penalizing it; no library solves that problem, so it runs the original's ADMM, ported to numpy.
OriTuning fits each cell's mean response per direction (the sum of the evoked bins of a trial)
with a two-peaked von Mises curve, f(φ) = a + b exp(κ(cos(φ − θ) − 1)) + c exp(κ(cos(φ − θ + π) − 1))
with b ≥ c ≥ 0, by the original's search (seven starting directions, then parabolic steps in θ
and a bisection in κ). R² is computed against single trials, and its p-value from 10,000 shuffles
of the direction-trial labels, the same for all cells; a cell is tuned if p < 0.05. At site 2031,
219 of 292 cells are tuned.
Each method is estimated on the training trials of 10 folds and evaluated on the held-out trials
with the normal (Gaussian) loss, in nats per cell per bin. Hyperparameters are chosen by an inner
10-fold cross-validation on the training trials: a random search seeds a pattern search on the
logarithmic grids in Method. The fold nfolds = 1 estimates from all trials, with hyperparameters
chosen the same way.
- Fold membership. Trial slots are assigned to folds by a fixed numpy permutation instead of
MATLAB's
rng(0); randperm. The structure is the same: each trial slot belongs to one fold for every direction, and all methods see the same folds. - Variance regularization.
find_best_deltareturns the value with the lowest loss on its grid.findBestDelta.mreturns the value one step past it. - Random search. Seeded from the key of each
CovMatrixrow, so results are reproducible. - Columns.
Method.conditionisdirection: 0 for all directions, 1 or 2 for the first or second direction of a two-direction set only (conditionis a reserved word in SQL).Method.hyperparam_spaceis JSON.CovMatrix.test_matrix(a copy of the test trials) is dropped, and so isQuadLoss, which is not in the paper.
The paper's Figures 3 to 7 and panel A of Figure S1, drawn from this pipeline's results, are pages of the companion dashboard, datajoint/neurocorr-dash (MIT License). Each page reads the pipeline as it is, so it includes any site added later.
# install
uv sync --group dev # or: pip install -e .
# start a local PostgreSQL database
cd local && docker compose up -d && cd ..
mkdir -p .secrets && echo datajoint > .secrets/database.password
mkdir -p data/store # the default object store in datajoint.jsonfrom neurocorr.pipeline import covest, ingest
ingest.SourceImport.populate() # one job per data file (218 MB in all): 15 sessions, 31 sites, 60 trial trace sets in covest
covest.CellSelection.populate() # 31 correlation sets; 27 keep more than 100 cells
covest.OriTuning.populate() # 29 tuning sets
covest.CovMatrix.populate() # every method, site, and foldThe data files are cached in ~/.cache/neurocorr, or in the directory named by
NEUROCORR_DATA_DIR. Schema names take the prefix in the DATABASE_PREFIX environment variable.
The bare <npy@> form needs datajoint-python after release 2.3.4, so the dependency points at its
master branch until the next release. For an S3-compatible store, install the s3 extra
(pip install -e ".[s3]"), which installs DataJoint's own s3 extra (s3fs).
datajoint.json holds the non-secret settings; credentials come from .secrets/ or the DJ_USER
and DJ_PASS environment variables. To use MySQL instead of PostgreSQL, set DJ_BACKEND=mysql,
DJ_HOST, and DJ_PORT.
uv sync --group dev
pytest tests/unit # fast, no database
pytest # all tests; starts PostgreSQL with testcontainers (requires Docker)
ruff check src tests && ruff format src testsneurocorr/
├── src/neurocorr/
│ ├── pipeline/covest.py # the pipeline schema, from Session on
│ ├── pipeline/ingest.py # the import of the published data files
│ ├── cove/ # estimators and cross-validation (port of +cove)
│ └── ingest/ # fetching, verifying, and reading the original data files
├── tests/unit, tests/integration
├── local/docker-compose.yml # PostgreSQL for development
├── datajoint.json
├── CITATION.cff
└── LICENSE # MIT, as atlab/cov-est
MIT, the license of the original code; see LICENSE. The article is CC-BY 4.0.