Kauvar, Richman, Liu et al. 2025: mouse Neuropixels claims (DANDI:001326)¶
Paper: Kauvar, Richman, Liu et al., "Conserved brain-wide emergence of emotional response from sensation in humans and mice", Science 388, eadt3971 (2025), doi:10.1126/science.adt3971.
Data: the mouse Neuropixels recordings are not in dandiset 001347; the paper deposited them in DANDI:001326 (ref. 104), so these claims are data: elsewhere. DANDI:001326 has one published version, 0.250528.1957, which is pinned; the paper cites no version. Files are streamed with remfile and h5py (dandi_io.open_h5); nothing is downloaded. The Neuropixels-rig behaviour files of DANDI:001347 at 0.250528.0702 are read only to link sessions.
Text read by the claims file: the version of record with its supplement (Materials and Methods, figure legends, Table S3), as extracted in paper/10.1126-science.adt3971.txt. Methods details not in the claims file come from that text and from the mouse analysis code, Zenodo record 14876694 (analysis.py, md5 6cc7aa5059737c3bd753aad1f116019d), which is read, not imported; mouse_npx.py reimplements what is needed and its docstring says where every parameter comes from.
| claim | tested here |
|---|---|
| fig6d-mouse-persistent-activity-reduced | tested |
| fig6f-mouse-fast-activity-unchanged | tested |
| fig5e-mouse-emotion-dimension-late-reduced | tested |
| fig6h-mouse-timescale-shorter | tested |
| fig6j-mouse-coupling-reduced | tested |
| fig5d-mouse-emotion-dimension-accumulates | tested |
| figs34-mouse-tone-persistence-reduced | tested |
| fig4hi-mouse-rise-unchanged-decay-shortened | tested |
| fig4g-mouse-slow-clusters-altered | tested |
| fig5g-mouse-first-order-model | tested (neural fit on DANDI:001326, dose fits on DANDI:001347) |
import os
os.environ.setdefault('OMP_NUM_THREADS', '1')
os.environ.setdefault('OPENBLAS_NUM_THREADS', '1')
import time, warnings
import numpy as np
import pandas as pd
import scipy.stats
import yaml
import matplotlib.pyplot as plt
import dandi_io
import mouse_npx as m
pd.set_option('display.width', 200)
pd.set_option('display.max_columns', 30)
pd.set_option('display.precision', 4)
with open('../../claims/001347/10.1126-science.adt3971.yaml') as fh:
CLAIMS = {c['id']: c for c in yaml.safe_load(fh)['claims']}
def expected(cid):
return [(e['quantity'], e['value'], e['tolerance']) for e in CLAIMS[cid].get('expected', [])]
def within(x, v, tol):
return abs(x - v) <= tol
Data check¶
Every NWB file of DANDI:001326 is streamed once. mouse_npx.process_session bins each unit's spikes at 10 ms, applies the 50-ms causal moving average, z-scores each unit over the recording, finds the puff-modulated cells, and keeps only the reduced arrays the claims need. This cell takes most of the run time.
t0 = time.time()
paths = m.npx_paths()
Rs = []
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
for p in paths:
Rs.append(m.process_session(p))
print(f'{len(paths)} files read and reduced in {(time.time() - t0) / 60:.1f} min')
rows = []
for R in Rs:
st = R['puff']
rows.append(dict(session=f"sub-{R['subject']}_{R['session']}", mouse=R['subject'], probes=R['n_probes'],
units=R['n_units'], modulated=R['n_sel'], puff_series=st.shape[0],
puff_isi_s=np.median(np.diff(st, axis=1)),
series_before_injection=int(np.sum(st[:, 0] < R['ket_time'])),
tone_series=0 if R['tone'] is None else R['tone'].shape[0],
light_series=0 if R['light'] is None else R['light'].shape[0],
units_sorted_spikes=int(R['spikes_sorted'].sum())))
SESS = pd.DataFrame(rows)
print(SESS.to_string())
print('mice:', SESS.mouse.nunique(), 'sessions:', len(SESS))
13 files read and reduced in 2.7 min
session mouse probes units modulated puff_series puff_isi_s series_before_injection tone_series light_series units_sorted_spikes
0 sub-0_20230901 0 4 942 754 60 3.0031 20 29 30 0
1 sub-1_20230218 1 4 583 404 50 3.0031 20 0 0 0
2 sub-1_20230226 1 4 1020 621 61 3.0031 20 30 31 0
3 sub-2_20230902 2 4 708 451 60 3.0031 20 29 30 0
4 sub-3_20230218 3 4 431 310 50 3.0031 20 0 0 0
5 sub-3_20230225 3 4 917 579 62 3.0031 20 31 31 0
6 sub-4_20230907 4 4 885 649 60 3.0031 20 30 30 0
7 sub-5_20230706 5 2 154 65 60 3.0031 20 30 30 0
8 sub-5_20230822 5 4 664 378 60 3.0031 20 29 30 0
9 sub-6_20230823 6 4 591 415 60 3.0031 20 29 30 0
10 sub-7_20220831 7 4 928 638 52 3.0137 20 0 0 0
11 sub-8_20220228 8 4 851 670 40 3.0100 20 0 0 0
12 sub-9_20230908 9 4 880 619 60 3.0031 20 30 30 0
mice: 10 sessions: 13
# Unit metadata against the Methods (Bombcell 'somatic' units, SNR > 2, > 500 spikes; ISI violations < 0.1)
q = np.concatenate([R['quality'] for R in Rs]); snr = np.concatenate([R['snr'] for R in Rs])
ns = np.concatenate([R['n_spikes'] for R in Rs])
print('units:', len(q), 'quality labels:', dict(zip(*np.unique(q, return_counts=True))))
print(f'SNR range {snr.min():.2f} to {snr.max():.1f}; n_spikes range {ns.min():.0f} to {ns.max():.0f}')
print('units with SNR <= 2:', int((snr <= 2).sum()), '; units with <= 500 spikes:', int((ns <= 500).sum()))
# stimulus durations and the injection time
dur = {}
for R in Rs:
f = dandi_io.open_h5(R['path'], dandi_io.NPX_DANDISET, dandi_io.NPX_VERSION)
st = f['intervals/Stimulus presentations/stimulus'][:].astype(str)
d = f['intervals/Stimulus presentations/stop_time'][:] - f['intervals/Stimulus presentations/start_time'][:]
for s in np.unique(st):
dur.setdefault(s, []).append(d[st == s])
f.close()
for s, v in dur.items():
v = np.concatenate(v)
print(f'{s}: {len(v)} presentations, duration {v.min():.3f} to {v.max():.3f} s')
gap = [R['ket_time'] - R['puff'][19, -1] for R in Rs]
nxt = [R['puff'][20, 0] - R['ket_time'] for R in Rs]
print(f"'Ketamine infusion start' falls {min(gap):.0f} to {max(gap):.0f} s after the last puff of series 20 "
f"and {min(nxt):.0f} to {max(nxt):.0f} s before series 21")
print('modulated cells per session:', SESS.modulated.tolist(), '; fraction of units', round(SESS.modulated.sum() / SESS.units.sum(), 3))
units: 9554 quality labels: {np.str_('good'): np.int64(9554)}
SNR range 2.00 to 825.8; n_spikes range 517 to 807685
units with SNR <= 2: 0 ; units with <= 500 spikes: 0
light: 2176 presentations, duration 0.250 to 0.250 s puff: 5880 presentations, duration 0.250 to 0.250 s tone: 2136 presentations, duration 0.250 to 0.250 s 'Ketamine infusion start' falls 124 to 450 s after the last puff of series 20 and 5 to 5 s before series 21 modulated cells per session: [754, 404, 621, 451, 310, 579, 649, 65, 378, 415, 638, 670, 619] ; fraction of units 0.686
Checks against the description:
- 13 sessions from 10 mice (subjects 0 to 9), as the paper states (Table S3, Fig. 4). Ten sessions have 4 probes; sub-5 ses-20230706 has 2 (the Methods say 3 to 4).
- Every unit is labelled 'good', with SNR above 2 and more than 500 spikes, the Methods' criteria; the ISI criterion cannot be checked because the files hold no ISI metric. The
n_spikescolumn equals the number of stored spike times. - Puffs, tones and lights last 250 ms, puffs come 3.0 s apart in series of 8 (the Methods give 250 ms and 3 s). Tones are in 9 sessions (the paper's tone analysis has 7 sessions from 7 mice).
- The 'Ketamine infusion start' timestamp falls between series 20 and 21 in every session (the Methods: infusion after the 20th trial).
- The spike times of every unit are stored unsorted (column
units_sorted_spikesis 0 for all sessions), which breaks the NWB expectation that a unit's spike times are ascending; they are sorted here before binning.
Matching sessions to the paper¶
The files carry subject numbers 0 to 9 and recording dates, not the paper's mouse or experiment identifiers. The Neuropixels-rig behaviour files in DANDI:001347 (named by experiment id, sub-m<mouse>-<experiment>) hold the trial onsets (events-on, 100 frames/s). The interval between series 20 and 21 contains the injection and differs from session to session; it identifies each behaviour file's Neuropixels session, and the whole interval sequence is then compared.
beh = {}
for p in sorted(dandi_io.list_assets()):
if p.startswith('sub-m') and not p.startswith('sub-mouse'):
f = dandi_io.open_h5(p)
beh[p.split('/')[0][4:]] = f['acquisition/events-on/data'][:] / f['acquisition/fps/data'][0]
f.close()
rows = []
for name, on in beh.items():
gb = np.diff(on)
best, err = None, np.inf
for R in Rs:
gn = np.diff(R['puff'][:, 0])
n = min(len(gb), len(gn))
e = np.max(np.abs(gb[:n] - gn[:n]))
if e < err:
best, err, nn = R, e, n
rows.append(dict(behaviour=name, trials=len(on), npx_session=f"sub-{best['subject']}_{best['session']}",
npx_series=best['puff'].shape[0], intervals_compared=nn, max_abs_diff_s=err,
gap_20_21_s=np.diff(on)[19]))
LINK = pd.DataFrame(rows)
print(LINK.to_string())
unmatched = sorted(set(SESS.session) - set(LINK.npx_session))
print('Neuropixels sessions without a behaviour file:', unmatched)
tone_sessions = SESS[SESS.tone_series > 0]
print('tone sessions:', len(tone_sessions), 'from', tone_sessions.mouse.nunique(), 'mice')
for mouse in ['1', '3', '5']:
rest = tone_sessions[tone_sessions.mouse != mouse]
rest_all = SESS[SESS.mouse != mouse]
print(f'excluding mouse {mouse}: population {len(rest_all)} sessions / {rest_all.mouse.nunique()} mice; '
f'tone {len(rest)} sessions / {rest.mouse.nunique()} mice')
behaviour trials npx_session npx_series intervals_compared max_abs_diff_s gap_20_21_s 0 m11208-959 48 sub-1_20230218 50 47 0.0129 154.82 1 m11208-962 59 sub-1_20230226 61 58 0.0151 434.09 2 m11210-960 48 sub-3_20230218 50 47 0.0131 180.81 3 m11219-1108 58 sub-6_20230823 60 57 0.0129 165.93 4 m11221-1119 58 sub-0_20230901 60 57 0.0120 150.00 5 m11223-1125 58 sub-2_20230902 60 57 0.0120 154.01 6 m11225-1131 58 sub-4_20230907 60 57 0.0123 149.91 7 m11227-1132 58 sub-9_20230908 60 57 0.0122 298.04 8 m304-876 50 sub-7_20220831 52 49 0.0186 475.72 9 m94-549 79 sub-8_20220228 40 39 0.0185 216.39 Neuropixels sessions without a behaviour file: ['sub-3_20230225', 'sub-5_20230706', 'sub-5_20230822'] tone sessions: 9 from 8 mice excluding mouse 1: population 11 sessions / 9 mice; tone 8 sessions / 7 mice excluding mouse 3: population 11 sessions / 9 mice; tone 8 sessions / 7 mice excluding mouse 5: population 11 sessions / 9 mice; tone 7 sessions / 7 mice
The ten behaviour files match ten Neuropixels sessions, each with all its inter-series intervals equal within 0.1 s, so the clocks of the two files agree up to an offset. Experiment ids m11208-959 and m11208-962 are subject 1, m11210-960 is subject 3's first session; sub-3 ses-20230225 and both sub-5 sessions have no behaviour file in 001347. The behaviour files name 9 mice (10 sessions), the Neuropixels files 10 mice (13 sessions).
The Methods exclude one animal "that exhibited notable stress response" from the cross-puff and population analyses, giving 11 sessions from 9 mice; neither the paper nor the code names it. Only mice 1, 3 and 5 have two sessions. Excluding mouse 5 is the only choice that also gives the tone analysis its 7 sessions from 7 mice (Table S3, fig. S34); excluding mouse 1 or 3 leaves 8 tone sessions. Mouse 5 is therefore taken as the excluded animal. This is an assumption: neither the paper, the code nor the dataset names the animal (mouse 5 also has the 2-probe session). Every claim that uses the 11 sessions (fig6d, fig6f, fig5e, fig5d, fig6h, fig6j, figs34) also prints its main reading on all 13 sessions (9 tone sessions for figs34); fig5g uses the 11 only. Fig. 4 analyses (clusters, rise and decay) use all 13 sessions, as Table S3 gives for Fig. 4.
# Assumption (see the cell above): mouse 5 is the animal the Methods exclude. fig6d, fig6f, fig5e, fig5d,
# fig6h, fig6j and figs34 also print their main reading on all 13 sessions (9 tone sessions for figs34).
POP = [R for R in Rs if R['subject'] != '5']
TONE = [R for R in POP if R['tone'] is not None]
print('population analyses:', len(POP), 'sessions from', len({R['subject'] for R in POP}), 'mice')
print('tone analyses:', len(TONE), 'sessions from', len({R['subject'] for R in TONE}), 'mice')
# Clock check: mean preinfusion z-rate of all units around the puff (window -1 to +2 s)
t = (np.arange(m.WIN) - 100) * m.BIN
rows = []
for R in Rs:
tr = R['pop_trace_all']
base = tr[:100].mean(); sd = tr[:100].std()
d = np.diff(tr)
rows.append(dict(session=f"sub-{R['subject']}_{R['session']}",
resp_0_200ms=tr[100:120].mean() - base,
resp_shifted_minus_500ms=tr[50:70].mean() - base,
largest_rise_at_s=t[1:][np.argmax(d)], peak_at_s=t[np.argmax(tr)]))
print(pd.DataFrame(rows).to_string())
population analyses: 11 sessions from 9 mice
tone analyses: 7 sessions from 7 mice
session resp_0_200ms resp_shifted_minus_500ms largest_rise_at_s peak_at_s
0 sub-0_20230901 0.3088 -0.0018 0.02 0.05
1 sub-1_20230218 0.2394 -0.0052 0.01 0.05
2 sub-1_20230226 0.1191 -0.0019 0.01 0.05
3 sub-2_20230902 0.1082 -0.0003 0.01 0.05
4 sub-3_20230218 0.0961 -0.0031 0.02 0.05
5 sub-3_20230225 0.0974 0.0011 0.02 0.06
6 sub-4_20230907 0.1844 0.0033 0.02 0.05
7 sub-5_20230706 0.0828 0.0011 0.02 0.05
8 sub-5_20230822 0.1051 0.0014 0.04 0.32
9 sub-6_20230823 0.2657 0.0006 0.02 0.05
10 sub-7_20220831 0.1401 -0.0056 0.01 0.05
11 sub-8_20220228 0.3416 0.0016 0.02 0.08
12 sub-9_20230908 0.1757 0.0005 0.01 0.05
In every session the mean z-scored rate of all units rises steeply within the first tens of milliseconds after the stored puff time (largest one-bin rise and peak shortly after 0 s), and the same 200-ms window placed 500 ms earlier shows no response. The spike times and the stimulus table share a clock in every session.
fig6d-mouse-persistent-activity-reduced¶
In mice, activity along the persistent dimension (1 to 2 s after the puff) falls during ketamine (0.86 to 0.02) and partly recovers after (0.52).
Persistent dimension built from 0.5 to 1 s after each puff against the 1 s before it, over the 20 preinfusion series; the trial-averaged projection is z-scored over its 3-s window and the 1-s baseline subtracted (code); quantified 1 to 2 s. Readings run: the coding-dimension variance (Methods, per neuron 'M'; code, pooled 'C'), the epochs (Methods trials 22 to 32 and 34 to 40; code and Table S3 23 to 32 and 35 to 40), and the first puff of each series (code) or all 8 puffs. The main reading is M, first puff, code/Table S3 epochs: the Methods and Table S3 disagree on the trial counts and the code agrees with Table S3 (an assumption).
cid = 'fig6d-mouse-persistent-activity-reduced'
def run_readings_3way(table_fn, readings, pop):
rows = []
for name, args in readings.items():
d = table_fn(pop, *args)
s = m.three_way_tests(d)
row = dict(reading=name)
for k in ('pre', 'ket', 'post'):
mu, se = m.mean_sem(d[k]); row[k] = f'{mu:.3f} ± {se:.3f}'
row.update({k: v for k, v in s.items()})
rows.append(row)
return pd.DataFrame(rows)
R6D = {}
for var in 'MC':
for first in (True, False):
for ename, (ket, post) in m.EPOCH_READINGS.items():
R6D[f"var {var}, {'first puff' if first else 'all puffs'}, {ename}"] = (var, first, ket, post)
tab6d = run_readings_3way(m.persist_table, R6D, POP)
print(tab6d.to_string())
MAIN6D = ('M', True, 'ket_C', 'post_C')
d6d = m.persist_table(POP, *MAIN6D)
print('\nmain reading, per session:'); print(d6d.to_string())
s6d = m.three_way_tests(d6d)
print({k: round(v, 5) for k, v in s6d.items()})
for (qname, v, tol), col in zip(expected(cid), ('pre', 'ket', 'post')):
mu, se = m.mean_sem(d6d[col])
print(f'{qname}: {mu:.3f} ± {se:.3f} (expected {v} ± {tol}; within: {within(mu, v, tol)})')
dm = m.by_mouse(d6d, ['pre', 'ket', 'post'])
print('\naveraged within mouse (9 mice):', {k: round(v, 4) for k, v in m.three_way_tests(dm).items()})
d13 = m.persist_table(Rs, *MAIN6D)
print('all 13 sessions:', {k: (round(m.mean_sem(d13[k])[0], 3)) for k in ('pre', 'ket', 'post')},
{k: round(v, 4) for k, v in m.three_way_tests(d13).items()})
reading pre ket post t_pre_ket p_pre_ket q_pre_ket t_ket_post p_ket_post q_ket_post
0 var M, first puff, Methods (22-32, 34-40) 0.867 ± 0.059 0.134 ± 0.121 0.431 ± 0.217 5.1777 0.0004 0.0008 -1.4839 0.1687 0.1687
1 var M, first puff, code / Table S3 (23-32, 35-40) 0.867 ± 0.059 -0.002 ± 0.133 0.515 ± 0.219 5.4830 0.0003 0.0005 -2.5396 0.0294 0.0294
2 var M, all puffs, Methods (22-32, 34-40) 0.455 ± 0.030 0.152 ± 0.060 0.424 ± 0.097 4.7183 0.0008 0.0016 -2.8514 0.0172 0.0172
3 var M, all puffs, code / Table S3 (23-32, 35-40) 0.455 ± 0.030 0.073 ± 0.068 0.474 ± 0.100 5.4100 0.0003 0.0006 -3.8678 0.0031 0.0031
4 var C, first puff, Methods (22-32, 34-40) 0.855 ± 0.050 0.172 ± 0.118 0.423 ± 0.213 5.0958 0.0005 0.0009 -1.3387 0.2103 0.2103
5 var C, first puff, code / Table S3 (23-32, 35-40) 0.855 ± 0.050 0.019 ± 0.127 0.513 ± 0.215 5.6441 0.0002 0.0004 -2.6148 0.0258 0.0258
6 var C, all puffs, Methods (22-32, 34-40) 0.441 ± 0.022 0.143 ± 0.061 0.346 ± 0.107 5.0523 0.0005 0.0010 -1.8372 0.0960 0.0960
7 var C, all puffs, code / Table S3 (23-32, 35-40) 0.441 ± 0.022 0.052 ± 0.070 0.394 ± 0.112 5.8753 0.0002 0.0003 -2.7357 0.0210 0.0210
main reading, per session:
session mouse pre ket post
0 sub-0_20230901 0 1.2058 0.6190 1.3187
1 sub-1_20230218 1 0.9087 -0.5183 0.6650
2 sub-1_20230226 1 0.8090 -0.0686 0.8411
3 sub-2_20230902 2 0.6080 0.2437 0.5366
4 sub-3_20230218 3 1.0712 -0.1232 1.6013
5 sub-3_20230225 3 0.7964 0.2237 -0.1297
6 sub-4_20230907 4 0.7771 0.5488 1.0859
7 sub-6_20230823 6 0.7335 -0.2761 0.5152
8 sub-7_20220831 7 0.8698 0.4405 0.0371
9 sub-8_20220228 8 1.1631 -0.8461 -1.1777
10 sub-9_20230908 9 0.5992 -0.2604 0.3733
{'t_pre_ket': 5.48299, 'p_pre_ket': 0.00027, 'q_pre_ket': np.float64(0.00054), 't_ket_post': -2.53957, 'p_ket_post': 0.02938, 'q_ket_post': np.float64(0.02938)}
preinfusion: 0.867 ± 0.059 (expected 0.86 ± 0.05; within: True)
ketamine: -0.002 ± 0.133 (expected 0.02 ± 0.13; within: True)
postinfusion: 0.515 ± 0.219 (expected 0.52 ± 0.21; within: True)
averaged within mouse (9 mice): {'t_pre_ket': 4.648, 'p_pre_ket': 0.0016, 'q_pre_ket': np.float64(0.0033), 't_ket_post': -2.6383, 'p_ket_post': 0.0298, 'q_ket_post': np.float64(0.0298)}
all 13 sessions: {'pre': np.float64(0.832), 'ket': np.float64(0.007), 'post': np.float64(0.445)} {'t_pre_ket': 6.0402, 'p_pre_ket': 0.0001, 'q_pre_ket': np.float64(0.0001), 't_ket_post': -2.4253, 'p_ket_post': 0.032, 'q_ket_post': np.float64(0.032)}
fig, ax = plt.subplots(1, 2, figsize=(8, 3))
for j, first in enumerate((True, False)):
for k, col in (('pre', 'C0'), ('ket', 'C3'), ('post', 'C4')):
tr = np.array([m.persist_projection(R, 'M', first, 'ket_C', 'post_C')[k] for R in POP])
mu, ci = tr.mean(0), 1.96 * tr.std(0) / np.sqrt(len(tr))
ax[j].plot(t, mu, color=col, label=k); ax[j].fill_between(t, mu - ci, mu + ci, color=col, alpha=.3)
ax[j].set_title('first puff' if first else 'all 8 puffs'); ax[j].axvline(0, color='k', ls='--')
ax[j].set_xlabel('time from puff (s)'); ax[j].set_ylabel('projection (z)')
ax[0].legend(); plt.tight_layout(); plt.show()
Result: reproduced, resting on the assumption about the epochs. Under the main reading (per-neuron variance, first puff, trials 23 to 32 and 35 to 40) the projection 1 to 2 s after the puff is 0.867 ± 0.059 before, -0.002 ± 0.133 during and 0.515 ± 0.219 after ketamine, each within its tolerance; pre against ketamine t = 5.48, corrected p = 0.0005; post against ketamine t = -2.54, corrected p = 0.029. The code's pooled variance gives the same picture (0.855, 0.019, 0.513; corrected p = 0.0004 and 0.026). With the Methods' trials (22 to 32, 34 to 40) the fall on ketamine holds (corrected p = 0.0008) but the recovery test is not significant with the first puff (corrected p = 0.17; 0.017 when all 8 puffs are averaged). Averaged within mouse (9 mice) the two tests give corrected p = 0.0033 and 0.030; with all 13 sessions 0.0001 and 0.032.
fig6f-mouse-fast-activity-unchanged¶
In mice, activity along the fast dimension (0 to 70 ms) is not changed by ketamine.
Fast dimension built from 0 to 70 ms after each puff (all 8 puffs of the 20 preinfusion series) against the 1 s before it (Methods, 'reference window (average activity in the 1s before each puff)'; the code uses the 70 ms before, a variant); projection averaged over all puffs, z-scored, quantified 0 to 70 ms; epochs as in fig6d.
cid = 'fig6f-mouse-fast-activity-unchanged'
R6F = {}
for var in 'MC':
for ref in ('1s', '70'):
for ename, (ket, post) in m.EPOCH_READINGS.items():
R6F[f"var {var}, reference {ref}, {ename}"] = (var, ref, ket, post)
tab6f = run_readings_3way(m.fast_table, R6F, POP)
print(tab6f.to_string())
MAIN6F = ('M', '1s', 'ket_C', 'post_C')
d6f = m.fast_table(POP, *MAIN6F)
print('\nmain reading, per session:'); print(d6f.to_string())
s6f = m.three_way_tests(d6f)
print({k: round(v, 5) for k, v in s6f.items()})
for (qname, v, tol), col in zip(expected(cid), ('pre', 'ket', 'post')):
mu, se = m.mean_sem(d6f[col])
print(f'{qname}: {mu:.3f} ± {se:.3f} (expected {v} ± {tol}; within: {within(mu, v, tol)})')
for a, b in (('pre', 'ket'), ('post', 'ket')):
diff = d6f[a] - d6f[b]
print(f'{a} - {b}: mean {diff.mean():.3f}, 95% CI {diff.mean() - 1.96 * diff.std(ddof=1) / np.sqrt(len(diff)):.3f} '
f'to {diff.mean() + 1.96 * diff.std(ddof=1) / np.sqrt(len(diff)):.3f}, Cohen dz {diff.mean() / diff.std(ddof=1):.2f}')
dm = m.by_mouse(d6f, ['pre', 'ket', 'post'])
print('averaged within mouse (9 mice):', {k: round(v, 4) for k, v in m.three_way_tests(dm).items()})
d13 = m.fast_table(Rs, *MAIN6F)
print('all 13 sessions:', {k: round(v, 4) for k, v in m.three_way_tests(d13).items()})
reading pre ket post t_pre_ket p_pre_ket q_pre_ket t_ket_post p_ket_post q_ket_post
0 var M, reference 1s, Methods (22-32, 34-40) 4.226 ± 0.097 4.439 ± 0.117 4.036 ± 0.130 -1.7489 0.1109 0.1109 3.5638 0.0051 0.0103
1 var M, reference 1s, code / Table S3 (23-32, 35-40) 4.226 ± 0.097 4.279 ± 0.137 3.957 ± 0.130 -0.4332 0.6741 0.6741 4.0718 0.0022 0.0045
2 var M, reference 70, Methods (22-32, 34-40) 4.195 ± 0.105 4.419 ± 0.123 4.049 ± 0.120 -1.8289 0.0973 0.0973 3.2877 0.0082 0.0164
3 var M, reference 70, code / Table S3 (23-32, 35-40) 4.195 ± 0.105 4.259 ± 0.141 3.966 ± 0.120 -0.5348 0.6045 0.6045 3.7333 0.0039 0.0078
4 var C, reference 1s, Methods (22-32, 34-40) 4.470 ± 0.129 4.622 ± 0.137 4.245 ± 0.171 -0.8816 0.3987 0.3987 3.4567 0.0062 0.0123
5 var C, reference 1s, code / Table S3 (23-32, 35-40) 4.470 ± 0.129 4.499 ± 0.180 4.179 ± 0.172 -0.1482 0.8851 0.8851 3.3375 0.0075 0.0150
6 var C, reference 70, Methods (22-32, 34-40) 4.433 ± 0.132 4.619 ± 0.132 4.231 ± 0.159 -1.1428 0.2797 0.2797 3.4602 0.0061 0.0122
7 var C, reference 70, code / Table S3 (23-32, 35-40) 4.433 ± 0.132 4.496 ± 0.171 4.161 ± 0.160 -0.3525 0.7318 0.7318 3.4743 0.0060 0.0120
main reading, per session:
session mouse pre ket post
0 sub-0_20230901 0 4.1139 4.2875 3.8410
1 sub-1_20230218 1 4.4091 4.5214 4.1380
2 sub-1_20230226 1 4.7271 4.6458 4.4360
3 sub-2_20230902 2 4.7947 5.1908 4.6432
4 sub-3_20230218 3 4.2579 3.4443 3.0446
5 sub-3_20230225 3 3.9577 3.7445 3.4471
6 sub-4_20230907 4 4.1599 4.3555 4.2030
7 sub-6_20230823 6 3.9252 4.3429 3.6183
8 sub-7_20220831 7 3.7797 4.5006 3.9631
9 sub-8_20220228 8 3.9184 3.8484 4.0680
10 sub-9_20230908 9 4.4377 4.1886 4.1267
{'t_pre_ket': -0.43321, 'p_pre_ket': 0.67406, 'q_pre_ket': np.float64(0.67406), 't_ket_post': 4.07184, 'p_ket_post': 0.00224, 'q_ket_post': np.float64(0.00449)}
fast projection preinfusion: 4.226 ± 0.097 (expected 4.48 ± 0.129; within: False)
fast projection ketamine: 4.279 ± 0.137 (expected 4.59 ± 0.185; within: False)
fast projection postinfusion: 3.957 ± 0.130 (expected 4.32 ± 0.184; within: False)
pre - ket: mean -0.054, 95% CI -0.296 to 0.189, Cohen dz -0.13
post - ket: mean -0.322, 95% CI -0.477 to -0.167, Cohen dz -1.23
averaged within mouse (9 mice): {'t_pre_ket': -0.9688, 'p_pre_ket': 0.361, 'q_pre_ket': np.float64(0.361), 't_ket_post': 3.3459, 'p_ket_post': 0.0101, 'q_ket_post': np.float64(0.0203)}
all 13 sessions: {'t_pre_ket': 0.5883, 'p_pre_ket': 0.5672, 'q_pre_ket': np.float64(0.5672), 't_ket_post': 3.428, 'p_ket_post': 0.005, 'q_ket_post': np.float64(0.01)}
Result: not reproduced. This is a claim of absence, judged on the paper's two tests, and one of them finds a difference. Pre against ketamine is not significant in any reading (main reading t = -0.43, corrected p = 0.67; difference -0.054, 95% CI -0.296 to 0.189, dz = -0.13), as the paper reports. Ketamine against post is significant in every reading (main reading t = 4.07, corrected p = 0.0045; post below ketamine by 0.322, dz = -1.23), where the paper reports corrected p = 0.066; the criterion asks both tests to be non-significant. The main reading's means (4.226, 4.279, 3.957) are below the paper's by more than its s.e.m.; the code's reading (pooled variance, 70-ms reference) gives 4.433, 4.496 and 4.161, within them.
fig5e-mouse-emotion-dimension-late-reduced¶
Along the mouse emotion-like dimension, ketamine leaves the rapid response to the first puff unchanged (0 to 0.25 s) but removes the late response (2 to 3 s).
Emotion-like dimension per session from the puff-modulated cells: target window 1 s starting 5 s after the last puff onset of each preinfusion series against the 750 ms ending 250 ms before the first puff, orthogonalised (QR) against a fast dimension (200 ms after against 200 ms before puffs 5 to 8). Signal: the projection of the trial-averaged series, baseline the 250 ms before the first puff; early 0 to 0.25 s, late 2 to 3 s after the first puff; paired t-tests with Benjamini-Hochberg over the two. Readings: variance M or C; fast dimension from puffs 5 to 8 or from the last 4 preinfusion series; ketamine trials 22 to 32 or 23 to 32. A further variant uses the code's target window, bins 400:500 of a window that, from the code's plotting offsets, starts 1 s before the last puff (3 to 4 s after the last puff onset). This offset is an inference: the code that builds the window is not in the record. The analysis.py lines it rests on are printed below the table. The criterion's 5 to 6 s target is the main reading and decides the outcome.
cid = 'fig5e-mouse-emotion-dimension-late-reduced'
rows = []
TAB5E = {}
for target in ('5-6', '3-4'):
for var in 'MC':
for eye in ('puffs5-8', 'last4series'):
for ket in ('ket_C', 'ket_M'):
d, sig, tt = m.affective_table(POP, target, var, eye, ket)
TAB5E[(target, var, eye, ket)] = (d, sig, tt)
te, pe = m.paired(d['pre_early'], d['ket_early']); tl, pl = m.paired(d['pre_late'], d['ket_late'])
q = m.bh([pe, pl])
r = dict(target=target + ' s', var=var, fast_dim=eye, ket=ket)
for c in ('pre_early', 'ket_early', 'pre_late', 'ket_late'):
mu, se = m.mean_sem(d[c]); r[c] = f'{mu:.3f} ± {se:.3f}'
r.update(t_early=te, q_early=q[0], t_late=tl, q_late=q[1])
rows.append(r)
print(pd.DataFrame(rows).to_string())
MAIN5E = ('5-6', 'M', 'puffs5-8', 'ket_C')
d5e = TAB5E[MAIN5E][0]
print('\nmain reading, per session:'); print(d5e.to_string())
te, pe = m.paired(d5e['pre_early'], d5e['ket_early']); tl, pl = m.paired(d5e['pre_late'], d5e['ket_late'])
print(f'early t = {te:.3f}, late t = {tl:.3f}, BH-corrected p = {np.round(m.bh([pe, pl]), 4)}')
for (qname, v, tol), col in zip(expected(cid), ('pre_early', 'ket_early', 'pre_late', 'ket_late')):
mu, se = m.mean_sem(d5e[col])
print(f'{qname}: {mu:.3f} ± {se:.3f} (expected {v} ± {tol}; within: {within(mu, v, tol)})')
CODE5E = ('3-4', 'C', 'puffs5-8', 'ket_C')
dc = TAB5E[CODE5E][0]
te, pe = m.paired(dc['pre_early'], dc['ket_early']); tl, pl = m.paired(dc['pre_late'], dc['ket_late'])
print(f'\ncode-window variant (3-4 s, pooled variance): early t = {te:.3f}, late t = {tl:.3f}, BH p = {np.round(m.bh([pe, pl]), 4)}')
for (qname, v, tol), col in zip(expected(cid), ('pre_early', 'ket_early', 'pre_late', 'ket_late')):
mu, se = m.mean_sem(dc[col]); print(f' {qname}: {mu:.3f} ± {se:.3f} (expected {v} ± {tol}; within: {within(mu, v, tol)})')
dm = m.by_mouse(d5e, ['pre_early', 'ket_early', 'pre_late', 'ket_late'])
print('main reading averaged within mouse: late t, p =', np.round(m.paired(dm['pre_late'], dm['ket_late']), 4),
'; early t, p =', np.round(m.paired(dm['pre_early'], dm['ket_early']), 4))
d13, _, _ = m.affective_table(Rs, *MAIN5E)
print('main reading, 13 sessions: late t, p =', np.round(m.paired(d13['pre_late'], d13['ket_late']), 4),
'; early t, p =', np.round(m.paired(d13['pre_early'], d13['ket_early']), 4))
target var fast_dim ket pre_early ket_early pre_late ket_late t_early q_early t_late q_late
0 5-6 s M puffs5-8 ket_C 0.899 ± 0.300 0.362 ± 0.160 0.343 ± 0.143 0.009 ± 0.076 1.9543 0.0859 1.9048 0.0859
1 5-6 s M puffs5-8 ket_M 0.899 ± 0.300 0.405 ± 0.159 0.343 ± 0.143 -0.001 ± 0.070 1.8864 0.0886 1.9620 0.0886
2 5-6 s M last4series ket_C 0.535 ± 0.382 0.292 ± 0.169 0.361 ± 0.146 0.014 ± 0.075 0.7480 0.4717 1.9552 0.1581
3 5-6 s M last4series ket_M 0.535 ± 0.382 0.321 ± 0.168 0.361 ± 0.146 0.004 ± 0.071 0.6760 0.5144 2.0025 0.1462
4 5-6 s C puffs5-8 ket_C 0.328 ± 0.236 0.258 ± 0.124 0.428 ± 0.165 -0.009 ± 0.053 0.3384 0.7421 2.4170 0.0725
5 5-6 s C puffs5-8 ket_M 0.328 ± 0.236 0.295 ± 0.122 0.428 ± 0.165 -0.010 ± 0.057 0.1599 0.8762 2.4671 0.0665
6 5-6 s C last4series ket_C 0.216 ± 0.359 0.249 ± 0.129 0.460 ± 0.178 -0.002 ± 0.053 -0.1100 0.9146 2.3879 0.0762
7 5-6 s C last4series ket_M 0.216 ± 0.359 0.278 ± 0.127 0.460 ± 0.178 -0.004 ± 0.058 -0.2080 0.8394 2.4383 0.0699
8 3-4 s M puffs5-8 ket_C 1.063 ± 0.326 0.472 ± 0.168 0.463 ± 0.157 0.005 ± 0.072 2.4211 0.0360 2.6519 0.0360
9 3-4 s M puffs5-8 ket_M 1.063 ± 0.326 0.519 ± 0.173 0.463 ± 0.157 -0.010 ± 0.057 2.3961 0.0376 2.7543 0.0376
10 3-4 s M last4series ket_C 0.732 ± 0.370 0.416 ± 0.170 0.455 ± 0.160 0.008 ± 0.072 1.0863 0.3028 2.5507 0.0577
11 3-4 s M last4series ket_M 0.732 ± 0.370 0.453 ± 0.173 0.455 ± 0.160 -0.004 ± 0.059 1.0136 0.3347 2.6124 0.0519
12 3-4 s C puffs5-8 ket_C 0.608 ± 0.286 0.308 ± 0.132 0.534 ± 0.177 0.004 ± 0.043 1.3035 0.2216 2.7736 0.0393
13 3-4 s C puffs5-8 ket_M 0.608 ± 0.286 0.343 ± 0.132 0.534 ± 0.177 -0.006 ± 0.038 1.2102 0.2540 2.8257 0.0360
14 3-4 s C last4series ket_C 0.435 ± 0.370 0.295 ± 0.140 0.538 ± 0.185 0.009 ± 0.042 0.4739 0.6457 2.6848 0.0458
15 3-4 s C last4series ket_M 0.435 ± 0.370 0.320 ± 0.139 0.538 ± 0.185 0.001 ± 0.039 0.4007 0.6970 2.7061 0.0442
main reading, per session:
session mouse pre_early pre_late ket_early ket_late
0 sub-0_20230901 0 2.0130 1.4125e+00 0.1120 -0.0680
1 sub-1_20230218 1 0.8847 2.9348e-01 -0.5212 -0.1733
2 sub-1_20230226 1 -0.5415 1.5574e-01 0.2181 0.2609
3 sub-2_20230902 2 0.3315 -1.9428e-05 0.3088 0.1390
4 sub-3_20230218 3 0.4905 5.7155e-02 0.3293 -0.1071
5 sub-3_20230225 3 0.4007 -2.1249e-01 -0.0124 -0.2308
6 sub-4_20230907 4 0.9958 6.2272e-01 1.4201 0.3710
7 sub-6_20230823 6 1.1002 7.9579e-01 0.4120 -0.1190
8 sub-7_20220831 7 -0.3734 -3.2727e-01 -0.0661 0.1113
9 sub-8_20220228 8 1.4949 4.0503e-01 0.5541 0.3693
10 sub-9_20230908 9 3.0886 5.7467e-01 1.2240 -0.4544
early t = 1.954, late t = 1.905, BH-corrected p = [0.0859 0.0859]
0 to 0.25 s, preinfusion: 0.899 ± 0.300 (expected 0.601 ± 0.28; within: False)
0 to 0.25 s, ketamine: 0.362 ± 0.160 (expected 0.332 ± 0.138; within: True)
2 to 3 s, preinfusion: 0.343 ± 0.143 (expected 0.546 ± 0.174; within: False)
2 to 3 s, ketamine: 0.009 ± 0.076 (expected 0.005 ± 0.043; within: True)
code-window variant (3-4 s, pooled variance): early t = 1.304, late t = 2.774, BH p = [0.2216 0.0393]
0 to 0.25 s, preinfusion: 0.608 ± 0.286 (expected 0.601 ± 0.28; within: True)
0 to 0.25 s, ketamine: 0.308 ± 0.132 (expected 0.332 ± 0.138; within: True)
2 to 3 s, preinfusion: 0.534 ± 0.177 (expected 0.546 ± 0.174; within: True)
2 to 3 s, ketamine: 0.004 ± 0.043 (expected 0.005 ± 0.043; within: True)
main reading averaged within mouse: late t, p = [1.8193 0.1064] ; early t, p = [2.0742 0.0718]
main reading, 13 sessions: late t, p = [2.3243 0.0385] ; early t, p = [2.6415 0.0215]
# The analysis.py lines (Zenodo 14876694) from which the 3 to 4 s variant is inferred
CODE = m.fetch_analysis_py()
print('analysis.py fetched,', len(CODE), 'lines, md5', m.CODE_NPX_MD5)
print('\n(a) get_windowed_puff_psth: a window starts int(1./bin_t) = 100 bins (1 s) before its puff')
m.show_code(CODE, 432, 432)
print('\n(b) get_affective_projections: the target window is bins 400:500 of the post-series window')
m.show_code(CODE, 1436, 1438)
print('\n(c) every line that mentions post_psth; none builds the array (the driver that does is not in the record)')
for i, line in enumerate(CODE, 1):
if 'post_psth' in line:
print(f'{i:5d}: {line.strip()[:150]}')
print('\n(d) plot_affective_projection_late: the post-series window from bin start = 500 is drawn from 25 s after the first puff')
m.show_code(CODE, 1482, 1484); m.show_code(CODE, 1490, 1494)
print('\n(e) the model of fig5g (time 0 is 1 s before the first puff, 10-ms bins): its second segment starts at bin 300*8 + 200, i.e. 25 s')
m.show_code(CODE, 2765, 2767); m.show_code(CODE, 2833, 2838)
print('\n(f) plot_affective_firstpuff_persist joins the post-series window from bin 400, but reads only bins 25:50 and 225:325, inside the series part')
m.show_code(CODE, 1586, 1595)
start_bin, drawn_at, last_puff = 500, 25.0, 7 * 3.0
win0 = drawn_at - start_bin * m.BIN
print(f'\nFrom (d) and (e): bin 0 of the post-series window lies {win0:.1f} s after the first puff, {last_puff - win0:.1f} s before the '
f'last puff onset ({last_puff:.0f} s), as (a) gives for a window built on the last puff; bins 400:500 then run '
f'{win0 + 4 - last_puff:.1f} to {win0 + 5 - last_puff:.1f} s after the last puff onset. The Methods and the criterion give 5 to 6 s.')
analysis.py fetched, 2958 lines, md5 6cc7aa5059737c3bd753aad1f116019d
(a) get_windowed_puff_psth: a window starts int(1./bin_t) = 100 bins (1 s) before its puff
432: sta = np.searchsorted(bins, data._events['puff_times']) - int(1./bin_t)
(b) get_affective_projections: the target window is bins 400:500 of the post-series window
1436: # cd_affective captures difference from post puff series to pre
1437: cd_affective = compute_cd_cov(post_psth_combined[r_cells,:,400:500].mean(-1).reshape((smfr.shape[0],-1)),
1438: psth_combined[r_cells,:,0,0:75].mean(-1).reshape((smfr.shape[0],-1)))
(c) every line that mentions post_psth; none builds the array (the driver that does is not in the record)
1408: def get_affective_projections(puff_psth, puff_psth_ket, puff_post_psth, puff_post_psth_ket, puff_post_psth_recovery, all_recording_names, recording_na
1420: post_psth_combined = puff_post_psth[:,:,:]
1437: cd_affective = compute_cd_cov(post_psth_combined[r_cells,:,400:500].mean(-1).reshape((smfr.shape[0],-1)),
1457: smfr = puff_post_psth[r_cells,:,:]
1463: smfr = puff_post_psth_ket[r_cells,:,:]
(d) plot_affective_projection_late: the post-series window from bin start = 500 is drawn from 25 s after the first puff
1482: start = 500
1483: ket_signal = np.concatenate([np.stack(ket_proj_normed).reshape((nsession, 8*300))[:,75:],
1484: np.stack(late_proj_normed_ket)[:,start:50*100]], axis=1) - np.expand_dims(np.stack(ket_proj_normed).reshape((nsession,
1490:
1491: sig = np.array([sosfiltfilt(sos,k) for k in ket_signal[:,3*8*100-75:]])
1492: ste = np.std(sig,0)/np.sqrt(sig.shape[0])
1493: plt.fill_between(np.arange(25, 25+45, 0.01), sig.mean(0)+1.96*ste, sig.mean(0)-1.96*ste, alpha=0.3, color=color_dict['Ket'], linewidth=0)
1494: plt.plot(np.arange(25, 25+45, 0.01), sig.mean(0), color=color_dict['Ket'])
(e) the model of fig5g (time 0 is 1 s before the first puff, 10-ms bins): its second segment starts at bin 300*8 + 200, i.e. 25 s
2765: bin_t = 0.01
2766: inputs = np.hstack([np.concatenate([np.zeros(100), np.ones(25), np.zeros(175)]) for _ in range(8)])
2767: inputs = np.hstack([inputs, np.zeros(int(47/bin_t))])
2833: def theoretical_persist(x, S, tB, tP, *kwargs):
2834: traces = []
2835: for c in kwargs:
2836: output = model_system(S, tB, tP*c)
2837: traces.append(np.hstack([output[75:300*8], output[(300*8)+200:(300*8)+int(49/bin_t)]]))
2838: return np.hstack(traces)
(f) plot_affective_firstpuff_persist joins the post-series window from bin 400, but reads only bins 25:50 and 225:325, inside the series part
1586: ket_signal = np.concatenate([np.stack(ket_proj_normed).reshape((nsession, 8*300))[:,75:],
1587: np.stack(late_proj_normed_ket)[:,400:50*100]], axis=1) - np.expand_dims(np.stack(ket_proj_normed).reshape((nsession, 8*30
1588:
1589: pre_signal = np.concatenate([np.stack(proj_normed).reshape((nsession, 8*300))[:,75:],
1590: np.stack(late_proj_normed)[:,400:50*100]], axis=1) - np.expand_dims(np.stack(proj_normed).reshape((nsession, 8*300))[:,75
1591:
1592: df = pd.DataFrame({"Projection (AU)": np.hstack([pre_signal[:,25:50].mean(-1),
1593: ket_signal[:,25:50].mean(-1),
1594: pre_signal[:,225:325].mean(-1),
1595: ket_signal[:,225:325].mean(-1)]),
From (d) and (e): bin 0 of the post-series window lies 20.0 s after the first puff, 1.0 s before the last puff onset (21 s), as (a) gives for a window built on the last puff; bins 400:500 then run 3.0 to 4.0 s after the last puff onset. The Methods and the criterion give 5 to 6 s.
Result: not reproduced. The outcome is decided under the criterion's target window (1 s starting 5 s after the last puff). Late activity falls from 0.343 ± 0.143 to 0.009 ± 0.076, but the late test gives t = 1.90, corrected p = 0.086, and the early test t = 1.95, corrected p = 0.086; the preinfusion means (0.899 early, 0.343 late) are outside their tolerances. No reading with the 5 to 6 s window gives a corrected late p below 0.05 (0.066 to 0.158). The 3 to 4 s window is not in the criterion and does not decide the outcome. It is inferred from the printed code lines: the target is bins 400:500 of a post-series window, and the plotting code and the fig5g model place bin 500 of that window 25 s after the first puff, so bin 0 lies 1 s before the last puff onset and bins 400:500 run 3 to 4 s after it. With that window and the pooled variance the paper's values return: early 0.608 and 0.308 (corrected p = 0.22), late 0.534 and 0.004 (t = 2.77, corrected p = 0.039), all four means within tolerance. With all 13 sessions the main reading gives uncorrected p = 0.039 for the late test and 0.022 for the early test.
fig5d-mouse-emotion-dimension-accumulates¶
In mice, activity along the emotion-like dimension rises and saturates across the eight-puff series and decays over tens of seconds after it; on ketamine it does not accumulate.
Signal of fig5e (each session on its own dimension, baseline the 250 ms before the first puff), mean and 95% interval (1.96 s.e.m.) over the 11 sessions; the 23 to 25 s stretch around the training window is left out, as in the code. Criterion values: the 1 s before the eighth puff (20 to 21 s) minus the 1 s before the second puff (2 to 3 s), and the mean 10 to 20 s after the last puff (31 to 41 s after the first).
cid = 'fig5d-mouse-emotion-dimension-accumulates'
rows = []
for key, (d, sig, tt) in TAB5E.items():
r = dict(target=key[0], var=key[1], fast_dim=key[2], ket=key[3])
for k in ('pre', 'ket'):
S = sig[k]
w = lambda a, b: np.nanmean(S[:, (tt >= a) & (tt < b)], 1)
acc = w(20, 21) - w(2, 3); late = w(31, 41)
r[f'{k}_accum'] = acc.mean(); r[f'{k}_accum_ci'] = 1.96 * acc.std() / np.sqrt(len(acc))
r[f'{k}_10_20s'] = late.mean(); r[f'{k}_10_20s_ci'] = 1.96 * late.std() / np.sqrt(len(late))
r['pass'] = (r['pre_accum'] > 0) and (r['pre_10_20s'] - r['pre_10_20s_ci'] > 0) and (r['ket_accum'] <= r['ket_accum_ci'])
rows.append(r)
tab5d = pd.DataFrame(rows)
print(tab5d.to_string())
d, sig, tt = TAB5E[MAIN5E]
for k in ('pre', 'ket'):
S = sig[k]
acc = np.nanmean(S[:, (tt >= 20) & (tt < 21)], 1) - np.nanmean(S[:, (tt >= 2) & (tt < 3)], 1)
print(k, 'per session accumulation:', np.round(acc, 3))
fig, ax = plt.subplots(figsize=(9, 3))
gapless = ~np.isnan(sig['pre'][0])
for k, col in (('pre', 'C0'), ('ket', 'C3')):
S = sig[k][:, gapless]; mu = S.mean(0); ci = 1.96 * S.std(0) / np.sqrt(len(S))
ax.plot(tt[gapless][tt[gapless] < 23], mu[tt[gapless] < 23], color=col, label=k)
ax.plot(tt[gapless][tt[gapless] >= 25], mu[tt[gapless] >= 25], color=col)
ax.fill_between(tt[gapless], mu - ci, mu + ci, color=col, alpha=.3)
for i in range(8):
ax.axvline(3 * i, color='k', ls='--', alpha=.3)
ax.axhline(0, color='k'); ax.set_xlabel('time from first puff (s)'); ax.set_ylabel('emotion-like dim. (a.u.)'); ax.legend()
plt.show()
print('plotted session means range from', min(sig[k][:, gapless].mean(0).min() for k in sig).round(3), 'to',
max(sig[k][:, gapless].mean(0).max() for k in sig).round(3), '; all within the axes')
# all 13 sessions (mouse 5 included), main reading
_, sig13, tt13 = m.affective_table(Rs, *MAIN5E)
for k in ('pre', 'ket'):
S = sig13[k]; w = lambda a, b: np.nanmean(S[:, (tt13 >= a) & (tt13 < b)], 1)
acc = w(20, 21) - w(2, 3); late = w(31, 41)
print(f'13 sessions, {k}: accumulation {acc.mean():.3f} ± {1.96 * acc.std() / np.sqrt(len(acc)):.3f}, '
f'10-20 s after the last puff {late.mean():.3f} ± {1.96 * late.std() / np.sqrt(len(late)):.3f}')
target var fast_dim ket pre_accum pre_accum_ci pre_10_20s pre_10_20s_ci ket_accum ket_accum_ci ket_10_20s ket_10_20s_ci pass 0 5-6 M puffs5-8 ket_C 0.8458 0.1697 0.8381 0.1692 -0.0477 0.1136 0.0453 0.1161 True 1 5-6 M puffs5-8 ket_M 0.8458 0.1697 0.8381 0.1692 -0.0199 0.1033 0.0395 0.1007 True 2 5-6 M last4series ket_C 0.8694 0.1799 0.8817 0.1852 -0.0568 0.1161 0.0426 0.1178 True 3 5-6 M last4series ket_M 0.8694 0.1799 0.8817 0.1852 -0.0285 0.1078 0.0363 0.1044 True 4 5-6 C puffs5-8 ket_C 1.0564 0.2316 1.0625 0.2337 -0.0596 0.0980 0.0130 0.0817 True 5 5-6 C puffs5-8 ket_M 1.0564 0.2316 1.0625 0.2337 -0.0403 0.0930 0.0141 0.0868 True 6 5-6 C last4series ket_C 1.0572 0.2359 1.0950 0.2507 -0.0655 0.0996 0.0131 0.0851 True 7 5-6 C last4series ket_M 1.0572 0.2359 1.0950 0.2507 -0.0465 0.0953 0.0121 0.0911 True 8 3-4 M puffs5-8 ket_C 0.8205 0.2124 0.7307 0.1636 -0.0491 0.0987 0.0295 0.1239 True 9 3-4 M puffs5-8 ket_M 0.8205 0.2124 0.7307 0.1636 -0.0115 0.0897 0.0327 0.0963 True 10 3-4 M last4series ket_C 0.8581 0.2515 0.7538 0.1711 -0.0567 0.0995 0.0262 0.1260 True 11 3-4 M last4series ket_M 0.8581 0.2515 0.7538 0.1711 -0.0195 0.0912 0.0323 0.0982 True 12 3-4 C puffs5-8 ket_C 1.0210 0.2739 0.9048 0.2160 -0.0551 0.0760 0.0164 0.0808 True 13 3-4 C puffs5-8 ket_M 1.0210 0.2739 0.9048 0.2160 -0.0240 0.0747 0.0215 0.0785 True 14 3-4 C last4series ket_C 1.0451 0.3054 0.9247 0.2250 -0.0612 0.0766 0.0154 0.0832 True 15 3-4 C last4series ket_M 1.0451 0.3054 0.9247 0.2250 -0.0314 0.0748 0.0215 0.0795 True pre per session accumulation: [0.667 0.703 1.126 1.002 0.354 0.88 0.4 0.741 1.052 1.291 1.087] ket per session accumulation: [-0.172 -0.248 -0.337 -0.075 0.161 0.349 0.07 0.111 -0.147 -0.176 -0.061]
plotted session means range from -0.408 to 2.219 ; all within the axes 13 sessions, pre: accumulation 0.779 ± 0.185, 10-20 s after the last puff 0.769 ± 0.177 13 sessions, ket: accumulation -0.044 ± 0.099, 10-20 s after the last puff 0.048 ± 0.098
Result: reproduced. Before ketamine the activity in the 1 s before the eighth puff exceeds that before the second by 0.846 (95% interval ± 0.170), every session positive (0.354 to 1.291), and 10 to 20 s after the last puff it is 0.838 ± 0.169, above zero; on ketamine the accumulation is -0.048 ± 0.114. All 16 readings pass the three conditions. With all 13 sessions the accumulation is 0.779 ± 0.185 before and -0.044 ± 0.099 on ketamine, and 0.769 ± 0.177 10 to 20 s after the last puff.
fig6h-mouse-timescale-shorter¶
In mice, the baseline intrinsic timescale of the persistent population is shorter during ketamine (0.14 s) than before (0.45 s).
The 30 s before each series projected on the persistent dimension (all puff-modulated cells, session z-scored rates); autocorrelation fitted with a*exp(-x/b)+c over 2.5 s from a = 1, b = 0.2, c = 0. Series: the code's 11 to 20 and 23 to 32 (Table S3: 10 and 10; the Methods do not say which) and, as a variant, 1 to 20 and 22 to 32. Readings: fit from lag 0 (Methods, "from the peak") or lag one bin (code); each series fitted and the timescales averaged (Methods, Table S3) or the series concatenated (code); variance M or C.
cid = 'fig6h-mouse-timescale-shorter'
SER6H = {'code 11-20 / 23-32': (np.arange(10, 20), np.arange(22, 32)),
'Methods epochs 1-20 / 22-32': (np.arange(0, 20), np.arange(21, 32))}
rows = []
with warnings.catch_warnings():
warnings.simplefilter('ignore')
for var in 'MC':
for sname, (sp, sk) in SER6H.items():
for lag in (0, 1):
for mode in ('per-series', 'concat'):
d, nf = m.timescale_table(POP, var, sp, sk, lag, mode)
tt_, pp_ = m.paired(d['pre'], d['ket'])
r = dict(var=var, series=sname, lag=lag, mode=mode, failed_fits=nf)
for k in ('pre', 'ket'):
mu, se = m.mean_sem(d[k]); r[k] = f'{mu:.3f} ± {se:.3f}'
r.update(t=tt_, p=pp_); rows.append(r)
print(pd.DataFrame(rows).to_string())
d6h, nf = m.timescale_table(POP, 'M', *SER6H['code 11-20 / 23-32'], 0, 'per-series')
print('\nmain reading (M, code series, lag 0, per series), per session:'); print(d6h.to_string())
print('failed fits:', nf, '; paired t, p =', np.round(m.paired(d6h['pre'], d6h['ket']), 4))
for (qname, v, tol), col in zip(expected(cid), ('pre', 'ket')):
mu, se = m.mean_sem(d6h[col]); print(f'{qname}: {mu:.3f} ± {se:.3f} s (expected {v} ± {tol}; within: {within(mu, v, tol)})')
dm = m.by_mouse(d6h, ['pre', 'ket'])
print('averaged within mouse:', np.round(m.paired(dm['pre'], dm['ket']), 4), '; sessions with ketamine shorter:', int((d6h.ket < d6h.pre).sum()), 'of', len(d6h))
with warnings.catch_warnings():
warnings.simplefilter('ignore')
d13, nf13 = m.timescale_table(Rs, 'M', *SER6H['code 11-20 / 23-32'], 0, 'per-series')
print('all 13 sessions, main reading:', {k: round(m.mean_sem(d13[k])[0], 3) for k in ('pre', 'ket')},
'failed fits', nf13, '; t, p =', np.round(m.paired(d13['pre'], d13['ket']), 4))
print('mouse 5 sessions:'); print(d13[d13.mouse == '5'].to_string())
ok13 = (d13.pre <= 2.5) & (d13.ket <= 2.5)
print(f'13 sessions without those whose fitted timescale exceeds the 2.5-s fitting window ({int((~ok13).sum())} dropped): '
f"pre {d13.pre[ok13].mean():.3f}, ket {d13.ket[ok13].mean():.3f}; t, p =", np.round(m.paired(d13.pre[ok13], d13.ket[ok13]), 4))
var series lag mode failed_fits pre ket t p
0 M code 11-20 / 23-32 0 per-series 0 0.437 ± 0.077 0.162 ± 0.079 3.5067 0.0057
1 M code 11-20 / 23-32 0 concat 0 0.360 ± 0.056 0.119 ± 0.062 3.6576 0.0044
2 M code 11-20 / 23-32 1 per-series 0 0.481 ± 0.082 0.189 ± 0.083 3.4419 0.0063
3 M code 11-20 / 23-32 1 concat 0 0.388 ± 0.058 0.119 ± 0.065 3.7140 0.0040
4 M Methods epochs 1-20 / 22-32 0 per-series 0 0.391 ± 0.061 0.179 ± 0.077 2.8856 0.0162
5 M Methods epochs 1-20 / 22-32 0 concat 0 0.323 ± 0.047 0.120 ± 0.062 3.1552 0.0102
6 M Methods epochs 1-20 / 22-32 1 per-series 0 0.440 ± 0.071 0.223 ± 0.079 2.5173 0.0305
7 M Methods epochs 1-20 / 22-32 1 concat 0 0.347 ± 0.049 0.156 ± 0.069 2.4291 0.0355
8 C code 11-20 / 23-32 0 per-series 0 0.493 ± 0.082 0.192 ± 0.073 3.2117 0.0093
9 C code 11-20 / 23-32 0 concat 0 0.391 ± 0.060 0.131 ± 0.060 3.4696 0.0060
10 C code 11-20 / 23-32 1 per-series 0 0.557 ± 0.097 0.221 ± 0.077 2.9413 0.0148
11 C code 11-20 / 23-32 1 concat 0 0.417 ± 0.063 0.161 ± 0.063 2.8845 0.0163
12 C Methods epochs 1-20 / 22-32 0 per-series 0 0.428 ± 0.064 0.214 ± 0.069 2.6921 0.0226
13 C Methods epochs 1-20 / 22-32 0 concat 0 0.345 ± 0.049 0.152 ± 0.060 2.7146 0.0218
14 C Methods epochs 1-20 / 22-32 1 per-series 0 0.476 ± 0.076 0.249 ± 0.072 2.4001 0.0373
15 C Methods epochs 1-20 / 22-32 1 concat 0 0.365 ± 0.051 0.171 ± 0.066 2.3877 0.0381
main reading (M, code series, lag 0, per series), per session:
session mouse pre ket
0 sub-0_20230901 0 0.8512 0.0384
1 sub-1_20230218 1 0.4721 0.0357
2 sub-1_20230226 1 0.2679 0.0593
3 sub-2_20230902 2 0.0686 0.0345
4 sub-3_20230218 3 0.2464 0.2478
5 sub-3_20230225 3 0.4019 0.0350
6 sub-4_20230907 4 0.1170 0.0672
7 sub-6_20230823 6 0.4729 0.0366
8 sub-7_20220831 7 0.4571 0.1531
9 sub-8_20220228 8 0.9028 0.9607
10 sub-9_20230908 9 0.5518 0.1154
failed fits: 0 ; paired t, p = [3.5067 0.0057]
timescale preinfusion: 0.437 ± 0.077 s (expected 0.45 ± 0.07; within: True)
timescale ketamine: 0.162 ± 0.079 s (expected 0.14 ± 0.06; within: True)
averaged within mouse: [3.142 0.0138] ; sessions with ketamine shorter: 9 of 11
all 13 sessions, main reading: {'pre': np.float64(0.397), 'ket': np.float64(12.049)} failed fits 0 ; t, p = [-0.9782 0.3473]
mouse 5 sessions:
session mouse pre ket
7 sub-5_20230706 5 0.1105 0.0240
8 sub-5_20230822 5 0.2464 154.8256
13 sessions without those whose fitted timescale exceeds the 2.5-s fitting window (1 dropped): pre 0.410, ket 0.151; t, p = [3.5379 0.0046]
Result: reproduced. Main reading: 0.437 ± 0.077 s before and 0.162 ± 0.079 s during ketamine, both within tolerance; t = 3.51, p = 0.0057; 9 of 11 sessions shorter on ketamine; no fit failed. All 16 readings give p < 0.05 (0.004 to 0.038) with ketamine shorter; their preinfusion means range from 0.323 to 0.557 s and some fall outside the tolerance. Averaged within mouse, t = 3.14, p = 0.014. With all 13 sessions the test fails (t = -0.98, p = 0.35) because one ketamine fit of mouse 5 (sub-5 ses-20230822) gives 154.8 s, far beyond the 2.5-s fitting window; without that session the 12 give 0.410 against 0.151 s, t = 3.54, p = 0.0046.
fig6j-mouse-coupling-reduced¶
In mice, baseline pairwise correlation among the most persistent neurons halves during ketamine (0.49 of preinfusion).
Pearson correlations of the session z-scored rates in the 30 s before each series, among the puff-modulated cells whose persistent weight is above the 90th percentile of |w|. Readings: selection and average (Methods: |w| above the percentile and the mean of r; code: signed w above it and the median of |r|); correlations within each 30-s window averaged over series (Methods, Table S3: "Values averaged within session across trials") or over the concatenated windows (code); variance M or C; series of the code (1 to 10, 23 to 32; Table S3: 10 and 10) or of the Methods epochs (1 to 20, 22 to 32). The test is the code's: paired t-test of the normalised values (1 against ketamine/preinfusion).
cid = 'fig6j-mouse-coupling-reduced'
SER6J = {'code 1-10 / 23-32': ('pre_C', 'ket_C', np.arange(0, 10), np.arange(22, 32)),
'Methods epochs 1-20 / 22-32': ('pre_M', 'ket_M', np.arange(0, 20), np.arange(21, 32))}
def coupling_table(pop, var, rule, sname, mode):
ep_pre, ep_ket, sp, sk = SER6J[sname]
rows = []
for R in pop:
if mode == 'per-series':
v = R['coupling_series'][(var, rule)]
rows.append(dict(pre=np.nanmean(v[sp]), ket=np.nanmean(v[sk]), n_cells=R['coupling'][(var, rule, ep_pre, 'n_cells')]))
else:
rows.append(dict(pre=R['coupling'][(var, rule, ep_pre)], ket=R['coupling'][(var, rule, ep_ket)],
n_cells=R['coupling'][(var, rule, ep_pre, 'n_cells')]))
d = m._df(pop, rows); d['ratio'] = d.ket / d.pre
return d
rows = []
for var in 'MC':
for rule in 'MC':
for sname in SER6J:
for mode in ('per-series', 'concat'):
d = coupling_table(POP, var, rule, sname, mode)
mu, se = m.mean_sem(d.ratio); tt_, pp_ = m.paired(np.ones(len(d)), d.ratio)
rows.append(dict(var=var, rule=rule, series=sname, mode=mode, ratio=f'{mu:.3f} ± {se:.3f}', t=tt_, p=pp_,
p_raw_pre_vs_ket=m.paired(d.pre, d.ket)[1]))
print(pd.DataFrame(rows).to_string())
d6j = coupling_table(POP, 'M', 'M', 'code 1-10 / 23-32', 'per-series')
print('\nmain reading (M weights, |w| and mean r, code series, per series), per session:'); print(d6j.to_string())
mu, se = m.mean_sem(d6j.ratio)
print(f'ratio {mu:.3f} ± {se:.3f}; paired t, p =', np.round(m.paired(np.ones(len(d6j)), d6j.ratio), 5))
for (qname, v, tol) in expected(cid):
print(f'{qname}: {mu:.3f} (expected {v} ± {tol}; within: {within(mu, v, tol)})')
dc = coupling_table(POP, 'C', 'C', 'code 1-10 / 23-32', 'concat')
print('code reading (C weights, signed w, median |r|, concatenated):', np.round(m.mean_sem(dc.ratio), 3), np.round(m.paired(np.ones(len(dc)), dc.ratio), 6))
dm = m.by_mouse(d6j, ['ratio'])
print('main reading averaged within mouse: mean ratio', round(dm.ratio.mean(), 3), np.round(m.paired(np.ones(len(dm)), dm.ratio), 4))
var rule series mode ratio t p p_raw_pre_vs_ket
0 M M code 1-10 / 23-32 per-series 0.852 ± 0.128 1.1028 2.9595e-01 0.1045
1 M M code 1-10 / 23-32 concat 0.922 ± 0.142 0.5234 6.1209e-01 0.1208
2 M M Methods epochs 1-20 / 22-32 per-series 0.858 ± 0.125 1.0833 3.0410e-01 0.0969
3 M M Methods epochs 1-20 / 22-32 concat 1.450 ± 0.210 -2.0404 6.8601e-02 0.4216
4 M C code 1-10 / 23-32 per-series 0.832 ± 0.045 3.5334 5.4158e-03 0.0073
5 M C code 1-10 / 23-32 concat 0.724 ± 0.086 3.0531 1.2188e-02 0.0112
6 M C Methods epochs 1-20 / 22-32 per-series 0.830 ± 0.046 3.5371 5.3823e-03 0.0078
7 M C Methods epochs 1-20 / 22-32 concat 1.020 ± 0.146 -0.1296 8.9947e-01 0.2091
8 C M code 1-10 / 23-32 per-series 0.576 ± 0.101 3.9923 2.5497e-03 0.0032
9 C M code 1-10 / 23-32 concat 0.619 ± 0.103 3.5438 5.3231e-03 0.0044
10 C M Methods epochs 1-20 / 22-32 per-series 0.599 ± 0.099 3.8660 3.1294e-03 0.0048
11 C M Methods epochs 1-20 / 22-32 concat 1.333 ± 0.255 -1.2448 2.4158e-01 0.8709
12 C C code 1-10 / 23-32 per-series 0.727 ± 0.038 6.7707 4.9151e-05 0.0003
13 C C code 1-10 / 23-32 concat 0.559 ± 0.067 6.3158 8.7280e-05 0.0011
14 C C Methods epochs 1-20 / 22-32 per-series 0.747 ± 0.039 6.1678 1.0582e-04 0.0006
15 C C Methods epochs 1-20 / 22-32 concat 0.923 ± 0.130 0.5635 5.8551e-01 0.1194
main reading (M weights, |w| and mean r, code series, per series), per session:
session mouse pre ket n_cells ratio
0 sub-0_20230901 0 0.0070 0.0038 76 0.5358
1 sub-1_20230218 1 0.0192 0.0117 41 0.6106
2 sub-1_20230226 1 0.0109 0.0100 62 0.9158
3 sub-2_20230902 2 0.0086 0.0085 45 0.9924
4 sub-3_20230218 3 0.0808 0.0625 31 0.7733
5 sub-3_20230225 3 0.0403 0.0286 58 0.7091
6 sub-4_20230907 4 0.0221 0.0161 65 0.7281
7 sub-6_20230823 6 0.0139 0.0047 42 0.3411
8 sub-7_20220831 7 0.0175 0.0266 64 1.5233
9 sub-8_20220228 8 0.0301 0.0134 67 0.4441
10 sub-9_20230908 9 0.0124 0.0223 62 1.7946
ratio 0.852 ± 0.128; paired t, p = [1.1028 0.29595]
correlation during ketamine divided by preinfusion: 0.852 (expected 0.49 ± 0.07; within: False)
code reading (C weights, signed w, median |r|, concatenated): [0.559 0.067] [6.315828e+00 8.700000e-05]
main reading averaged within mouse: mean ratio 0.874 [0.7736 0.4614]
# The Methods rule averages signed r. Mean signed r per session and condition, for each reading of the Methods rule
print("analysis.py, plot_coding_dimension_early_late: the code's rule (signed w above the 90th percentile of |w|, |r|, median per session; ratio ket/pre)")
m.show_code(CODE, 1981, 1999); m.show_code(CODE, 2010, 2019); m.show_code(CODE, 2238, 2239)
SIGNED6J = {}
for var in 'MC':
cols = {}
for sname in SER6J:
for mode in ('per-series', 'concat'):
d = coupling_table(POP, var, 'M', sname, mode)
lab = f"{sname.split(' ')[0]} series, {mode}"
cols[(lab, 'pre r')] = d.pre.values; cols[(lab, 'ket r')] = d.ket.values; cols[(lab, 'ratio')] = d.ratio.values
SIGNED6J[(var, sname, mode)] = d
print(f'\nvariance {var}, Methods rule (|w| above the 90th percentile, mean signed r), per session:')
print(pd.DataFrame(cols, index=d.session).round(4).to_string())
allpre = np.concatenate([d.pre.values for d in SIGNED6J.values()])
allket = np.concatenate([d.ket.values for d in SIGNED6J.values()])
allratio = np.concatenate([d.ratio.values for d in SIGNED6J.values()])
print(f'\nover the 8 Methods-rule readings and {len(POP)} sessions: preinfusion mean signed r {allpre.min():.4f} to {allpre.max():.4f} '
f'(median {np.median(allpre):.4f}); ketamine {allket.min():.4f} to {allket.max():.4f} (median {np.median(allket):.4f}); '
f'session values with |pre r| < 0.02: {int((np.abs(allpre) < 0.02).sum())} of {len(allpre)}')
print(f'per-session ratio {allratio.min():.2f} to {allratio.max():.2f}; session-mean ratio over the 8 readings '
f'{min(d.ratio.mean() for d in SIGNED6J.values()):.3f} to {max(d.ratio.mean() for d in SIGNED6J.values()):.3f}')
cr = np.concatenate([coupling_table(POP, var, 'C', sname, mode).pre.values for var in 'MC' for sname in SER6J for mode in ('per-series', 'concat')])
print(f'for comparison, the code rule (median |r|) preinfusion: {cr.min():.4f} to {cr.max():.4f} (median {np.median(cr):.4f})')
# all 13 sessions (mouse 5 included)
for name, args in (('Methods reading', ('M', 'M', 'code 1-10 / 23-32', 'per-series')), ("code reading", ('C', 'C', 'code 1-10 / 23-32', 'concat'))):
d13 = coupling_table(Rs, *args)
bad = d13[~np.isfinite(d13.ratio)]
d13 = d13[np.isfinite(d13.ratio)]
print(f'all 13 sessions, {name}: {len(d13)} sessions with a finite value (none for {bad.session.tolist()}, top cells {bad.n_cells.tolist()}); '
f'ratio {m.mean_sem(d13.ratio)[0]:.3f} ± {m.mean_sem(d13.ratio)[1]:.3f}; t, p =', np.round(m.paired(np.ones(len(d13)), d13.ratio), 5))
analysis.py, plot_coding_dimension_early_late: the code's rule (signed w above the 90th percentile of |w|, |r|, median per session; ratio ket/pre)
1981: pre_bins = []
1982: tminus = 30
1983: tplus = 0
1984:
1985: for pid in range(10):
1986: t = data._events['puff_times'].reshape((-1,8))[pid,0]
1987: pre_bins.append(np.arange(int((t-tminus)/0.01), int((t+tplus)/0.01)))
1988:
1989: ket_bins = []
1990: for pid in range(10):
1991: t = data._events['puff_times'].reshape((-1,8))[22+pid,0]
1992: ket_bins.append(np.arange(int((t-tminus)/0.01), int((t+tplus)/0.01)))
1993:
1994: recording_cell_ids = all_cell_ix[(all_recording_names==data._recording)]
1995: top_cd = recording_cell_ids[np.argwhere(cd > np.percentile(np.abs(cd), 90)).flatten()]
1996: bottom_cd = recording_cell_ids[np.argwhere(cd < np.percentile(np.abs(cd), 10)).flatten()]
1997:
1998: corr_pre = np.abs(np.corrcoef(all_zfr[ix][:,np.hstack(pre_bins)]))
1999: corr_ket = np.abs(np.corrcoef(all_zfr[ix][:,np.hstack(ket_bins)]))
2010:
2011: top_edges_pre = []
2012: top_edges_ket = []
2013: for i_cid in top_cd:
2014: for j_cid in top_cd:
2015: if i_cid != j_cid:
2016: top_edges_pre.append(corr_pre[i_cid, j_cid])
2017: top_edges_ket.append(corr_ket[i_cid, j_cid])
2018: top_cd_edgemeans_pre.append(np.nanmedian(top_edges_pre))
2019: top_cd_edgemeans_ket.append(np.nanmedian(top_edges_ket))
2238: cd_edges = pd.DataFrame({"cd": np.hstack([2*len(top_cd_edgemeans_pre)*["Late\nneurons"]]),
2239: "corr": np.hstack([top_cd_edgemeans_pre, top_cd_edgemeans_ket]) / np.hstack([top_cd_edgemeans_pre, top_cd_edgemeans_pre]),
variance M, Methods rule (|w| above the 90th percentile, mean signed r), per session:
code series, per-series code series, concat Methods series, per-series Methods series, concat
pre r ket r ratio pre r ket r ratio pre r ket r ratio pre r ket r ratio
session
sub-0_20230901 0.0070 0.0038 0.5358 0.0080 0.0081 1.0124 0.0069 0.0036 0.5198 0.0088 0.0199 2.2680
sub-1_20230218 0.0192 0.0117 0.6106 0.0259 0.0115 0.4431 0.0156 0.0111 0.7108 0.0196 0.0132 0.6727
sub-1_20230226 0.0109 0.0100 0.9158 0.0113 0.0110 0.9792 0.0111 0.0100 0.8931 0.0114 0.0214 1.8679
sub-2_20230902 0.0086 0.0085 0.9924 0.0091 0.0158 1.7348 0.0074 0.0082 1.1148 0.0078 0.0223 2.8795
sub-3_20230218 0.0808 0.0625 0.7733 0.0833 0.0640 0.7684 0.0795 0.0636 0.7996 0.0813 0.0674 0.8286
sub-3_20230225 0.0403 0.0286 0.7091 0.0412 0.0286 0.6949 0.0398 0.0278 0.6989 0.0413 0.0309 0.7487
sub-4_20230907 0.0221 0.0161 0.7281 0.0227 0.0153 0.6750 0.0222 0.0160 0.7232 0.0230 0.0264 1.1476
sub-6_20230823 0.0139 0.0047 0.3411 0.0202 0.0072 0.3550 0.0138 0.0047 0.3388 0.0172 0.0261 1.5211
sub-7_20220831 0.0175 0.0266 1.5233 0.0205 0.0242 1.1793 0.0183 0.0253 1.3843 0.0219 0.0264 1.2067
sub-8_20220228 0.0301 0.0134 0.4441 0.0297 0.0140 0.4716 0.0309 0.0135 0.4363 0.0334 0.0255 0.7625
sub-9_20230908 0.0124 0.0223 1.7946 0.0155 0.0283 1.8274 0.0115 0.0209 1.8169 0.0148 0.0303 2.0502
variance C, Methods rule (|w| above the 90th percentile, mean signed r), per session:
code series, per-series code series, concat Methods series, per-series Methods series, concat
pre r ket r ratio pre r ket r ratio pre r ket r ratio pre r ket r ratio
session
sub-0_20230901 0.0215 0.0031 0.1461 0.0241 0.0104 0.4292 0.0211 0.0035 0.1664 0.0257 0.0452 1.7599
sub-1_20230218 0.0296 0.0075 0.2551 0.0286 0.0051 0.1781 0.0221 0.0076 0.3433 0.0266 0.0080 0.3007
sub-1_20230226 0.0107 0.0072 0.6760 0.0114 0.0125 1.0984 0.0097 0.0073 0.7476 0.0115 0.0402 3.5020
sub-2_20230902 0.0272 0.0086 0.3161 0.0279 0.0181 0.6500 0.0245 0.0091 0.3711 0.0261 0.0451 1.7299
sub-3_20230218 0.0659 0.0550 0.8343 0.0693 0.0485 0.6991 0.0647 0.0568 0.8783 0.0668 0.0614 0.9198
sub-3_20230225 0.0789 0.0450 0.5707 0.0828 0.0404 0.4882 0.0766 0.0449 0.5861 0.0807 0.0510 0.6317
sub-4_20230907 0.0209 0.0146 0.6961 0.0232 0.0139 0.5981 0.0210 0.0151 0.7186 0.0232 0.0327 1.4136
sub-6_20230823 0.0266 0.0096 0.3623 0.0439 0.0149 0.3400 0.0269 0.0092 0.3426 0.0363 0.0442 1.2162
sub-7_20220831 0.0306 0.0225 0.7350 0.0347 0.0204 0.5893 0.0319 0.0220 0.6896 0.0357 0.0280 0.7847
sub-8_20220228 0.0647 0.0225 0.3473 0.0638 0.0207 0.3251 0.0652 0.0226 0.3472 0.0684 0.0408 0.5959
sub-9_20230908 0.0152 0.0212 1.3956 0.0208 0.0293 1.4089 0.0143 0.0199 1.3927 0.0191 0.0344 1.8058
over the 8 Methods-rule readings and 11 sessions: preinfusion mean signed r 0.0069 to 0.0833 (median 0.0231); ketamine 0.0031 to 0.0674 (median 0.0206); session values with |pre r| < 0.02: 31 of 88
per-session ratio 0.15 to 3.50; session-mean ratio over the 8 readings 0.576 to 1.450
for comparison, the code rule (median |r|) preinfusion: 0.0169 to 0.0841 (median 0.0382)
all 13 sessions, Methods reading: 13 sessions with a finite value (none for [], top cells []); ratio 0.779 ± 0.124; t, p = [1.70968 0.11303]
all 13 sessions, code reading: 12 sessions with a finite value (none for ['sub-5_20230706'], top cells [7]); ratio 0.560 ± 0.061; t, p = [6.89909e+00 3.00000e-05]
Result: partially reproduced: not reproduced under the Methods' reading, reproduced under the code's reading that the criterion names as a variant. With the per-neuron-variance dimension, cells with |w| above the 90th percentile and the mean Pearson r, the ketamine/preinfusion ratio is 0.852 ± 0.128 (t = 1.10, p = 0.30), outside 0.49 ± 0.07; 9 of 11 sessions are below 1. This rule averages signed correlations that are close to zero: over its 8 readings and 11 sessions the preinfusion mean r runs from 0.0069 to 0.0833 (median 0.0231), and 31 of 88 session values are below 0.02. Ratios of such small numbers are unstable: per session they run from 0.15 to 3.50, and the session mean from 0.576 to 1.450 across the 8 readings. The code's rule (median |r|) gives preinfusion values of 0.0169 to 0.0841. With the code's pooled variance every reading over the code's series gives a reduction (the code's full reading 0.559 ± 0.067, t = 6.32, p = 8.7e-05, within tolerance; with the Methods' selection and averaging 0.576 ± 0.101, p = 0.0025). With per-neuron variance only the code's selection (signed w, median |r|) gives a reduction (0.832, p = 0.005; 0.724, p = 0.012), outside tolerance. Over the Methods' series (1 to 20, 22 to 32) the concatenated readings show no reduction. With all 13 sessions the Methods reading gives 0.779 ± 0.124 (p = 0.11) and the code reading 0.560 ± 0.061 (12 sessions with a value, p = 3e-05).
figs34-mouse-tone-persistence-reduced¶
For a neutral tone, mice do not close their eyes; fast tone activity is unchanged by ketamine, while persistent tone activity, its baseline timescale and its baseline coupling are all reduced.
Seven tone sessions from seven mice (mouse 5 excluded). Tone dimensions are built from the tone series before ketamine on the puff-modulated cells (code): fast 0 to 70 ms against the 1 s before (Methods; code 70 ms, a variant), persistent 150 to 350 ms against the 1 s before, quantified 350 to 500 ms, all 8 tones averaged. Epochs: up to 10 tone series before the split, tone series 2 to 6 after it (ketamine), 8 to 10 after it (post) (code; Table S3 10, 5, 3). Split at the code's ketamine start (30 s after the last puff of series 20: 9 preinfusion tone series) or at the file's injection time (10, as Table S3). Timescale and coupling use the 30 s before each tone series (Methods, Table S3). Coupling is run under the Methods' rule and under the code's rule (pooled variance, signed w, median |r|) over these tone windows, and over the windows before puff series 1 to 10 and 23 to 32 that the code uses for every stimulus (printed code lines), as a variant. The eye-closure part of the statement is not part of the criterion and is not tested: the behaviour files hold puff trials only.
cid = 'figs34-mouse-tone-persistence-reduced'
def tone_tables(pop, split, var, ref):
fr, pr = [], []
for R in pop:
n = R['n_sel']; W = R['tone_win'][split]
w = m.cd_weights(W['f70'].reshape(n, -1), (W['b1s'] if ref == '1s' else W['b70']).reshape(n, -1), var)
wp = R['cd_tone_persist'][split][var]
fr.append({k: m.zproj(np.tensordot(w, R['tone_ep'][split][k], axes=([0], [0])).mean(0))[100:107].mean() for k in ('pre', 'ket', 'post')})
pr.append({k: m.zproj(np.tensordot(wp, R['tone_ep'][split][k], axes=([0], [0])).mean(0))[135:150].mean() for k in ('pre', 'ket', 'post')})
return m._df(pop, fr), m._df(pop, pr)
rows = []
for split in ('code', 'inj'):
for var in 'MC':
for ref in ('1s', '70'):
df_, dp_ = tone_tables(TONE, split, var, ref)
sf, sp_ = m.three_way_tests(df_), m.three_way_tests(dp_)
r = dict(split=split, var=var, fast_ref=ref)
for k in ('pre', 'ket', 'post'):
r['fast_' + k] = '%.3f ± %.3f' % m.mean_sem(df_[k])
r.update(fast_q_pre_ket=sf['q_pre_ket'], fast_q_ket_post=sf['q_ket_post'])
for k in ('pre', 'ket', 'post'):
r['pers_' + k] = '%.3f ± %.3f' % m.mean_sem(dp_[k])
r.update(pers_t_pre_ket=sp_['t_pre_ket'], pers_q_pre_ket=sp_['q_pre_ket'], pers_t_ket_post=sp_['t_ket_post'], pers_q_ket_post=sp_['q_ket_post'])
rows.append(r)
print(pd.DataFrame(rows).to_string())
MAINS34 = ('inj', 'M', '1s')
dfs34, dps34 = tone_tables(TONE, *MAINS34)
print('\nmain reading (injection split, M, 1-s reference): fast', {k: round(v, 4) for k, v in m.three_way_tests(dfs34).items()})
print('persistent', {k: round(v, 4) for k, v in m.three_way_tests(dps34).items()})
print(dps34.to_string())
split var fast_ref fast_pre fast_ket fast_post fast_q_pre_ket fast_q_ket_post pers_pre pers_ket pers_post pers_t_pre_ket pers_q_pre_ket pers_t_ket_post pers_q_ket_post 0 code M 1s 2.778 ± 0.275 1.522 ± 0.493 1.363 ± 0.246 0.2322 0.6134 1.639 ± 0.056 0.501 ± 0.153 1.506 ± 0.310 8.5258 0.0003 -3.2201 0.0181 1 code M 70 2.268 ± 0.229 1.516 ± 0.379 1.330 ± 0.282 0.3114 0.3114 1.639 ± 0.056 0.501 ± 0.153 1.506 ± 0.310 8.5258 0.0003 -3.2201 0.0181 2 code C 1s 2.812 ± 0.305 1.691 ± 0.423 1.391 ± 0.229 0.2498 0.2524 1.522 ± 0.081 0.403 ± 0.169 1.375 ± 0.337 7.0995 0.0008 -3.5244 0.0125 3 code C 70 2.393 ± 0.227 1.650 ± 0.358 1.306 ± 0.301 0.1903 0.1689 1.522 ± 0.081 0.403 ± 0.169 1.375 ± 0.337 7.0995 0.0008 -3.5244 0.0125 4 inj M 1s 2.716 ± 0.272 1.455 ± 0.433 1.248 ± 0.246 0.1562 0.4323 1.639 ± 0.053 0.481 ± 0.162 1.509 ± 0.324 7.3013 0.0007 -2.9133 0.0269 5 inj M 70 2.283 ± 0.233 1.579 ± 0.258 1.225 ± 0.216 0.1407 0.1407 1.639 ± 0.053 0.481 ± 0.162 1.509 ± 0.324 7.3013 0.0007 -2.9133 0.0269 6 inj C 1s 2.780 ± 0.305 1.673 ± 0.346 1.301 ± 0.228 0.1011 0.1011 1.512 ± 0.076 0.356 ± 0.177 1.364 ± 0.344 6.0192 0.0019 -3.0306 0.0231 7 inj C 70 2.411 ± 0.232 1.703 ± 0.252 1.271 ± 0.208 0.1041 0.0838 1.512 ± 0.076 0.356 ± 0.177 1.364 ± 0.344 6.0192 0.0019 -3.0306 0.0231
main reading (injection split, M, 1-s reference): fast {'t_pre_ket': 2.1219, 'p_pre_ket': 0.0781, 'q_pre_ket': np.float64(0.1562), 't_ket_post': 0.8415, 'p_ket_post': 0.4323, 'q_ket_post': np.float64(0.4323)}
persistent {'t_pre_ket': 7.3013, 'p_pre_ket': 0.0003, 'q_pre_ket': np.float64(0.0007), 't_ket_post': -2.9133, 'p_ket_post': 0.0269, 'q_ket_post': np.float64(0.0269)}
session mouse pre ket post
0 sub-0_20230901 0 1.8638 0.4595 2.0575
1 sub-1_20230226 1 1.3616 0.1121 1.2597
2 sub-2_20230902 2 1.6924 1.3686 0.9996
3 sub-3_20230225 3 1.6785 0.5775 2.8295
4 sub-4_20230907 4 1.5878 0.2697 1.7793
5 sub-6_20230823 6 1.6119 -0.0477 -0.1157
6 sub-9_20230908 9 1.6755 0.6247 1.7519
rows = []
with warnings.catch_warnings():
warnings.simplefilter('ignore')
for split in ('code', 'inj'):
for var in 'MC':
for lag in (0, 1):
for mode in ('per-series', 'concat'):
d, nf = m.timescale_table(TONE, var, lambda R: R['tone_sets'][split]['pre'],
lambda R: R['tone_sets'][split]['ket'], lag, mode,
key=lambda R: R['tone_base_proj'][split])
r = dict(split=split, var=var, lag=lag, mode=mode, failed=nf, pre='%.3f ± %.3f' % m.mean_sem(d.pre),
ket='%.3f ± %.3f' % m.mean_sem(d.ket))
r['t'], r['p'] = m.paired(d.pre, d.ket); rows.append(r)
print(pd.DataFrame(rows).to_string())
dts34, _ = m.timescale_table(TONE, 'M', lambda R: R['tone_sets']['inj']['pre'], lambda R: R['tone_sets']['inj']['ket'], 0,
'per-series', key=lambda R: R['tone_base_proj']['inj'])
print('\nmain reading timescales:'); print(dts34.to_string()); print('t, p =', np.round(m.paired(dts34.pre, dts34.ket), 4))
rows = []
for split in ('code', 'inj'):
for var in 'MC':
for rule in 'MC':
for mode in ('per-series', 'concat', 'code puff windows'):
vals = []
for R in TONE:
ts = R['tone_sets'][split]
if mode == 'per-series':
v = R['tone_coupling_series'][split][(var, rule)]; vals.append((np.nanmean(v[ts['pre']]), np.nanmean(v[ts['ket']])))
elif mode == 'concat':
c = R['tone_coupling'][split]; vals.append((c[(var, rule, 'pre')], c[(var, rule, 'ket')]))
else:
c = R['tone_coupling_puffwin'][split]; vals.append((c[(var, rule, 'pre_C')], c[(var, rule, 'ket_C')]))
vals = np.array(vals); ratio = vals[:, 1] / vals[:, 0]
tt_, pp_ = m.paired(np.ones(len(ratio)), ratio)
rows.append(dict(split=split, var=var, rule=rule, windows=mode, ratio='%.3f ± %.3f' % m.mean_sem(ratio), t=tt_, p=pp_))
COUP34 = pd.DataFrame(rows)
print(COUP34.to_string())
ex = expected(cid)
print('\nmain reading against expected values:')
vals = [m.mean_sem(dps34.pre)[0], m.mean_sem(dps34.ket)[0], m.mean_sem(dts34.pre)[0], m.mean_sem(dts34.ket)[0]]
main_c = COUP34[(COUP34['split'] == 'inj') & (COUP34['var'] == 'M') & (COUP34['rule'] == 'M') & (COUP34['windows'] == 'per-series')].iloc[0]
vals.append(float(main_c.ratio.split(' ± ')[0]))
for (qname, v, tol), x in zip(ex, vals):
print(f'{qname}: {x:.3f} (expected {v} ± {tol}; within: {within(x, v, tol)})')
print('main coupling test:', main_c.to_dict())
split var lag mode failed pre ket t p
0 code M 0 per-series 0 0.291 ± 0.035 0.141 ± 0.035 3.7715 0.0093
1 code M 0 concat 0 0.236 ± 0.027 0.087 ± 0.011 5.1847 0.0020
2 code M 1 per-series 0 0.310 ± 0.040 0.165 ± 0.049 2.6385 0.0386
3 code M 1 concat 0 0.243 ± 0.028 0.089 ± 0.014 5.2932 0.0018
4 code C 0 per-series 0 0.312 ± 0.040 0.146 ± 0.034 3.6121 0.0112
5 code C 0 concat 0 0.233 ± 0.023 0.088 ± 0.009 5.9904 0.0010
6 code C 1 per-series 0 0.350 ± 0.061 0.161 ± 0.041 2.6641 0.0373
7 code C 1 concat 0 0.239 ± 0.024 0.089 ± 0.011 6.2436 0.0008
8 inj M 0 per-series 0 0.304 ± 0.043 0.115 ± 0.027 5.4250 0.0016
9 inj M 0 concat 0 0.245 ± 0.032 0.080 ± 0.010 5.0144 0.0024
10 inj M 1 per-series 0 0.321 ± 0.046 0.144 ± 0.048 3.7930 0.0090
11 inj M 1 concat 0 0.252 ± 0.033 0.082 ± 0.014 5.2863 0.0019
12 inj C 0 per-series 0 0.319 ± 0.043 0.122 ± 0.031 4.8732 0.0028
13 inj C 0 concat 0 0.239 ± 0.028 0.082 ± 0.007 5.7036 0.0013
14 inj C 1 per-series 0 0.351 ± 0.056 0.145 ± 0.046 3.2364 0.0178
15 inj C 1 concat 0 0.246 ± 0.029 0.083 ± 0.010 6.1182 0.0009
main reading timescales:
session mouse pre ket
0 sub-0_20230901 0 0.3765 0.2277
1 sub-1_20230226 1 0.2862 0.0741
2 sub-2_20230902 2 0.1191 0.0864
3 sub-3_20230225 3 0.3317 0.0351
4 sub-4_20230907 4 0.1910 0.0654
5 sub-6_20230823 6 0.4940 0.2248
6 sub-9_20230908 9 0.3289 0.0942
t, p = [5.425e+00 1.600e-03]
split var rule windows ratio t p
0 code M M per-series 1.252 ± 0.457 -0.5100 0.6282
1 code M M concat 1.565 ± 0.462 -1.1321 0.3008
2 code M M code puff windows 1.220 ± 0.443 -0.4591 0.6623
3 code M C per-series 0.872 ± 0.111 1.0696 0.3259
4 code M C concat 1.066 ± 0.211 -0.2886 0.7826
5 code M C code puff windows 0.647 ± 0.124 2.6401 0.0385
6 code C M per-series 0.917 ± 0.329 0.2328 0.8237
7 code C M concat 1.173 ± 0.318 -0.5026 0.6332
8 code C M code puff windows 0.751 ± 0.234 0.9852 0.3626
9 code C C per-series 0.759 ± 0.090 2.4919 0.0470
10 code C C concat 1.019 ± 0.241 -0.0730 0.9442
11 code C C code puff windows 0.540 ± 0.106 4.0385 0.0068
12 inj M M per-series 1.238 ± 0.451 -0.4894 0.6419
13 inj M M concat 1.296 ± 0.508 -0.5389 0.6094
14 inj M M code puff windows 1.180 ± 0.410 -0.4063 0.6986
15 inj M C per-series 0.834 ± 0.097 1.5837 0.1643
16 inj M C concat 0.809 ± 0.179 0.9890 0.3608
17 inj M C code puff windows 0.671 ± 0.129 2.3636 0.0560
18 inj C M per-series 0.906 ± 0.344 0.2542 0.8078
19 inj C M concat 0.918 ± 0.381 0.1996 0.8484
20 inj C M code puff windows 0.732 ± 0.224 1.1099 0.3095
21 inj C C per-series 0.741 ± 0.092 2.5978 0.0408
22 inj C C concat 0.664 ± 0.163 1.9075 0.1051
23 inj C C code puff windows 0.535 ± 0.104 4.1387 0.0061
main reading against expected values:
persistent tone activity preinfusion: 1.639 (expected 1.531 ± 0.081; within: False)
persistent tone activity ketamine: 0.481 (expected 0.392 ± 0.177; within: True)
tone timescale preinfusion: 0.304 (expected 0.268 ± 0.037; within: True)
tone timescale ketamine: 0.115 (expected 0.089 ± 0.012; within: False)
tone coupling ketamine divided by preinfusion: 1.238 (expected 0.46 ± 0.08; within: False)
main coupling test: {'split': 'inj', 'var': 'M', 'rule': 'M', 'windows': 'per-series', 'ratio': '1.238 ± 0.451', 't': -0.48942156410616156, 'p': 0.6419200387661228}
# Coupling over the 30 s before each TONE series (Methods, Table S3), under the Methods rule and under the code's rule
# (pooled variance C, signed w above the 90th percentile of |w|, median |r|), next to the code's windows before puff series.
print("analysis.py, plot_coding_dimension_early_late: the tone branch, and the coupling windows, taken before puff series 1 to 10 and 23 to 32 for every stimulus")
m.show_code(CODE, 1954, 1955); m.show_code(CODE, 1981, 1992)
def tone_coupling(pop, split, var, rule, windows):
rows = []
for R in pop:
ts = R['tone_sets'][split]
if windows == 'tone series, per-series':
v = R['tone_coupling_series'][split][(var, rule)]; rows.append(dict(pre=np.nanmean(v[ts['pre']]), ket=np.nanmean(v[ts['ket']])))
elif windows == 'tone series, concat':
c = R['tone_coupling'][split]; rows.append(dict(pre=c[(var, rule, 'pre')], ket=c[(var, rule, 'ket')]))
else:
c = R['tone_coupling_puffwin'][split]; rows.append(dict(pre=c[(var, rule, 'pre_C')], ket=c[(var, rule, 'ket_C')]))
d = m._df(pop, rows); d['ratio'] = d.ket / d.pre
return d
def tone_coupling_rows(pop):
rows = []
for split in ('inj', 'code'):
for var, rule, rname in (('M', 'M', 'Methods: M, |w|, mean r'), ('C', 'C', 'code: C, signed w, median |r|')):
for windows in ('tone series, per-series', 'tone series, concat', 'puff series 1-10 / 23-32, concat'):
d = tone_coupling(pop, split, var, rule, windows)
mu, se = m.mean_sem(d.ratio); t1, p1 = m.paired(np.ones(len(d)), d.ratio); t2, p2 = m.paired(d.pre, d.ket)
rows.append(dict(split=split, reading=rname, windows=windows, ratio=f'{mu:.3f} ± {se:.3f}', within_tol=within(mu, 0.46, 0.08),
t_1_vs_ratio=t1, p_1_vs_ratio=p1, t_pre_vs_ket=t2, p_pre_vs_ket=p2, sessions_below_1=int((d.ratio < 1).sum())))
return pd.DataFrame(rows)
print(tone_coupling_rows(TONE).to_string())
for split, var, rule, windows in (('inj', 'M', 'M', 'tone series, per-series'), ('code', 'C', 'C', 'tone series, concat'), ('inj', 'C', 'C', 'tone series, concat')):
d = tone_coupling(TONE, split, var, rule, windows)
print(f'\nper session, split {split}, variance {var}, rule {rule}, {windows} (pre and ket are {"mean signed r" if rule == "M" else "median |r|"}):')
print(d.round(4).to_string())
# all 9 tone sessions (mouse 5 included), main readings
TONE13 = [R for R in Rs if R['tone'] is not None]
_, dp9 = tone_tables(TONE13, *MAINS34)
print('\n9 tone sessions (mouse 5 included): persistent', {k: round(m.mean_sem(dp9[k])[0], 3) for k in ('pre', 'ket', 'post')},
{k: round(v, 4) for k, v in m.three_way_tests(dp9).items() if k.startswith('q') or k.startswith('t')})
with warnings.catch_warnings():
warnings.simplefilter('ignore')
dt9, _ = m.timescale_table(TONE13, 'M', lambda R: R['tone_sets']['inj']['pre'], lambda R: R['tone_sets']['inj']['ket'], 0,
'per-series', key=lambda R: R['tone_base_proj']['inj'])
print('9 tone sessions: timescale', {k: round(m.mean_sem(dt9[k])[0], 3) for k in ('pre', 'ket')}, 't, p =', np.round(m.paired(dt9.pre, dt9.ket), 4))
print(tone_coupling_rows(TONE13)[['split', 'reading', 'windows', 'ratio', 't_1_vs_ratio', 'p_1_vs_ratio']].to_string())
analysis.py, plot_coding_dimension_early_late: the tone branch, and the coupling windows, taken before puff series 1 to 10 and 23 to 32 for every stimulus
1954: if stimulus=='tone':
1955: late_start, late_stop = (115, 135)
1981: pre_bins = []
1982: tminus = 30
1983: tplus = 0
1984:
1985: for pid in range(10):
1986: t = data._events['puff_times'].reshape((-1,8))[pid,0]
1987: pre_bins.append(np.arange(int((t-tminus)/0.01), int((t+tplus)/0.01)))
1988:
1989: ket_bins = []
1990: for pid in range(10):
1991: t = data._events['puff_times'].reshape((-1,8))[22+pid,0]
1992: ket_bins.append(np.arange(int((t-tminus)/0.01), int((t+tplus)/0.01)))
split reading windows ratio within_tol t_1_vs_ratio p_1_vs_ratio t_pre_vs_ket p_pre_vs_ket sessions_below_1
0 inj Methods: M, |w|, mean r tone series, per-series 1.238 ± 0.451 False -0.4894 0.6419 0.1439 0.8903 4
1 inj Methods: M, |w|, mean r tone series, concat 1.296 ± 0.508 False -0.5389 0.6094 0.2924 0.7798 4
2 inj Methods: M, |w|, mean r puff series 1-10 / 23-32, concat 1.180 ± 0.410 False -0.4063 0.6986 0.4475 0.6702 4
3 inj code: C, signed w, median |r| tone series, per-series 0.741 ± 0.092 False 2.5978 0.0408 2.5615 0.0428 6
4 inj code: C, signed w, median |r| tone series, concat 0.664 ± 0.163 False 1.9075 0.1051 2.2789 0.0629 5
5 inj code: C, signed w, median |r| puff series 1-10 / 23-32, concat 0.535 ± 0.104 True 4.1387 0.0061 3.1241 0.0205 7
6 code Methods: M, |w|, mean r tone series, per-series 1.252 ± 0.457 False -0.5100 0.6282 -0.0283 0.9783 4
7 code Methods: M, |w|, mean r tone series, concat 1.565 ± 0.462 False -1.1321 0.3008 -0.6397 0.5460 3
8 code Methods: M, |w|, mean r puff series 1-10 / 23-32, concat 1.220 ± 0.443 False -0.4591 0.6623 0.2962 0.7771 4
9 code code: C, signed w, median |r| tone series, per-series 0.759 ± 0.090 False 2.4919 0.0470 2.4377 0.0506 6
10 code code: C, signed w, median |r| tone series, concat 1.019 ± 0.241 False -0.0730 0.9442 0.8339 0.4363 5
11 code code: C, signed w, median |r| puff series 1-10 / 23-32, concat 0.540 ± 0.106 True 4.0385 0.0068 3.1489 0.0198 6
per session, split inj, variance M, rule M, tone series, per-series (pre and ket are mean signed r):
session mouse pre ket ratio
0 sub-0_20230901 0 0.0263 0.0087 0.3319
1 sub-1_20230226 1 0.0280 0.0160 0.5727
2 sub-2_20230902 2 0.0136 0.0508 3.7260
3 sub-3_20230225 3 0.0304 0.0068 0.2236
4 sub-4_20230907 4 0.0298 0.0382 1.2838
5 sub-6_20230823 6 0.0516 0.0211 0.4088
6 sub-9_20230908 9 0.0250 0.0530 2.1208
per session, split code, variance C, rule C, tone series, concat (pre and ket are median |r|):
session mouse pre ket ratio
0 sub-0_20230901 0 0.0460 0.0321 0.6992
1 sub-1_20230226 1 0.0330 0.0282 0.8535
2 sub-2_20230902 2 0.0245 0.0463 1.8881
3 sub-3_20230225 3 0.0668 0.0269 0.4032
4 sub-4_20230907 4 0.0381 0.0253 0.6648
5 sub-6_20230823 6 0.0575 0.0295 0.5129
6 sub-9_20230908 9 0.0222 0.0469 2.1111
per session, split inj, variance C, rule C, tone series, concat (pre and ket are median |r|):
session mouse pre ket ratio
0 sub-0_20230901 0 0.0467 0.0194 0.4155
1 sub-1_20230226 1 0.0324 0.0183 0.5629
2 sub-2_20230902 2 0.0240 0.0323 1.3486
3 sub-3_20230225 3 0.0677 0.0167 0.2464
4 sub-4_20230907 4 0.0361 0.0194 0.5368
5 sub-6_20230823 6 0.0583 0.0142 0.2436
6 sub-9_20230908 9 0.0227 0.0294 1.2934
9 tone sessions (mouse 5 included): persistent {'pre': np.float64(1.642), 'ket': np.float64(0.377), 'post': np.float64(1.316)} {'t_pre_ket': 5.8321, 'q_pre_ket': np.float64(0.0008), 't_ket_post': -3.1682, 'q_ket_post': np.float64(0.0132)}
9 tone sessions: timescale {'pre': np.float64(0.263), 'ket': np.float64(0.134)} t, p = [2.5792 0.0327]
split reading windows ratio t_1_vs_ratio p_1_vs_ratio
0 inj Methods: M, |w|, mean r tone series, per-series 1.053 ± 0.369 -0.1350 0.8959
1 inj Methods: M, |w|, mean r tone series, concat 1.058 ± 0.424 -0.1298 0.8999
2 inj Methods: M, |w|, mean r puff series 1-10 / 23-32, concat 0.967 ± 0.347 0.0890 0.9313
3 inj code: C, signed w, median |r| tone series, per-series 0.701 ± 0.088 3.2105 0.0124
4 inj code: C, signed w, median |r| tone series, concat 0.613 ± 0.142 2.5813 0.0325
5 inj code: C, signed w, median |r| puff series 1-10 / 23-32, concat 0.491 ± 0.092 5.1947 0.0008
6 code Methods: M, |w|, mean r tone series, per-series 1.060 ± 0.381 -0.1477 0.8862
7 code Methods: M, |w|, mean r tone series, concat 1.302 ± 0.400 -0.7121 0.4966
8 code Methods: M, |w|, mean r puff series 1-10 / 23-32, concat 1.006 ± 0.374 -0.0156 0.9880
9 code code: C, signed w, median |r| tone series, per-series 0.721 ± 0.081 3.2287 0.0121
10 code code: C, signed w, median |r| tone series, concat 0.889 ± 0.210 0.4996 0.6308
11 code code: C, signed w, median |r| puff series 1-10 / 23-32, concat 0.493 ± 0.094 5.0996 0.0009
Result: partially reproduced. Fast tone activity is not significantly changed (main reading corrected p = 0.16 pre against ketamine and 0.43 post against ketamine; every reading above 0.08). Persistent tone activity falls from 1.639 ± 0.053 to 0.481 ± 0.162 (t = 7.30, corrected p = 0.0007) and recovers (t = -2.91, corrected p = 0.027); the preinfusion mean is above 1.531 ± 0.081, while the code's pooled variance gives 1.512 to 1.522, within it. The tone timescale falls from 0.304 ± 0.043 to 0.115 ± 0.027 s (t = 5.43, p = 0.0016); the preinfusion value is within 0.268 ± 0.037 and the ketamine value outside 0.089 ± 0.012. The coupling part fails over the tone windows of the Methods and Table S3. Under the Methods' rule the ratio is 1.238 ± 0.451 (p = 0.64), with preinfusion signed r of 0.0136 to 0.0516 per session. Under the code's rule over the same windows the concatenated ratio is 0.664 ± 0.163 (p = 0.105) with the injection split and 1.019 ± 0.241 (p = 0.94) with the code's split; averaged per series it is 0.741 ± 0.092 (p = 0.041) and 0.759 ± 0.090 (p = 0.047), outside 0.46 ± 0.08. The ratio comes within tolerance only over the windows before puff series (0.535 ± 0.104, p = 0.006; 0.540 ± 0.106, p = 0.0068 with the code's split), which are not the tone windows. With the 9 tone sessions (mouse 5 included) persistent activity falls (corrected p = 0.0008), the timescale falls (p = 0.033), and the code's rule over tone windows gives ratios of 0.613 to 0.889.
fig4g-mouse-slow-clusters-altered¶
Ketamine leaves neuron clusters with rapid preinfusion puff responses little changed and substantially alters clusters with slower, more persistent responses.
All 13 sessions. Each puff-modulated cell's preinfusion per-puff trace (mean over the 20 series and 8 puffs, baseline subtracted) enters Scanpy (PCA arpack, 30 neighbours on 20 PCs, UMAP, Leiden with leidenalg at resolution 1). Per cluster, cells' traces before and during ketamine (trials as in fig6d), z-scored within cluster and condition, mean and 95% interval over cells; the measure is the fraction of 0 to 2 s after the puff in which the two intervals do not overlap; clusters ordered by the time their preinfusion average decays to 25% of its peak (code find_decay).
cid = 'fig4g-mouse-slow-clusters-altered'
t0 = time.time()
with warnings.catch_warnings():
warnings.simplefilter('ignore')
LAB, SESSIX = m.cluster_cells(Rs)
print(f'clustering {len(LAB)} cells: {len(np.unique(LAB))} clusters in {time.time() - t0:.0f} s; sizes {np.bincount(LAB).tolist()}')
mice_of = np.array([Rs[i]['subject'] for i in SESSIX])
print('mice contributing to each cluster:', [len(np.unique(mice_of[LAB == c])) for c in np.unique(LAB)])
def cluster_measure(lab, ket):
tr = m.cluster_condition_traces(Rs, lab, ket)
rows = []
for c, d in tr.items():
mp, cp, n = d['pre']; mk, ck, _ = d['ket']
rows.append(dict(cluster=c, cells=n, decay_s=m.find_decay(mp), nonoverlap_0_2s=np.mean((np.abs(mp - mk) > (cp + ck))[100:300])))
d = pd.DataFrame(rows).sort_values('decay_s').reset_index(drop=True)
rho, p = scipy.stats.spearmanr(d.decay_s, d.nonoverlap_0_2s)
return d, rho, p, d.nonoverlap_0_2s.iloc[-3:].mean(), d.nonoverlap_0_2s.iloc[:3].mean(), tr
for ket in ('ket_C', 'ket_M'):
d, rho, p, slow, fast, tr = cluster_measure(LAB, ket)
print(f'\nketamine trials {ket}:'); print(d.to_string())
print(f'Spearman rho {rho:.3f} (p {p:.3f}); slowest 3 {slow:.3f}, fastest 3 {fast:.3f}')
if ket == 'ket_C':
D4G, TR4G = d, tr
with warnings.catch_warnings():
warnings.simplefilter('ignore')
for res in (1.2, 1.5):
lab2, _ = m.cluster_cells(Rs, resolution=res)
d, rho, p, slow, fast, _ = cluster_measure(lab2, 'ket_C')
print(f'resolution {res}: {len(np.unique(lab2))} clusters, Spearman rho {rho:.3f} (p {p:.3f}); slowest 3 {slow:.3f}, fastest 3 {fast:.3f}')
clustering 6553 cells: 11 clusters in 30 s; sizes [1054, 800, 797, 717, 573, 551, 520, 461, 382, 359, 339]
mice contributing to each cluster: [10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10]
ketamine trials ket_C:
cluster cells decay_s nonoverlap_0_2s
0 8 382 0.09 0.255
1 4 573 0.25 0.335
2 2 797 0.30 0.515
3 5 551 0.35 0.280
4 7 461 0.37 0.380
5 3 717 0.47 0.520
6 6 520 0.62 0.375
7 9 359 0.79 0.425
8 1 800 1.02 0.685
9 10 339 1.20 0.540
10 0 1054 1.24 0.295
Spearman rho 0.473 (p 0.142); slowest 3 0.507, fastest 3 0.368
ketamine trials ket_M:
cluster cells decay_s nonoverlap_0_2s
0 8 382 0.09 0.255
1 4 573 0.25 0.330
2 2 797 0.30 0.480
3 5 551 0.35 0.335
4 7 461 0.37 0.375
5 3 717 0.47 0.545
6 6 520 0.62 0.365
7 9 359 0.79 0.370
8 1 800 1.02 0.665
9 10 339 1.20 0.525
10 0 1054 1.24 0.235
Spearman rho 0.273 (p 0.417); slowest 3 0.475, fastest 3 0.355
resolution 1.2: 13 clusters, Spearman rho 0.537 (p 0.058); slowest 3 0.612, fastest 3 0.187
resolution 1.5: 17 clusters, Spearman rho 0.331 (p 0.194); slowest 3 0.398, fastest 3 0.267
fig, axs = plt.subplots(2, 6, figsize=(14, 4.5), sharex=True)
for ax, (_, row) in zip(axs.ravel(), D4G.iterrows()):
d = TR4G[row.cluster]
for k, col in (('pre', 'C0'), ('ket', 'C3')):
mu, ci, _ = d[k]; ax.plot(t, mu, color=col); ax.fill_between(t, mu - ci, mu + ci, color=col, alpha=.3)
ax.set_title(f"cl {row.cluster}: decay {row.decay_s:.2f} s", fontsize=8); ax.set_xlim(-0.5, 2)
plt.tight_layout(); plt.show()
Result: reproduced. The 6553 modulated cells of 13 sessions form 11 clusters (the paper reports 12), each with cells from all 10 mice. With ketamine trials 23 to 32, the fraction of 0 to 2 s in which the two intervals separate rises with the preinfusion decay time (Spearman rho = 0.47, p = 0.14), and the three slowest clusters average 0.507 against 0.368 for the three fastest. With trials 22 to 32 rho = 0.27 (0.475 against 0.355); at resolutions giving 13 and 17 clusters rho = 0.54 and 0.33, slowest above fastest in both. The correlation is not significant in any version, and the three fastest clusters still separate over 0.255 to 0.515 of the window.
fig4hi-mouse-rise-unchanged-decay-shortened¶
Ketamine leaves the rise times of puff-evoked regional activity unchanged but shortens regional decay times; decay times run from subcortical to frontal regions.
All 13 sessions, clusters of fig4g. Per region (layer suffixes merged, code reduce_region; at least 30 modulated cells; 'NA', 'STR', 'MB', 'HY' left out as in the code's plots), cluster averages of the first two puffs of each series (40 preinfusion replicates, 20 ketamine), d' against the 1-s baseline, rise to and decay to 25% of the normalised peak, peaks after 750 ms dropped, cluster-weighted; 1000 bootstrap resamples of the puff replicates (seed 0). Region filter on peak d' 1.3 (Methods, code) or 2 (Table S3). Two readings of "its bootstrap 95% interval": the interval of the median over regions across bootstrap samples, or the 2.5 to 97.5 percentiles of the code's pooled distribution of differences over regions and bootstrap samples.
cid = 'fig4hi-mouse-rise-unchanged-decay-shortened'
rows = []
t0 = time.time()
with warnings.catch_warnings():
warnings.simplefilter('ignore')
for ket_name, ket_series in (('code 23-32', np.arange(22, 32)), ('Methods 22-32', np.arange(21, 32))):
for thr in (1.3, 2.0):
res = m.region_rise_decay(Rs, LAB, ket_series, dprime_min=thr, n_boot=1000, seed=0)
regs = list(res)
D = np.stack([res[r]['ket'][1] - res[r]['pre'][1] for r in regs])
RR = np.stack([res[r]['ket'][0] - res[r]['pre'][0] for r in regs])
md_, mr_ = np.nanmedian(D, 0), np.nanmedian(RR, 0)
rows.append(dict(ket=ket_name, dprime_min=thr, regions=len(regs), decay_median=np.nanmedian(D),
decay_ci_of_median=np.round(np.percentile(md_, [2.5, 97.5]), 4).tolist(),
decay_pooled_95=np.round(np.nanpercentile(D, [2.5, 97.5]), 4).tolist(),
rise_median=np.nanmedian(RR), rise_ci_of_median=np.round(np.percentile(mr_, [2.5, 97.5]), 4).tolist(),
rise_pooled_95=np.round(np.nanpercentile(RR, [2.5, 97.5]), 4).tolist()))
if ket_name == 'code 23-32' and thr == 1.3:
RES4H = res
print(f'{time.time() - t0:.0f} s')
TAB4H = pd.DataFrame(rows); print(TAB4H.to_string())
reg_rows = [dict(region=r, cells=v['n'], peak_dprime=v['dmax'], pre_rise=np.nanmedian(v['pre'][0]), ket_rise=np.nanmedian(v['ket'][0]),
pre_decay=np.nanmedian(v['pre'][1]), ket_decay=np.nanmedian(v['ket'][1])) for r, v in RES4H.items()]
REG4H = pd.DataFrame(reg_rows).sort_values('pre_decay')
print(REG4H.to_string())
print('regions with shorter decay on ketamine:', int((REG4H.ket_decay < REG4H.pre_decay).sum()), 'of', len(REG4H),
'; with earlier rise:', int((REG4H.ket_rise < REG4H.pre_rise).sum()), ', later rise:', int((REG4H.ket_rise > REG4H.pre_rise).sum()))
116 s
ket dprime_min regions decay_median decay_ci_of_median decay_pooled_95 rise_median rise_ci_of_median rise_pooled_95
0 code 23-32 1.3 36 -0.1885 [-0.2312, -0.1495] [-0.4859, 0.0884] -0.0094 [-0.0184, -0.0017] [-0.0827, 0.0675]
1 code 23-32 2.0 30 -0.1857 [-0.2283, -0.1485] [-0.4755, 0.0289] -0.0094 [-0.0177, -0.0009] [-0.0736, 0.047]
2 Methods 22-32 1.3 36 -0.1828 [-0.2235, -0.1432] [-0.465, 0.1078] -0.0089 [-0.0177, -0.0009] [-0.0816, 0.0617]
3 Methods 22-32 2.0 30 -0.1804 [-0.2236, -0.1412] [-0.4555, 0.0397] -0.0091 [-0.0174, -0.0008] [-0.0749, 0.0457]
region cells peak_dprime pre_rise ket_rise pre_decay ket_decay
30 RT 36 3.6780 0.0258 0.0242 0.2252 0.1400
24 PAG 137 3.1705 0.0101 0.0117 0.2420 0.1340
14 LSv 44 1.3059 0.0347 0.0683 0.2508 0.3123
19 NPC 84 3.5220 0.0116 0.0176 0.2537 0.1129
6 CA3 35 2.7657 0.0309 0.0226 0.2779 0.1874
17 MOs 73 2.0278 0.0470 0.0467 0.3109 0.2582
10 FF 51 3.1287 0.0506 0.0473 0.3124 0.2323
31 SCig 42 2.2804 0.0927 0.0449 0.3159 0.1510
18 MRN 79 2.3465 0.0160 0.0076 0.3410 0.1978
28 PoT 68 2.1875 0.0341 0.0256 0.3417 0.2239
11 ILA 61 2.4526 0.0420 0.0241 0.3553 0.1406
33 TH 53 3.5176 0.0656 0.0361 0.3691 0.2376
16 MOp 91 1.7602 0.0572 0.1086 0.3785 0.3901
35 ZI 169 5.1983 0.0267 0.0179 0.4044 0.1649
12 LP 81 2.4257 0.0756 0.0404 0.4536 0.2406
32 SUB 46 2.7501 0.0382 0.0489 0.4589 0.2719
8 DG-mo 56 2.8655 0.0441 0.0340 0.4632 0.2038
27 PO 275 3.5447 0.0380 0.0281 0.4842 0.3704
22 ORBm 44 2.0584 0.0895 0.0377 0.4869 0.1755
5 CA1 97 3.6950 0.0318 0.0410 0.4924 0.2276
1 ACAv 153 3.4339 0.0630 0.0540 0.4947 0.2460
29 RSPv 292 3.4583 0.0711 0.0785 0.5014 0.3779
15 MD 185 3.2888 0.0873 0.0488 0.5251 0.2911
7 CP 237 2.9285 0.0516 0.0498 0.5261 0.3985
3 AIv 131 1.9149 0.0677 0.0559 0.5371 0.2757
34 VPM 131 2.9765 0.0534 0.0408 0.5498 0.2306
0 ACAd 147 3.8787 0.0504 0.0429 0.5522 0.2751
2 ACB 189 3.0656 0.0526 0.0420 0.5630 0.3450
9 EPd 95 1.3038 0.0729 0.0465 0.5772 0.2065
21 ORBl 383 4.3024 0.0966 0.0898 0.5839 0.4037
26 PL 186 3.6307 0.0426 0.0289 0.5865 0.2127
23 ORBvl 125 2.0040 0.1005 0.0656 0.6467 0.3319
4 AON 62 1.4667 0.1199 0.0556 0.6539 0.4219
20 OLF 176 2.3368 0.1010 0.0714 0.6840 0.3870
25 PIR 153 1.7013 0.1444 0.0961 0.7204 0.3825
13 LSr 398 3.9839 0.0341 0.0266 0.7278 0.2395
regions with shorter decay on ketamine: 34 of 36 ; with earlier rise: 29 , later rise: 7
# Each part of the criterion under each reading of 'its bootstrap 95% interval'
chk = TAB4H[['ket', 'dprime_min', 'decay_median', 'rise_median']].copy()
for name, dcol, rcol in (('median', 'decay_ci_of_median', 'rise_ci_of_median'), ('pooled', 'decay_pooled_95', 'rise_pooled_95')):
chk[f'{name}: decay CI'] = TAB4H[dcol]; chk[f'{name}: rise CI'] = TAB4H[rcol]
chk[f'{name}: decay below 0'] = (TAB4H.decay_median < 0) & np.array([ci[1] < 0 for ci in TAB4H[dcol]])
chk[f'{name}: rise includes 0'] = np.array([ci[0] <= 0 <= ci[1] for ci in TAB4H[rcol]])
chk[f'{name}: both'] = chk[f'{name}: decay below 0'] & chk[f'{name}: rise includes 0']
print(chk.to_string())
print('readings (ketamine trials x d\' threshold x interval) satisfying both parts:',
int(chk['median: both'].sum() + chk['pooled: both'].sum()), 'of', 2 * len(chk))
ket dprime_min decay_median rise_median median: decay CI median: rise CI median: decay below 0 median: rise includes 0 median: both pooled: decay CI pooled: rise CI pooled: decay below 0 pooled: rise includes 0 pooled: both 0 code 23-32 1.3 -0.1885 -0.0094 [-0.2312, -0.1495] [-0.0184, -0.0017] True False False [-0.4859, 0.0884] [-0.0827, 0.0675] False True False 1 code 23-32 2.0 -0.1857 -0.0094 [-0.2283, -0.1485] [-0.0177, -0.0009] True False False [-0.4755, 0.0289] [-0.0736, 0.047] False True False 2 Methods 22-32 1.3 -0.1828 -0.0089 [-0.2235, -0.1432] [-0.0177, -0.0009] True False False [-0.465, 0.1078] [-0.0816, 0.0617] False True False 3 Methods 22-32 2.0 -0.1804 -0.0091 [-0.2236, -0.1412] [-0.0174, -0.0008] True False False [-0.4555, 0.0397] [-0.0749, 0.0457] False True False readings (ketamine trials x d' threshold x interval) satisfying both parts: 0 of 8
Result: partially reproduced. Over 36 regions (peak d' of at least 1.3) decay times shorten on ketamine: median change -0.189 s, interval of the median -0.231 to -0.150 s, 34 of 36 regions shorter; with d' of at least 2 (30 regions) -0.186 s. Rise times change by a median of -0.009 s, a twentieth of the decay change, and 29 of 36 regions rise earlier. No single reading of the interval satisfies both parts of the criterion: 0 of the 8 combinations of ketamine trials, d' threshold and interval reading do. Under the interval of the median the decay part holds in all 4 combinations and the rise part in none (rise interval -0.018 to -0.002 s at d' 1.3). Under the code's pooled distribution over regions and bootstrap samples the rise part holds in all 4 (-0.083 to 0.068 s at d' 1.3) and the decay part in none (-0.486 to 0.088 s). The subcortical-to-frontal ordering of decay times is not part of the criterion and is not judged.
fig5g-mouse-first-order-model¶
A two-phase first-order model in which only the persistence time constant differs between preinfusion and ketamine fits the emotion-like activity closely, and varying only that constant captures the dose dependence of affective eye closure for ketamine and PCP; alternative models varying the input scale or the broadcast constant fit worse.
Neural part: the model of the Methods (code model_system) fitted with curve_fit (parameters in [0, 100], condition multipliers in [0, 1]) to the across-session mean emotion-like signal of fig5e, preinfusion and ketamine concatenated (-0.25 to 23 s and 25 to 70 s after the first puff); variance explained 1 - MSE(model)/MSE(null), the null being the same model with no parameter varying across conditions. Alternatives: the broadcast constant or the input scale varies instead. The ketamine-epoch variance explained is the same ratio on the ketamine half.
cid = 'fig5g-mouse-first-order-model'
rows = []
for key in (MAIN5E, CODE5E):
d, sig, tt = TAB5E[key]
ok = ~np.isnan(sig['pre'][0])
data = np.hstack([sig['pre'][:, ok].mean(0), sig['ket'][:, ok].mean(0)])
n = len(data) // 2
for kind in ('persist', 'broadcast', 'saturation'):
r = m.compare_model(kind, data)
ve_ket = 1 - np.mean((data[n:] - r['pred'][n:]) ** 2) / np.mean((data[n:] - r['pred_null'][n:]) ** 2)
rows.append(dict(dimension=f'{key[0]} s, var {key[1]}', model=kind, varexp=r['varexp'], varexp_ketamine_half=ve_ket,
S=r['params'][0], tB=r['params'][1], tP=r['params'][2], c_pre=r['params'][3], c_ket=r['params'][4]))
TAB5G = pd.DataFrame(rows); print(TAB5G.to_string())
v, tol = 0.954, 0.095
print(f"neural fit, main dimension: {TAB5G.varexp.iloc[0]:.3f} (expected {v} ± {tol}; within: {within(TAB5G.varexp.iloc[0], v, tol)})")
dimension model varexp varexp_ketamine_half S tB tP c_pre c_ket 0 5-6 s, var M persist 0.8815 0.9180 2.1236 1.1878 31.7484 0.9981 0.0099 1 5-6 s, var M broadcast 0.8709 0.8962 2.1475 47.1032 31.6427 0.0257 1.0000 2 5-6 s, var M saturation 0.8712 0.8967 2.1352 1.2010 31.6897 1.0000 0.0449 3 3-4 s, var C persist 0.9468 0.9574 2.1497 0.7619 27.7834 1.0000 0.0067 4 3-4 s, var C broadcast 0.9410 0.9453 2.0819 86.8257 27.9601 0.0081 1.0000 5 3-4 s, var C saturation 0.9412 0.9457 2.5225 0.7013 27.9761 0.8241 0.0171 neural fit, main dimension: 0.881 (expected 0.954 ± 0.095; within: True)
# analysis.py lines that define the models compared above: the fit and its null (compare_model), the persistence model and its
# two alternatives (theoretical_persist, theoretical_broadcast, theoretical_saturation), and the bounds (compare_models)
m.show_code(CODE, 2813, 2852)
m.show_code(CODE, 2867, 2882)
2813: def compare_model(model, data, p0, p0_null, bounds, bounds_null):
2814: # Fit model and its null, return varexp
2815: popt_m, pcov_m = curve_fit(model,
2816: np.arange(len(data)),
2817: data,
2818: p0=p0,
2819: maxfev=5000,
2820: bounds=bounds)
2821: pred = model(1, *popt_m)
2822: pred_err = np.mean((data-pred)**2)
2823: popt_null, pcov_null = curve_fit(model,
2824: np.arange(len(data)),
2825: data,
2826: p0=p0_null,
2827: maxfev=5000,
2828: bounds=bounds_null)
2829: pred_null = model(1, *popt_null)
2830: null_err = np.mean((data-pred_null)**2)
2831: return {"varexp": 1-(pred_err/null_err), "fit": popt_m, "trace": pred, "model null": pred_null}
2832:
2833: def theoretical_persist(x, S, tB, tP, *kwargs):
2834: traces = []
2835: for c in kwargs:
2836: output = model_system(S, tB, tP*c)
2837: traces.append(np.hstack([output[75:300*8], output[(300*8)+200:(300*8)+int(49/bin_t)]]))
2838: return np.hstack(traces)
2839:
2840: def theoretical_broadcast(x, S, tB, tP, *kwargs):
2841: traces = []
2842: for c in kwargs:
2843: output = model_system(S, tB*c, tP)
2844: traces.append(np.hstack([output[75:300*8], output[(300*8)+200:(300*8)+int(49/bin_t)]]))
2845: return np.hstack(traces)
2846:
2847: def theoretical_saturation(x, S, tB, tP, *kwargs):
2848: traces = []
2849: for c in kwargs:
2850: output = model_system(S*c, tB, tP)
2851: traces.append(np.hstack([output[75:300*8], output[(300*8)+200:(300*8)+int(49/bin_t)]]))
2852: return np.hstack(traces)
2867: def compare_models(model_dict, data, conditions=2):
2868: """ Compare different models against the data.
2869: - Model dict is a dictionary of model functions by name
2870: - Returns the variance explained from each model fit to the data, vs a null model with no variable fit across conditions.
2871: """
2872: var_exp_dict = {}
2873: model_names = model_dict.keys()
2874:
2875: from tqdm import tqdm
2876: for model_name in tqdm(model_names):
2877: p0 = [1,1,1] + [1]*conditions
2878: p0_null = [1,1,1] + [1]*conditions
2879: bounds = [[0,0,0] + [0]*conditions,[100,100,100]+[1]*conditions]
2880: bounds_null = [[0,0,0]+[0.99999]*conditions,
2881: [100,100,100]+[1]*conditions]
2882:
Behavioural part: eye closure of the dose experiments in DANDI:001347 (sub-mouse_ses-<id>_behavior.nwb). Which sessions received which dose is not in the dataset; the lists are read at run time from the paper's code record (Zenodo 14920570, mouse_eyepuff/scripts/fig1_dose.py, experiments), and that script's just_ket and just_pcp orderings (ketamine 6, 13, 25, 35 mg/kg; PCP 1, 5, 10 mg/kg) are taken as the conditions of Fig. 5H and 5I. Variants add saline and the highest dose.
import io, zipfile, hashlib, ast, urllib.request
url = 'https://zenodo.org/api/records/14920570/files/kauvar-richman-liu-2025_zenodo.zip/content'
raw = urllib.request.urlopen(url, timeout=120).read()
print('md5', hashlib.md5(raw).hexdigest(), 'expected bd1fb567121aaf0291337783df57bcfc')
assert hashlib.md5(raw).hexdigest() == 'bd1fb567121aaf0291337783df57bcfc'
zf = zipfile.ZipFile(io.BytesIO(raw))
src = zf.read([n for n in zf.namelist() if n.endswith('mouse_eyepuff/scripts/fig1_dose.py')][0]).decode()
DOSE = None
for node in ast.walk(ast.parse(src)):
if isinstance(node, ast.Assign) and getattr(node.targets[0], 'id', None) == 'experiments':
DOSE = ast.literal_eval(node.value)
print(DOSE)
t0 = time.time()
DTR, DMICE = {}, {}
for k, ids in DOSE.items():
if k != 'kx':
DTR[k], DMICE[k] = m.dose_trace(ids)
print(f'{sum(len(v) for k, v in DOSE.items() if k != "kx")} sessions read in {time.time() - t0:.0f} s')
print('late eye closure 0.3-0.8 s after the last puff (mean of mice):',
{k: round(float(v[100 + 7 * 300 + 30:100 + 7 * 300 + 80].mean()), 3) for k, v in DTR.items()})
rows = []
for name, keys in {'ketamine (code just_ket)': ['ket6', 'ket13', 'ket25', 'ket35'],
'PCP (code just_pcp)': ['pcp1', 'pcp5', 'pcp10'],
'ketamine + 50': ['ket6', 'ket13', 'ket25', 'ket35', 'ket50'],
'PCP + 20': ['pcp1', 'pcp5', 'pcp10', 'pcp20'],
'ketamine + saline + 50': ['sal', 'ket6', 'ket13', 'ket25', 'ket35', 'ket50'],
'PCP + saline + 20': ['sal', 'pcp1', 'pcp5', 'pcp10', 'pcp20']}.items():
r = m.fit_behaviour([DTR[k] for k in keys])
rows.append(dict(fit=name, varexp=r['varexp'], S=r['params'][0], tB=r['params'][1], tP=r['params'][2],
multipliers=np.round(r['params'][3:], 3).tolist()))
TAB5H = pd.DataFrame(rows); print(TAB5H.to_string())
vk, vp = TAB5H.varexp.iloc[0], TAB5H.varexp.iloc[1]
for a, (ek, ep) in {'ketamine 0.927, PCP 0.870': (0.927, 0.870), 'ketamine 0.870, PCP 0.927': (0.870, 0.927)}.items():
print(f'assignment {a}: ketamine {vk:.3f} within {within(vk, ek, 0.1 * ek)}, PCP {vp:.3f} within {within(vp, ep, 0.1 * ep)}')
md5 bd1fb567121aaf0291337783df57bcfc expected bd1fb567121aaf0291337783df57bcfc
{'ket50': [275, 276, 277, 280, 281], 'pcp20': [287, 288, 289, 290, 291], 'sal': [282, 283, 284, 285, 286], 'kx': [316, 338, 339, 343, 344], 'ket25': [355, 356, 357, 362, 363], 'pcp10': [358, 359, 360, 361, 364], 'ket13': [401, 400, 399, 398, 397], 'ket6': [403, 405, 404, 406, 426], 'pcp5': [370, 368, 366, 365, 434], 'pcp1': [430, 432, 431, 435, 436], 'ket35': [1020, 1021, 1022, 1023, 1024, 1025, 1026, 1027]}
53 sessions read in 61 s
late eye closure 0.3-0.8 s after the last puff (mean of mice): {'ket50': -0.014, 'pcp20': -0.004, 'sal': 0.595, 'ket25': 0.187, 'pcp10': 0.103, 'ket13': 0.275, 'ket6': 0.388, 'pcp5': 0.287, 'pcp1': 0.229, 'ket35': 0.079}
fit varexp S tB tP multipliers 0 ketamine (code just_ket) 0.9298 25.1549 100.0000 21.1188 [0.996, 0.359, 0.182, 0.024] 1 PCP (code just_pcp) 0.8704 0.6455 3.0734 16.5658 [1.0, 0.582, 0.116] 2 ketamine + 50 0.9413 25.1524 99.9998 21.0395 [1.0, 0.361, 0.182, 0.024, 0.0] 3 PCP + 20 0.9188 0.6454 3.0728 16.5607 [1.0, 0.582, 0.116, 0.001] 4 ketamine + saline + 50 0.8850 44.8148 99.9999 13.7863 [1.0, 0.757, 0.291, 0.15, 0.027, 0.001] 5 PCP + saline + 20 0.8838 54.7878 99.9998 11.4342 [0.986, 0.378, 0.239, 0.066, 0.0] assignment ketamine 0.927, PCP 0.870: ketamine 0.930 within True, PCP 0.870 within True assignment ketamine 0.870, PCP 0.927: ketamine 0.930 within True, PCP 0.870 within True
Result: reproduced. The persistence model explains 0.881 of the variance of the across-session emotion-like signal of the main dimension (0.947 with the code's window), within 0.954 ± 0.095; letting the broadcast constant or the input scale vary instead explains less (0.871 and 0.871 overall, 0.896 and 0.897 against 0.918 on the ketamine half; 0.941, 0.941 against 0.947 with the code's window). The dose fits give 0.930 for ketamine (6, 13, 25, 35 mg/kg) and 0.870 for PCP (1, 5, 10 mg/kg), within 10% of 0.927 and 0.870 under the assignment of 0.927 to ketamine (Fig. 5H) and 0.870 to PCP (Fig. 5I); the 10% tolerance also admits the other assignment (0.930 against 0.870, 0.870 against 0.927), so the panels cannot be told apart this way. Adding saline and the highest dose lowers both (0.885 and 0.884). The printed code lines define the models compared: the null is the same model with both condition multipliers bounded to [0.99999, 1] (compare_models), and the alternatives apply the multiplier to the broadcast constant tB or to the input S in place of tP (theoretical_broadcast, theoretical_saturation).
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| claim | outcome | main numbers |
|---|---|---|
| fig6d-mouse-persistent-activity-reduced | reproduced (assumed epochs) | 0.867, -0.002, 0.515; corrected p 0.0005 and 0.029 |
| fig6f-mouse-fast-activity-unchanged | not reproduced (absence claim; the ketamine part holds) | pre vs ketamine corrected p 0.67; ketamine vs post 0.0045 |
| fig5e-mouse-emotion-dimension-late-reduced | not reproduced (decided on the criterion's 5-6 s target) | late corrected p 0.086 (5-6 s); 0.039 with the 3-4 s window inferred from the code |
| fig6h-mouse-timescale-shorter | reproduced | 0.437 to 0.162 s, p 0.0057 |
| fig6j-mouse-coupling-reduced | partially reproduced | ratio 0.852, p 0.30 (Methods; mean signed r median 0.023, ratio 0.576 to 1.450 across readings); 0.559, p 8.7e-05 (code) |
| fig5d-mouse-emotion-dimension-accumulates | reproduced | accumulation 0.846 pre, -0.048 ketamine; 10-20 s 0.838 |
| figs34-mouse-tone-persistence-reduced | partially reproduced | fast n.s., persistent and timescale reduced; coupling over tone windows not within tolerance under either rule (1.238 Methods; 0.664 to 1.019 code), 0.535 only over puff windows |
| fig4hi-mouse-rise-unchanged-decay-shortened | partially reproduced | decay -0.189 s (CI below 0); rise -0.009 s (CI excludes 0); 0 of 8 readings satisfy both parts |
| fig4g-mouse-slow-clusters-altered | reproduced | Spearman rho 0.47; slowest 3 0.507 vs fastest 3 0.368 |
| fig5g-mouse-first-order-model | reproduced | neural 0.881; ketamine 0.930, PCP 0.870 |
Differences from the paper's method:
- Session selection (assumption): mouse 5 is taken as the animal the Methods exclude, the only choice that gives both 11 sessions from 9 mice and 7 tone sessions from 7 mice. The 13-session main readings are printed for fig6d, fig6f, fig5e, fig5d, fig6h, fig6j and figs34 (9 tone sessions).
- Session z-scoring of the rates is assumed (the code's driver that makes
all_zfris not in the record). It does not affect the z-scored projections of Figs. 6D, 6F and S34. - Coding-dimension variance: the Methods' per-neuron variance is the main reading, the code's pooled scalar a variant. It decides fig6j and the coupling part of figs34, and moves the means of fig6f and fig5e.
- Epochs: the Methods' trials 22 to 32 and 34 to 40 disagree with Table S3's counts (10 and 6) and the code's slices (23 to 32, 35 to 40); the code's slices are adopted for fig6d, fig6f, fig5e, fig4g and fig4hi.
- Emotion-like dimension: the Methods put the target 5 s after the last puff, the code (as its plotting offsets place it; the code lines are printed in the fig5e section) 3 to 4 s after; the criterion's 5 to 6 s is the main reading and decides the outcome. Only the code's window reproduces the paper's Fig. 5E values.
- Tone coupling: the code takes the windows before puff series 1 to 10 and 23 to 32 for every stimulus (printed in the figs34 section); the Methods' and Table S3's tone windows are the main reading, run under both the Methods' and the code's rule.
- Timescale and coupling: the code's series (11 to 20 and 23 to 32 for the timescale, 1 to 10 and 23 to 32 for coupling) are used, as Table S3's trial counts agree with them; fits from lag 0 and per-series averaging follow the Methods, the code's lag one and concatenation are variants.
- Tone epochs: split at the injection time (10 preinfusion tone series, as Table S3) as the main reading; the code's split gives 9.
- Clustering: Scanpy 1.12 with leidenalg gives 11 clusters, the paper 12; no cluster is dropped (the code drops its cluster 11 in Fig. 4H, I).
- fig5g dose fits: the conditions are the code's
just_ketandjust_pcpdose lists; the dose session lists come from the code record.
What was found about the dataset:
- DANDI:001326 holds 13 sessions from 10 mice, as the paper states; units are pre-curated (all 'good', SNR above 2, more than 500 spikes); no ISI metric is stored.
- Spike times are stored unsorted within every unit.
- The subjects are numbered 0 to 9 with no link to the paper's mouse or experiment identifiers; the link to DANDI:001347's Neuropixels-rig behaviour files was made here from the inter-series intervals (10 of 13 sessions; sub-3 ses-20230225 and both sub-5 sessions have no behaviour file).
- The excluded animal is not named in the paper, the code or the dataset.
- sub-5 ses-20230706 has 2 probes, where the Methods say 3 to 4.