Testing Kutschireiter et al. (2023) with dandiset 000289¶
Paper. Kutschireiter A, Basnak MA, Wilson RI, Drugowitsch J (2023). Bayesian inference in ring attractor networks. PNAS 120, e2210622120. doi:10.1073/pnas.2210622120. The claims file was written from the version of record (Europe PMC, PMC9992764) with its SI Appendix.
Dandiset. DANDI:000289, version 0.240801.1040, "Two-photon calcium imaging of EPG cells in Drosophila Melanogaster presented with cues of varying reliability" (Basnak and Wilson). The paper cites no dandiset version. This is the only published version and the most recent; the draft has the same 221 assets and size, so the published version is pinned.
What the dandiset holds, and what the paper uses. The paper is a theory paper: every result is a simulation, apart from a connectome comparison and one prediction. The dandiset holds two-photon imaging of EPG neurons in the protocerebral bridge, with the cue position and trial tables, from a separate study (Basnak et al., "Uncertainty and learning in a biological ring attractor", named in the dandiset's description). So, as the claims file sets out:
- twelve claims are
data: elsewhere: they are tested on the paper's own code and processed outputs, Zenodo doi:10.5281/zenodo.7615975 (BayesianRingAttractor v1.0, commit a9db6f9), and, for Fig. 4B, on the hemibrain connectome; - one claim, the Discussion's prediction that bump amplitude is higher with a reliable cue, is tested on the dandiset: 15 flies shown a visual cue at zero, low and high contrast.
The dandiset's other 206 files (subjects 16 to 68) hold other experiments of that study: a cue rotating in open loop, a cue of inverted gain, visual and wind cues in closed loop, and cues shifted every 5 minutes (from the session descriptions and trial tables, read file by file while preparing this attempt). No claim uses them.
All sources are read into memory; nothing is downloaded to disk.
| Claim | Data | Tested |
|---|---|---|
| fig1-conflict-decreases-certainty | paper's code | tested |
| fig2h-bump-tracks-circkf | paper's code (Figure2h.ipynb, run) | tested |
| fig2i-accuracy-and-certainty-vs-reliability | figure2i.npz | tested |
| fig3b-bayesian-beats-conventional | figure2i.npz | tested |
| fig3c-update-weights | figure2i.npz | tested |
| fig3d-best-network-slow-decay | figure3d.npz | tested |
| fig3f-noise-drop-below-70 | figure3f/sigma=*.npz | tested |
| fig3f-retuning-reduces-drop | figure3f/sigma=*.npz | tested |
| noise-effects-on-mu-and-kappa | paper's code (RNN_filter_Run) | tested |
| fig4b-connectome-motifs | hemibrain v1.2 export | tested (sorting reading (a) only; reading (b) cannot be formed) |
| fig4gh-drosophila-network-matches | figure4g.npz | tested |
| figs2-performance-independent-of-n | figureS2.npz | tested |
| prediction-bump-amplitude-cue-reliability | dandiset 000289 | tested, then blocked (methods_underspecified): its readings disagree |
Every claim is tested. One, the prediction, is then blocked as methods_underspecified under the
claims file's own rule, because its readings give different outcomes (see its section).
import time
import numpy as np
import pandas as pd
from scipy import stats
from concurrent.futures import ThreadPoolExecutor
import bra
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
np.set_printoptions(linewidth=200, suppress=True)
started = time.time()
code = bra.Code() # the paper's Zenodo deposit, in memory, md5 checked
flt = code.flt # circularFiltering.py, with the patches listed below
R = bra.RATES
print("Zenodo deposit doi:10.5281/zenodo.7615975, files:", len(code.names))
print("\n".join(n for n in code.names if n.endswith((".py", ".ipynb", ".npz"))))
print("\nPatches to the released code:", *code.patches, sep="\n- ")
Zenodo deposit doi:10.5281/zenodo.7615975, files: 49 Figure2a,c_FigureS1.ipynb Figure2h.ipynb Figure2i.ipynb Figure3b,c,e,d_FigureS3.ipynb Figure3f_FigureS4.ipynb Figure4g,h.ipynb FigureS2.ipynb circularFiltering.py data_processed/figure2i.npz data_processed/figure3d.npz data_processed/figure3f/sigma=0.0.npz data_processed/figure3f/sigma=1.0.npz data_processed/figure3f/sigma=10.0.npz data_processed/figure3f/sigma=15.0.npz data_processed/figure3f/sigma=2.0.npz data_processed/figure3f/sigma=20.0.npz data_processed/figure3f/sigma=3.0.npz data_processed/figure3f/sigma=4.0.npz data_processed/figure3f/sigma=5.0.npz data_processed/figure3f/sigma=6.0.npz data_processed/figure3f/sigma=7.0.npz data_processed/figure3f/sigma=8.0.npz data_processed/figure3f/sigma=9.0.npz data_processed/figure4g.npz data_processed/figureS2.npz data_processed/vonMises_q_stats.npz network_filter.py performance_scan_figure2i.py performance_scan_figure3d.py performance_scan_figure3f.py performance_scan_figure4g.py performance_scan_figureS2.py Patches to the released code: - RNN_filter_Run: added z=None to the signature (it uses z and has no such parameter; claims file header) - import of isscalar from numpy, not numpy.core.numeric (numpy 2)
Data check¶
The paper's processed outputs, the connectome export, and the 15 imaging sessions of the cue-contrast experiment are loaded and described. Each check of a description against the values is stated below the output.
# --- the paper's processed outputs
for name in ["figure2i.npz", "figure3d.npz", "figure4g.npz", "figureS2.npz", "vonMises_q_stats.npz",
"figure3f/sigma=0.0.npz"]:
z = code.npz(name)
print(name, {k: z[k].shape for k in z.files})
f2i = code.npz("figure2i.npz")
vm, vq, pf, nu = f2i["vonMises_stats"], f2i["vonMises_q_stats"], f2i["PF_stats"], f2i["noUncert_stats"]
print("\nfigure2i.npz: axis 0 is the 33 kappa* of the fixed-certainty filter, axis 1 the 23 information rates",
len(bra.KAPPA_FIXED), len(R))
print("accuracy |m1| range over all filters and rates: %.3f to %.3f (must lie in 0..1)" %
(min(vm[0, :, 1].min(), vq[0, :, 1].min(), nu[:, :, 1].min()),
max(vm[0, :, 1].max(), vq[0, :, 1].max(), nu[:, :, 1].max())))
print("mean final certainty, circKF: %.2f to %.2f; Bayesian ring attractor: %.2f to %.2f" %
(vm[0, :, 2].min(), vm[0, :, 2].max(), vq[0, :, 2].min(), vq[0, :, 2].max()))
ref_pf = bra.objective(pf[0, :, 1])
print("particle filter reference, objective(PF accuracy) = %.10f; Figure3f_FigureS4.ipynb hard-codes %.10f"
% (ref_pf, bra.REF_PF_3F))
print("rows of figure2i.npz beyond the first, vonMises_stats all zero:", bool(np.all(vm[1:] == 0)))
figure2i.npz {'vonMises_stats': (33, 23, 3), 'vonMises_q_stats': (33, 23, 3), 'PF_stats': (33, 23, 3), 'noUncert_stats': (33, 23, 2)}
figure3d.npz {'stats': (23, 3), 'MSE_array': (25, 13)}
figure4g.npz {'vonMises_stats': (23, 3), 'vonMises_q_stats': (23, 3), 'CX_stats': (23, 3)}
figureS2.npz {'MSE_array': (7, 23)}
vonMises_q_stats.npz {'stats': (23, 3)}
figure3f/sigma=0.0.npz {'MSE_array': (40, 13), 'avg_network_activity': (40, 13)}
figure2i.npz: axis 0 is the 33 kappa* of the fixed-certainty filter, axis 1 the 23 information rates 33 23
accuracy |m1| range over all filters and rates: 0.004 to 0.965 (must lie in 0..1)
mean final certainty, circKF: 0.18 to 13.60; Bayesian ring attractor: 1.00 to 13.58
particle filter reference, objective(PF accuracy) = 0.6912447327; Figure3f_FigureS4.ipynb hard-codes 0.6912447327
rows of figure2i.npz beyond the first, vonMises_stats all zero: True
# Where the arrays and parameters used below are defined in the released code (lines printed from
# the deposit): the cue rate of Figure 2H, the columns of the stats arrays, the grids of Fig. 3, and
# the number of runs per setting of each scan.
for name, pat in [("Figure2h.ipynb", r"^kappa_z = |alpha = flt.xi_fun_inv|kappa_z_array\[ind_off"),
("circularFiltering.py", r"def circular_mean|return phi_hat, r_hat"),
("Figure2i.ipynb", r"_stats\[i,j,|_stats\[i,j\] =|noUncert_stats = np.zeros|vonMises_stats = np.zeros"),
("Figure3b,c,e,d_FigureS3.ipynb", r"^kappa_array = |^\s+1\.6,|^\s+2\.5,|^\s+7\. |^kappa_star_array|^beta_array|MSE_array\[i,j\] =|MSE_array = np.array\(\[ objective\(noUncert"),
("Figure3f_FigureS4.ipynb", r"^kappa_star_array|^\s+1\.1, 1\.2|^\s+2\.0, 2\.1|^\s+3\.0, 3\.2|^beta_array|^sigmaN_array|^N = |MSE_array\[i,j\] =|^ref = |max_MSE\[ind_sig\] =|noretuning_MSE\[ind_sig\] =|idx_ref = idx"),
("FigureS2.ipynb", r"^N_array|MSE_array\[i\] ="),
("Figure4g,h.ipynb", r"_stats\[i,\[0,1\]\]|_stats\[i,2\]")]:
print("== " + name)
print("\n".join(" " + ln for ln in code.source_lines(name, pat)))
for name in ["performance_scan_figure2i.py", "performance_scan_figure3d.py", "performance_scan_figure3f.py",
"performance_scan_figure4g.py", "performance_scan_figureS2.py"]:
print("== %s: %s" % (name, "; ".join(ln.strip() for ln in code.source_lines(name, r"^nIter =|^T = |^dt = "))))
== Figure2h.ipynb
kappa_z = 15 # precision of absolute heading info (called gamma_z in manuscript)
alpha = flt.xi_fun_inv(kappa_z*dt)
kappa_z_array[ind_off:ind_on] = 0 # maximally uninformative while light is off
== circularFiltering.py
def circular_mean(phi,w=None):
return phi_hat, r_hat
== Figure2i.ipynb
vonMises_stats = np.zeros([len_k,len_kz,3]) # 1st (bias + variance) and 2nd order statistics for each kappa_z,
noUncert_stats = np.zeros([len_k,len_kz,2])
vonMises_stats[i,j,2] = np.mean(vonMises[:,1])
vonMises_q_stats[i,j,2] = np.mean(vonMises_q[:,1])
PF_stats[i,j,2] = np.mean(PF[:,1])
vonMises_stats[i,j,[0,1]] = flt.circular_mean(flt.backToCirc(vonMises[:,0]-phi_final))
vonMises_q_stats[i,j,[0,1]] = flt.circular_mean(flt.backToCirc(vonMises_q[:,0]-phi_final))
PF_stats[i,j,[0,1]] = flt.circular_mean(flt.backToCirc(PF[:,0]-phi_final))
noUncert_stats[i,j] = flt.circular_mean(flt.backToCirc(noUncert - phi_final))
noUncert_stats[i,j] = flt.circular_mean(flt.backToCirc(noUncert - phi_final))
== Figure3b,c,e,d_FigureS3.ipynb
kappa_array = np.array([ 0.1, 0.3, 0.9, 1. , 1.1, 1.2, 1.3, 1.4, 1.5,
1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2, 2.3, 2.4,
2.5, 2.6, 2.7, 2.8, 2.9, 3. , 4. , 5. , 6. ,
7. , 8. , 9. , 10. , 30. , 100. ])
MSE_array = np.array([ objective(noUncert_stats[i,:,1],mu=mu_perf,sigma=sig_perf) for i in np.arange(len(kappa_array))])
kappa_array = np.arange(0.2,10,0.01) # own uncertainty right before update
kappa_array = np.arange(0.2,10,0.01) # own uncertainty right before update
kappa_star_array = np.array([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5])
beta_array = np.array([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 1, 2 ,4 ,8])
MSE_array[i,j] = objective( stats[:,1] ,mu=mu_perf,sigma=sig_perf)
== Figure3f_FigureS4.ipynb
ref = 0.6912447327304075
kappa_star_array = np.array([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,
1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9,
3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0])
beta_array = np.array([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 1, 2 ,4 ,8])
N = 64
sigmaN_array = np.array([0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,15,20])
MSE_array[i,j] = objective( stats[:,1] ,mu=mu_perf,sigma=sig_perf)
idx_ref = idx
max_MSE[ind_sig] = np.nanmax(MSE_array)
noretuning_MSE[ind_sig] = MSE_array[idx_ref]
== FigureS2.ipynb
N_array = np.array([3,4,8,16,32,64,128])
MSE_array[i] = stats[:,1]
== Figure4g,h.ipynb
vonMises_stats[i,2] = np.mean(vonMises[:,1])
vonMises_q_stats[i,2] = np.mean(vonMises_q[:,1])
CX_stats[i,2] = np.mean(data['CX'][:,1])
vonMises_stats[i,[0,1]] = flt.circular_mean(flt.backToCirc(vonMises[:,0]-phi_final))
vonMises_q_stats[i,[0,1]] = flt.circular_mean(flt.backToCirc(vonMises_q[:,0]-phi_final))
CX_stats[i,[0,1]] = flt.circular_mean(flt.backToCirc(CX[:,0]-phi_final))
== performance_scan_figure2i.py: nIter = 10000; T = 20 # simulation time; dt = 0.01 # step size
== performance_scan_figure3d.py: nIter = 10000; T = 20 # simulation time; dt = 0.01 # step size
== performance_scan_figure3f.py: nIter = 10000; T = 20 # simulation time; dt = 0.01 # step size
== performance_scan_figure4g.py: nIter = 10000; T = 20 # simulation time; dt = 0.001 # step size
== performance_scan_figureS2.py: nIter = 10000; T = 20 # simulation time; dt = 0.01 # step size
The processed outputs have the shapes the figure notebooks expect, accuracies lie in 0 to 1, and the
particle filter reference computed from figure2i.npz equals the constant hard-coded in
Figure3f_FigureS4.ipynb, so the two notebooks normalize alike. The certainty stored for the
Bayesian ring attractor never falls below 1, its baseline kappa* (Methods), while the circKF's falls
to 0.18.
The lines printed from the released code fix what the arrays mean. In Figure2h.ipynb,
kappa_z = 15 is commented "called gamma_z in manuscript" and enters as xi_fun_inv(kappa_z*dt).
Every *_stats array holds circular_mean of the end errors in columns 0 and 1 (phase, then
r_hat, the accuracy |m1|) and, where present, the mean certainty in column 2; noUncert_stats
has the first two only, with axis 0 the 33 kappa* of kappa_array. In figure3d.npz and the
figure3f files, axis 0 is kappa_star_array and axis 1 beta_array, as the loops
MSE_array[i,j] show; the rows of figureS2.npz follow N_array. Every scan averages 10,000 runs
of 20 s per setting (dt 0.01; 0.001 for Fig. 4G), which sets the Monte Carlo error used in the
sections below: the standard error of an accuracy near 0.2 (0.1/s) is about 0.007, and that of a
difference between two independent scans about 0.010.
# --- the hemibrain connectome export (v1.2): neurons of the three types and their names
neurons, conns = bra.hemibrain()
types = {"EPG": "EPG", "PEN_a(PEN1)": "PEN1", "Delta7": "Delta7"}
cx = neurons[neurons.type.isin(types)].copy()
cx["t"] = cx.type.map(types)
cx["g"] = cx.instance.map(bra.pb_glomerulus)
print("traced neurons:", len(neurons), "; connections:", len(conns), "; total synapses:", conns.weight.sum())
print(cx.groupby("t").agg(n=("bodyId", "size"), with_glomerulus=("g", lambda s: s.notna().sum())).to_string())
for t in ["EPG", "PEN1", "Delta7"]:
print(t, "instances:", sorted(cx[cx.t == t].instance.unique())[:6], "...")
print("glomeruli named: EPG", sorted(cx[cx.t == 'EPG'].g.unique(), key=bra.PB_ORDER.index),
"\n PEN1", sorted(cx[cx.t == 'PEN1'].g.unique(), key=bra.PB_ORDER.index))
traced neurons: 21739 ; connections: 3550403 ; total synapses: 14329229
n with_glomerulus
t
Delta7 42 0
EPG 46 46
PEN1 20 20
EPG instances: ['EPG(PB08)_L1', 'EPG(PB08)_L2', 'EPG(PB08)_L3', 'EPG(PB08)_L4', 'EPG(PB08)_L5', 'EPG(PB08)_L6'] ...
PEN1 instances: ['PEN_a(PB06a)_L2', 'PEN_a(PB06a)_L3', 'PEN_a(PB06a)_L4', 'PEN_a(PB06a)_L5', 'PEN_a(PB06a)_L6', 'PEN_a(PB06a)_L7'] ...
Delta7 instances: ['Delta7(PB15)_L1L9R8_R', 'Delta7(PB15)_L2R7_R', 'Delta7(PB15)_L3R6_R', 'Delta7(PB15)_L4R5_R', 'Delta7(PB15)_L4R6_R', 'Delta7(PB15)_L5R4_L'] ...
glomeruli named: EPG ['L8', 'L7', 'L6', 'L5', 'L4', 'L3', 'L2', 'L1', 'R1', 'R2', 'R3', 'R4', 'R5', 'R6', 'R7', 'R8']
PEN1 ['L9', 'L8', 'L7', 'L6', 'L5', 'L4', 'L3', 'L2', 'R2', 'R3', 'R4', 'R5', 'R6', 'R7', 'R8', 'R9']
46 EPG, 20 PEN1 and 42 Delta7 neurons. Every EPG and PEN1 instance names one protocerebral bridge glomerulus (EPG L8 to R8, PEN1 L9 to L2 and R2 to R9), which is reading (a) of the claim's sort. No instance names an ellipsoid body or gall position, so reading (b) cannot be formed from the names. Delta7 names list the glomeruli the neuron spans, not one position.
# --- dandiset 000289: every asset, and the cue-contrast sessions
assets = bra.assets()
print("dandiset", bra.DANDISET, "version", bra.VERSION, ": assets", len(assets),
"; total size %.0f GB" % (sum(a["size"] for a in assets) / 1e9))
subj = lambda a: int(a["path"].split("/")[0][4:])
contrast_assets = sorted([a for a in assets if subj(a) <= 15], key=subj)
def load(a):
f = bra.open_nwb(a)
s = bra.read_session(f)
s["path"] = a["path"]
return s
with ThreadPoolExecutor(8) as ex:
sessions = list(ex.map(load, contrast_assets))
print("loaded %d sessions in %.0f s" % (len(sessions), time.time() - started))
rows = []
for s in sessions:
tr = s["trials"]
dt = np.diff(s["t"])
rows.append(dict(sub=s["subject"], description=s["description"][:45], rois=s["dff"].shape[1],
frames=s["dff"].shape[0], planes=s["tps_shape"][1], dur_s=round(s["t"][-1] - s["t"][0], 1),
frame_hz=round(1 / dt.mean(), 2), meta_hz=s["rate_meta"],
dff_min=round(s["dff"].min(), 2), dff_max=round(s["dff"].max(), 2),
cue_min=round(s["cue"].min(), 3), cue_max=round(s["cue"].max(), 3),
same_clock=bool(np.array_equal(s["t"], s["t_cue"]) and np.array_equal(s["t"], s["t_tps"])),
blocks="".join(str(int(c)) for c in tr["cue_contrast"]),
block_s=" ".join(str(round(b - a)) for a, b in zip(tr["start_time"], tr["stop_time"]))))
info = pd.DataFrame(rows)
print(info.to_string(index=False))
print("\nsession descriptions:", info.description.unique())
print("trial column descriptions:", sessions[0]["trial_desc"])
print("cue position unit:", sessions[0]["cue_unit"])
# recordings deposited twice? compare the dF/F of every pair of sessions with the same shape
dups = [(a["sub"], b["sub"]) for i, a in enumerate(rows) for b in rows[i + 1:]
if sessions[i]["dff"].shape == sessions[[r["sub"] for r in rows].index(b["sub"])]["dff"].shape
and np.allclose(sessions[i]["dff"], sessions[[r["sub"] for r in rows].index(b["sub"])]["dff"])]
npix = [int((m > 0).sum()) for s in sessions for m in s["masks"]]
print("ROI mask sizes, pixels with nonzero weight per ROI: min %d, median %g, max %d over %d ROIs; mask array shape of sub-1:"
% (min(npix), np.median(npix), max(npix), len(npix)), sessions[0]["masks"].shape)
print("identical dF/F pairs among the 15 sessions:", dups)
print("identical cue traces:", [(rows[i]["sub"], rows[j]["sub"]) for i in range(15) for j in range(i + 1, 15)
if sessions[i]["cue"].shape == sessions[j]["cue"].shape
and np.allclose(sessions[i]["cue"], sessions[j]["cue"])])
dandiset 000289 version 0.240801.1040 : assets 221 ; total size 761 GB
loaded 15 sessions in 35 s
sub description rois frames planes dur_s frame_hz meta_hz dff_min dff_max cue_min cue_max same_clock blocks block_s
1 fly presented with visual cues of varying con 37 11020 12 1200.0 9.18 9.18 -0.47 3.90 0.002 6.281 True 213123 195 200 200 200 200 205
2 fly presented with visual cues of varying con 38 11020 12 1200.0 9.18 9.18 -0.81 4.33 0.000 6.283 True 31212 195 200 200 200 405
3 fly presented with visual cues of varying con 39 11020 12 1200.0 9.18 9.18 -0.40 4.12 0.000 6.283 True 121323 195 200 200 200 200 205
4 fly presented with visual cues of varying con 41 9446 14 1200.0 7.87 9.18 -0.37 3.64 0.000 6.283 True 121323 195 200 200 200 200 205
5 fly presented with visual cues of varying con 42 11020 12 1200.0 9.18 9.18 -0.82 3.72 0.001 6.281 True 312123 195 200 200 200 200 205
6 fly presented with visual cues of varying con 43 11020 12 1200.0 9.18 9.18 -0.44 4.25 0.000 6.283 True 213123 195 200 200 200 200 205
7 fly presented with visual cues of varying con 47 11020 12 1200.0 9.18 9.18 -0.84 4.36 0.001 6.282 True 121323 196 200 200 200 200 204
8 fly presented with visual cues of varying con 41 11020 12 1200.0 9.18 9.18 -0.40 3.30 0.000 6.283 True 312123 196 200 200 200 200 204
9 fly presented with visual cues of varying con 44 11020 12 1200.0 9.18 9.18 -0.53 5.24 0.001 6.281 True 213123 195 200 200 200 200 205
10 fly presented with visual cues of varying con 44 9446 14 1200.0 7.87 9.18 -0.59 7.87 0.002 6.283 True 213123 196 200 200 200 200 204
11 fly presented with visual cues of varying con 42 11020 12 1200.0 9.18 9.18 -0.89 3.67 0.000 6.282 True 312123 195 200 200 200 200 205
12 fly presented with visual cues of varying con 43 11020 12 1200.0 9.18 9.18 -0.87 3.99 0.001 6.283 True 121323 195 200 200 200 200 205
13 fly presented with visual cues of varying con 42 11020 12 1200.0 9.18 9.18 -0.46 5.02 0.000 6.283 True 312123 196 200 200 200 200 204
14 fly presented with visual cues of varying con 44 11020 12 1200.0 9.18 9.18 -0.53 5.84 0.000 6.283 True 121323 195 200 200 200 200 205
15 fly presented with visual cues of varying con 41 11020 12 1200.0 9.18 9.18 -0.51 5.72 0.000 6.283 True 213123 195 200 200 200 200 205
session descriptions: <StringArray>
['fly presented with visual cues of varying con']
Length: 1, dtype: str
trial column descriptions: {'cue_contrast': 'visual cue contrast (where 1 = zero contrast, 2 = low contrast, 3 = high contrast', 'start_time': 'start time of trial', 'stop_time': 'end of each trial'}
cue position unit: radians
ROI mask sizes, pixels with nonzero weight per ROI: min 3, median 7, max 11 over 628 ROIs; mask array shape of sub-1: (37, 256, 128)
identical dF/F pairs among the 15 sessions: []
identical cue traces: []
The cue-contrast experiment is subjects 1 to 15, one session each ("fly presented with visual cues
of varying contrast"), with a trials table whose cue_contrast is 1 = zero, 2 = low and
3 = high contrast. Each session is 1200 s of 37 to 47 dF/F traces along the PB midline. The dF/F,
the cue position and the raw imaging share one timestamp array. The cue position lies in 0 to 2 pi,
as its unit (radians) says. No session duplicates another. The ROI masks are 2D (one plane of 256 x 128 pixels) and
have 3 to 11 pixels each (median 7).
Two disagreements with the metadata: sub-4 and sub-10 have a frame interval of 0.127 s (7.87 Hz)
and 14 planes, while imaging_rate says 9.18 Hz for every file; and sub-2 has five blocks (high,
zero, low, zero, then 405 s labelled low), so it has one high-contrast block where every other fly
has two, and its last block is twice as long as the protocol's 200 s.
# --- does the bump follow the cue on the session's clock? Bump phase (frequency 2 across the PB)
# against the cue position, with the cue trace shifted in time, in each contrast.
def contrast_mask(s, c, shift_blocks=0):
tr = s["trials"]
labels = np.roll(tr["cue_contrast"], shift_blocks)
m = np.zeros(len(s["t"]), bool)
for a, b, l in zip(tr["start_time"], tr["stop_time"], labels):
if l == c:
m |= (s["t"] >= a) & (s["t"] < b)
return m
lags = np.arange(-30, 31, 5)
coh = {c: [] for c in (1, 2, 3)}
peak = []
for s in sessions:
phase, _ = bra.fourier(s["dff"])
for c in (1, 2, 3):
m = contrast_mask(s, c)
coh[c].append([np.abs(np.mean(np.exp(1j * (phase - np.roll(s["cue"], L)))[m])) for L in lags])
peak.append(lags[np.argmax(coh[3][-1])])
for c, name in [(1, "zero"), (2, "low"), (3, "high")]:
print("%-4s contrast, phase-cue coherence by lag (frames):" % name,
dict(zip(lags.tolist(), np.round(np.median(coh[c], 0), 3))))
print("lag of peak coherence in high contrast, per fly (frames):", peak)
# the trial table's clock: coherence in high-contrast frames with true and shifted block labels
true = [np.abs(np.mean(np.exp(1j * (bra.fourier(s["dff"])[0] - s["cue"]))[contrast_mask(s, 3)])) for s in sessions]
shift = [np.abs(np.mean(np.exp(1j * (bra.fourier(s["dff"])[0] - s["cue"]))[contrast_mask(s, 3, 1)])) for s in sessions]
print("median coherence in frames labelled high contrast: true labels %.3f, labels shifted one block %.3f"
% (np.median(true), np.median(shift)))
print("flies where the true labels give the higher coherence: %d of 15" % np.sum(np.array(true) > np.array(shift)))
# bra.fourier on a synthetic bump at ROI angle 1.0 rad over 40 ROIs (two periods): phase sign and amplitude
_n = 40
_ph, _amp = bra.fourier((3 * np.cos(2 * np.pi * 2 * np.arange(_n) / _n - 1.0))[None, :])
print("synthetic bump at +1.000 rad, amplitude 3: phase %+.3f, amplitude %.3f (the phase sign is a convention, see bra.fourier)"
% (_ph[0], _amp[0]))
print("ROIs per fly (frequency 2 assumes the ROIs span two bump periods):", [s["dff"].shape[1] for s in sessions])
zero contrast, phase-cue coherence by lag (frames): {-30: np.float64(0.078), -25: np.float64(0.084), -20: np.float64(0.093), -15: np.float64(0.103), -10: np.float64(0.106), -5: np.float64(0.116), 0: np.float64(0.13), 5: np.float64(0.124), 10: np.float64(0.112), 15: np.float64(0.115), 20: np.float64(0.102), 25: np.float64(0.104), 30: np.float64(0.104)}
low contrast, phase-cue coherence by lag (frames): {-30: np.float64(0.375), -25: np.float64(0.38), -20: np.float64(0.396), -15: np.float64(0.42), -10: np.float64(0.436), -5: np.float64(0.463), 0: np.float64(0.533), 5: np.float64(0.532), 10: np.float64(0.483), 15: np.float64(0.457), 20: np.float64(0.442), 25: np.float64(0.428), 30: np.float64(0.403)}
high contrast, phase-cue coherence by lag (frames): {-30: np.float64(0.513), -25: np.float64(0.538), -20: np.float64(0.568), -15: np.float64(0.602), -10: np.float64(0.658), -5: np.float64(0.714), 0: np.float64(0.79), 5: np.float64(0.762), 10: np.float64(0.694), 15: np.float64(0.659), 20: np.float64(0.621), 25: np.float64(0.577), 30: np.float64(0.542)}
lag of peak coherence in high contrast, per fly (frames): [np.int64(0), np.int64(5), np.int64(0), np.int64(0), np.int64(0), np.int64(5), np.int64(5), np.int64(0), np.int64(0), np.int64(5), np.int64(0), np.int64(5), np.int64(5), np.int64(5), np.int64(0)]
median coherence in frames labelled high contrast: true labels 0.790, labels shifted one block 0.390
flies where the true labels give the higher coherence: 12 of 15
synthetic bump at +1.000 rad, amplitude 3: phase -1.000, amplitude 3.000 (the phase sign is a convention, see bra.fourier)
ROIs per fly (frequency 2 assumes the ROIs span two bump periods): [37, 38, 39, 41, 42, 43, 47, 41, 44, 44, 42, 43, 42, 44, 41]
The bump phase locks to the cue position at high contrast (median coherence 0.79) and barely in darkness (0.13), and in every fly the coupling peaks at lag 0 or 5 frames (5 frames is 0.54 s at 9.18 Hz and 0.64 s at 7.87 Hz), the lags closest to 0 that were tested, so the dF/F and the cue position share a clock. Frames the trials table labels high contrast show stronger coupling than the same frames with the labels shifted by one block, (median 0.79 against 0.39, 12 of 15 flies), which is consistent with the trials table being on the same clock. This is consistent with the data, not a proof.
fig1-conflict-decreases-certainty¶
In the circular Kalman filter, an HD observation consistent with the current estimate increases certainty, while one that strongly conflicts with it (opposite direction) decreases certainty.
The code's circKF step, vM_Projection (kappa_phi = 1, kappa_v = 1, dt = 0.01), from mu = 0 at
kappa_t = 1, 2, 5, with the per-step reliability xi_fun_inv(gamma_z dt) for gamma_z = 1/s and 10/s.
(a) is the vector addition alone (the code's update on the natural parameters); (b) is the full
step against the same step with no observation.
rows = []
for gz in (1, 10):
alpha = flt.xi_fun_inv(gz * 0.01)
for k in (1, 2, 5):
upd_conf = abs(k * np.exp(0j) + alpha * np.exp(0j)) # observation at z = mu
upd_opp = abs(k * np.exp(0j) + alpha * np.exp(1j * np.pi)) # observation at z = mu + pi
_, full_conf = flt.vM_Projection(0., k, 1., z=0., kappa_z=alpha, dy=0., kappa_v=1., dt=0.01)
_, full_opp = flt.vM_Projection(0., k, 1., z=np.pi, kappa_z=alpha, dy=0., kappa_v=1., dt=0.01)
_, full_none = flt.vM_Projection(0., k, 1., z=0., kappa_z=0, dy=0., kappa_v=1., dt=0.01)
rows.append(dict(gamma_z=gz, alpha=round(alpha, 4), kappa_t=k, a_conf=upd_conf, a_opp=upd_opp,
b_conf=full_conf, b_none=full_none, b_opp=full_opp,
a_holds=upd_conf > k > upd_opp, b_holds=full_conf > full_none > full_opp))
fig1 = pd.DataFrame(rows)
print(fig1.round(4).to_string(index=False))
print("\n(a) holds in %d of 6 cases, (b) in %d of 6" % (fig1.a_holds.sum(), fig1.b_holds.sum()))
gamma_z alpha kappa_t a_conf a_opp b_conf b_none b_opp a_holds b_holds
1 0.1416 1 1.1416 0.8584 1.1378 0.9969 0.8558 True True
1 0.1416 2 2.1416 1.8584 2.1293 1.9894 1.8493 True True
1 0.1416 5 5.1416 4.8584 5.0388 4.9037 4.7684 True True
10 0.4529 1 1.4529 0.5471 1.4472 0.9969 0.5457 True True
10 0.4529 2 2.4529 1.5471 2.4361 1.9894 1.5408 True True
10 0.4529 5 5.4529 4.5471 5.3349 4.9037 4.4703 True True
(a) holds in 6 of 6 cases, (b) in 6 of 6
Result: reproduced. In all six cases (three certainties, two reliabilities) a confirmatory observation raises the certainty and an opposite one lowers it, both for the observation update alone and for the full step against the step without observation. Comparison (a) is a weak test: adding an observation of weight alpha at the estimate or opposite to it gives kappa + alpha and |kappa - alpha|, so it holds for any 0 < alpha < 2 kappa; (b) is the informative comparison.
fig2h-bump-tracks-circkf¶
In a simulated Bayesian ring attractor of 80 neurons, the bump position tracks HD about as accurately as the circular Kalman filter, the bump amplitude tracks the filter's certainty, and when the visual cue is removed the amplitude decays toward its baseline.
The simulation cells of Figure2h.ipynb are executed as released (seed 42, T = 15, N = 80, cue off
from 5 to 10 s, rate 15) and with the Methods' rate of 10/s. (a) mean absolute circular difference
network to circKF against the circKF's own mean absolute error; (b) Pearson correlation of bump
amplitude and circKF certainty over the cue periods; decay over the darkness.
f2h = {}
for rate in (15, 10):
ns = bra.run_figure2h(code, rate)
t, phi = ns["t"], ns["phi"]
muK, kK, muN, kN = ns["mu_ProF"], ns["kappa_ProF"], ns["mu_NN"], ns["kappa_NN"]
cue = (t < 5) | (t >= 10)
i0, i1 = ns["ind_off"], ns["ind_on"] - 1
f2h[rate] = dict(
net_vs_kf=np.mean(np.abs(bra.circ_diff(muN, muK))), kf_err=np.mean(np.abs(bra.circ_diff(muK, phi))),
net_err=np.mean(np.abs(bra.circ_diff(muN, phi))), r_cue=np.corrcoef(kN[cue], kK[cue])[0, 1],
amp_dark_start=kN[i0], amp_dark_end=kN[i1], kf_dark_start=kK[i0], kf_dark_end=kK[i1],
amp_range=f"{kN.min():.2f} to {kN.max():.2f}")
d = f2h[rate]
d["a"] = d["net_vs_kf"] < d["kf_err"]
d["b"] = d["r_cue"] >= 0.9
d["decay"] = (d["amp_dark_end"] < d["amp_dark_start"]) and abs(d["amp_dark_end"] - 1) < abs(d["amp_dark_start"] - 1)
print(pd.DataFrame(f2h).T.to_string())
net_vs_kf kf_err net_err r_cue amp_dark_start amp_dark_end kf_dark_start kf_dark_end amp_range a b decay 15 0.016266 0.483045 0.487024 0.993624 7.214148 1.220141 7.270836 0.881076 1.22 to 10.37 True True True 10 0.02131 0.490626 0.488095 0.99313 5.8882 1.211796 5.955798 0.86244 0.55 to 10.00 True True True
Result: reproduced. Under both cue strengths the network's position differs from the circKF's by a mean of 0.016 to 0.021 rad, against a circKF error to the true HD of 0.48 to 0.49 rad. Bump amplitude and circKF certainty correlate at r = 0.99 over the cue periods. In darkness the amplitude falls from 7.2 (rate 15) or 5.9 (rate 10) to 1.22 or 1.21, near the baseline of 1.
fig2i-accuracy-and-certainty-vs-reliability¶
The Bayesian ring attractor tracks HD with the same accuracy as the circular Kalman filter when HD observations are reliable and with slightly lower accuracy when they are unreliable; in the unreliable regime its certainty stays near its baseline and overestimates the filter's certainty.
accK, accB, kK, kB = vm[0, :, 1], vq[0, :, 1], vm[0, :, 2], vq[0, :, 2]
print(pd.DataFrame(dict(rate=R, acc_circKF=accK, acc_attractor=accB, diff=accB - accK,
cert_circKF=kK, cert_attractor=kB)).round(4).to_string(index=False))
res = {}
for lo in (1, 3):
for bound in (0.01, 0.02):
res[(lo, bound)] = np.abs(accB - accK)[R >= lo].max() <= bound
print("reliable = rates >= %g/s: max |diff| = %.4f; within %.2f: %s" %
(lo, np.abs(accB - accK)[R >= lo].max(), bound, res[(lo, bound)]))
low = R <= 0.1
print("rates 0.01 to 0.1: attractor accuracy below circKF at all:", bool(np.all(accB[low] < accK[low])),
(accB - accK)[low].round(4))
print("rates 0.01 to 0.1: attractor certainty above circKF at all:", bool(np.all(kB[low] > kK[low])),
(kB - kK)[low].round(3))
rate acc_circKF acc_attractor diff cert_circKF cert_attractor 0.01 0.1020 0.0292 -0.0729 0.1763 1.0006 0.03 0.1529 0.0694 -0.0835 0.3076 1.0013 0.10 0.2732 0.2159 -0.0573 0.5602 1.0525 0.20 0.3552 0.3214 -0.0337 0.7931 1.1625 0.30 0.4217 0.3966 -0.0251 0.9756 1.2657 0.40 0.4677 0.4588 -0.0089 1.1226 1.3594 0.50 0.5063 0.4947 -0.0116 1.2674 1.4644 0.60 0.5359 0.5286 -0.0073 1.3900 1.5546 0.70 0.5680 0.5657 -0.0023 1.4781 1.6158 0.80 0.5809 0.5796 -0.0013 1.5643 1.6804 0.90 0.6094 0.6081 -0.0012 1.6750 1.7668 1.00 0.6308 0.6297 -0.0010 1.7613 1.8350 2.00 0.7272 0.7277 0.0005 2.4186 2.3919 3.00 0.7785 0.7786 0.0002 2.8609 2.7976 4.00 0.8142 0.8143 0.0001 3.2569 3.1772 5.00 0.8321 0.8318 -0.0003 3.5707 3.4859 6.00 0.8470 0.8468 -0.0002 3.8873 3.8012 7.00 0.8606 0.8606 -0.0000 4.1506 4.0659 8.00 0.8680 0.8677 -0.0002 4.3998 4.3180 9.00 0.8802 0.8802 0.0000 4.6229 4.5446 10.00 0.8845 0.8844 -0.0000 4.8197 4.7451 30.00 0.9366 0.9366 -0.0000 7.8513 7.8164 100.00 0.9654 0.9654 -0.0000 13.5953 13.5798 reliable = rates >= 1/s: max |diff| = 0.0010; within 0.01: True reliable = rates >= 1/s: max |diff| = 0.0010; within 0.02: True reliable = rates >= 3/s: max |diff| = 0.0003; within 0.01: True reliable = rates >= 3/s: max |diff| = 0.0003; within 0.02: True rates 0.01 to 0.1: attractor accuracy below circKF at all: True [-0.0729 -0.0835 -0.0573] rates 0.01 to 0.1: attractor certainty above circKF at all: True [0.824 0.694 0.492]
Result: reproduced. At rates of 1/s and above the accuracies differ by at most 0.0010, so all four combinations hold. At 0.01 to 0.1/s the attractor is less accurate (by 0.057 to 0.084) and its mean certainty is higher than the circKF's (by 0.49 to 0.82).
fig3b-bayesian-beats-conventional¶
A conventional ring attractor (fast decay to a fixed bump amplitude, with that amplitude optimized) is less accurate than the Bayesian ring attractor at every HD observation reliability, and reaches its accuracy only when observations are highly unreliable.
score = np.array([bra.objective(nu[i, :, 1]) for i in range(len(bra.KAPPA_FIXED))])
ib = int(np.argmax(score))
print("weighted accuracy of the fixed-certainty filter by kappa*:",
dict(zip(bra.KAPPA_FIXED[:15].tolist(), np.round(score[:15], 4))), "...")
print("optimal kappa* = %g (weighted accuracy %.4f; Bayesian attractor %.4f)" %
(bra.KAPPA_FIXED[ib], score[ib], bra.objective(accB)))
dB = accB - nu[ib, :, 1]
print(pd.DataFrame(dict(rate=R, bayes=accB, conventional=nu[ib, :, 1], diff=dB)).round(4).to_string(index=False))
mid = (R >= 0.1) & (R <= 10)
print("Bayesian >= conventional at all 23 rates:", bool(np.all(dB >= 0)), "; min diff %.4f" % dB.min())
print("strictly above at every rate 0.1 to 10:", bool(np.all(dB[mid] > 0)), "; min there %.4f" % dB[mid].min())
print("three smallest differences at rates:", R[np.argsort(dB)[:3]], np.sort(dB)[:3].round(4))
weighted accuracy of the fixed-certainty filter by kappa*: {0.1: np.float64(0.1329), 0.3: np.float64(0.2732), 0.9: np.float64(0.5616), 1.0: np.float64(0.5771), 1.1: np.float64(0.5943), 1.2: np.float64(0.6043), 1.3: np.float64(0.6126), 1.4: np.float64(0.6201), 1.5: np.float64(0.6229), 1.6: np.float64(0.6246), 1.7: np.float64(0.6257), 1.8: np.float64(0.6253), 1.9: np.float64(0.6247), 2.0: np.float64(0.6248), 2.1: np.float64(0.6218)} ...
optimal kappa* = 1.7 (weighted accuracy 0.6257; Bayesian attractor 0.6871)
rate bayes conventional diff
0.01 0.0292 0.0177 0.0115
0.03 0.0694 0.0336 0.0358
0.10 0.2159 0.1252 0.0907
0.20 0.3214 0.1988 0.1226
0.30 0.3966 0.2786 0.1179
0.40 0.4588 0.3415 0.1173
0.50 0.4947 0.4002 0.0945
0.60 0.5286 0.4499 0.0787
0.70 0.5657 0.4862 0.0795
0.80 0.5796 0.5128 0.0668
0.90 0.6081 0.5473 0.0609
1.00 0.6297 0.5674 0.0623
2.00 0.7277 0.6934 0.0343
3.00 0.7786 0.7437 0.0349
4.00 0.8143 0.7662 0.0481
5.00 0.8318 0.7804 0.0515
6.00 0.8468 0.7905 0.0563
7.00 0.8606 0.7971 0.0635
8.00 0.8677 0.8001 0.0676
9.00 0.8802 0.8127 0.0675
10.00 0.8844 0.8083 0.0762
30.00 0.9366 0.8275 0.1091
100.00 0.9654 0.8344 0.1310
Bayesian >= conventional at all 23 rates: True ; min diff 0.0115
strictly above at every rate 0.1 to 10: True ; min there 0.0343
three smallest differences at rates: [0.01 2. 3. ] [0.0115 0.0343 0.0349]
Result: reproduced. With the optimal kappa* = 1.7, the Bayesian attractor is more accurate than the conventional one at all 23 rates (smallest difference 0.0115, at 0.01/s), and strictly so from 0.1 to 10/s. The reported, not judged, part: the difference is smallest at 0.01/s, but at 0.03/s (0.036) it is larger than at 2 and 3/s (0.034 and 0.035), so "only once highly unreliable" holds at the lowest rate only.
fig3c-update-weights¶
The weight given to a new HD observation by the Bayesian ring attractor depends on both the observation's reliability and the current certainty and is close to that of the circular Kalman filter, deviating only for very uninformative observations; a conventional attractor weights observations by their reliability alone.
weights = {"i": (0.01, lambda a, k: np.arctan(a / k) / np.pi),
"ii": (0.1, lambda a, k: np.arctan2(a, k) / (np.pi / 2))}
ok = {}
for name, (dt, w) in weights.items():
alpha = np.array([flt.xi_fun_inv(r * dt) for r in R])
wK, wB = w(alpha, kK), w(alpha, kB)
d = np.abs(wK - wB)
print("reading (%s), dt = %g: weights circKF %.4f to %.4f, attractor %.4f to %.4f" % (name, dt, wK.min(), wK.max(), wB.min(), wB.max()))
print(" |diff| by rate:", dict(zip(R.tolist(), d.round(4))))
for bound in (0.02, 0.05):
ok[(name, bound)] = d[R >= 1].max() <= bound
print(" max |diff| at rates >= 1: %.4f, within %.2f: %s" % (d[R >= 1].max(), bound, ok[(name, bound)]))
print(" largest |diff| %.4f at rate %g (below 1/s: %s)" % (d.max(), R[d.argmax()], R[d.argmax()] < 1))
reading (i), dt = 0.01: weights circKF 0.0251 to 0.0375, attractor 0.0045 to 0.0375
|diff| by rate: {0.01: np.float64(0.021), 0.03: np.float64(0.0175), 0.1: np.float64(0.0118), 0.2: np.float64(0.008), 0.3: np.float64(0.0058), 0.4: np.float64(0.0044), 0.5: np.float64(0.0034), 0.6: np.float64(0.0026), 0.7: np.float64(0.0022), 0.8: np.float64(0.0018), 0.9: np.float64(0.0013), 1.0: np.float64(0.001), 2.0: np.float64(0.0003), 3.0: np.float64(0.0006), 4.0: np.float64(0.0007), 5.0: np.float64(0.0007), 6.0: np.float64(0.0006), 7.0: np.float64(0.0006), 8.0: np.float64(0.0005), 9.0: np.float64(0.0005), 10.0: np.float64(0.0005), 30.0: np.float64(0.0001), 100.0: np.float64(0.0)}
max |diff| at rates >= 1: 0.0010, within 0.02: True
max |diff| at rates >= 1: 0.0010, within 0.05: True
largest |diff| 0.0210 at rate 0.01 (below 1/s: True)
reading (ii), dt = 0.1: weights circKF 0.1566 to 0.4191, attractor 0.0284 to 0.4194
|diff| by rate: {0.01: np.float64(0.1297), 0.03: np.float64(0.108), 0.1: np.float64(0.0725), 0.2: np.float64(0.0489), 0.3: np.float64(0.035), 0.4: np.float64(0.0266), 0.5: np.float64(0.0204), 0.6: np.float64(0.016), 0.7: np.float64(0.0131), 0.8: np.float64(0.0107), 0.9: np.float64(0.008), 1.0: np.float64(0.0062), 2.0: np.float64(0.0018), 3.0: np.float64(0.0037), 4.0: np.float64(0.0042), 5.0: np.float64(0.0042), 6.0: np.float64(0.004), 7.0: np.float64(0.0038), 8.0: np.float64(0.0035), 9.0: np.float64(0.0032), 10.0: np.float64(0.003), 30.0: np.float64(0.0011), 100.0: np.float64(0.0004)}
max |diff| at rates >= 1: 0.0062, within 0.02: True
max |diff| at rates >= 1: 0.0062, within 0.05: True
largest |diff| 0.1297 at rate 0.01 (below 1/s: True)
Result: reproduced. At rates of 1/s and above the weights under the two models' certainties differ by at most 0.0010 (Methods' weight) and 0.0062 (notebook's weight), inside both bounds, and the largest difference (0.021 and 0.130) falls at the lowest rate, 0.01/s.
fig3d-best-network-slow-decay¶
Across ring attractors with different bump decay speeds and baseline amplitudes, the best network has a slow decay (low beta) and performs about as well as the Bayesian ring attractor; performance drops as decay speed increases, which a larger baseline amplitude partly offsets.
m3d = code.npz("figure3d.npz")["MSE_array"] / ref_pf
i, j = np.unravel_index(np.nanargmax(m3d), m3d.shape)
kb, bb = bra.KSTAR_3D[i], bra.BETA_GRID[j]
bay = m3d[list(bra.KSTAR_3D).index(1.0), list(bra.BETA_GRID).index(0.5)]
print("normalized performance range %.4f to %.4f, NaN cells %d" % (np.nanmin(m3d), np.nanmax(m3d), np.isnan(m3d).sum()))
print("best network: kappa* = %g, beta = %g, performance %.4f; Bayesian (1, 0.5): %.4f, %.2f%% below the best"
% (kb, bb, m3d[i, j], bay, 100 * (m3d[i, j] - bay) / m3d[i, j]))
above = bra.BETA_GRID > bb
best_per_beta = np.nanmax(m3d, 0)
kstar_per_beta = bra.KSTAR_3D[np.nanargmax(m3d, 0)]
print(pd.DataFrame(dict(beta=bra.BETA_GRID, max_over_kstar=best_per_beta, best_kstar=kstar_per_beta)).round(4).T.to_string())
steps = np.diff(best_per_beta[j:])
c = {1: bb <= 1, 2: (kb, bb) != (1.0, 0.5), 3: (m3d[i, j] - bay) / m3d[i, j] <= 0.10,
4: bool(np.all(steps < 0)), 5: bool(np.all(kstar_per_beta[above] >= kb))}
print("conditions:", c)
print("steps of the maximum for beta above the best:", steps.round(4), "; rises:", int((steps > 0).sum()), "of", len(steps))
rho = stats.spearmanr(bra.BETA_GRID[j:], best_per_beta[j:])
print("variant, not judged: Spearman rho of beta and the maximum, beta >= best: %.3f (p = %.2g); drop from beta %g to 8: %.4f"
% (rho.statistic, rho.pvalue, bb, best_per_beta[j] - best_per_beta[-1]))
# Monte Carlo error of one grid cell: each cell is its own scan of 10,000 runs per rate, a separate
# process of performance_scan_figure3d.py started with (kappa*, beta) as arguments. The script sets
# no seed, so numpy's global generator is seeded from the operating system in each process and the
# cells' trajectories are independent. The objective is linear in the per-rate accuracies.
print("performance_scan_figure3d.py, lines that set the grid cell or a seed:",
code.source_lines("performance_scan_figure3d.py", r"sys\.argv|seed|nIter ="))
# The per-rate SE is taken at the Bayesian attractor's accuracies (figure2i.npz), close to the networks'.
se_cell = bra.objective_se(bra.resultant_se(accB)[0]) / ref_pf
print("SE of one cell's normalized performance %.4f (distribution-free bound %.4f); SE of a difference of two cells %.4f"
% (se_cell, bra.objective_se(bra.resultant_se(accB)[1]) / ref_pf, np.sqrt(2) * se_cell))
print("rises of the maximum over kappa*, in SEs of a difference:", (steps[steps > 0] / (np.sqrt(2) * se_cell)).round(2))
normalized performance range 0.6520 to 1.0012, NaN cells 15
best network: kappa* = 0.7, beta = 0.3, performance 1.0012; Bayesian (1, 0.5): 0.9955, 0.57% below the best
0 1 2 3 4 5 6 7 8 9 10 11 12
beta 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.900 1.0000 2.0000 4.0000 8.0000
max_over_kstar 0.9992 0.9986 1.0012 0.9981 0.9968 0.9947 0.9934 0.9942 0.988 0.9893 0.9725 0.9525 0.9322
best_kstar 0.3000 0.5000 0.7000 0.8000 0.8000 0.9000 1.1000 1.1000 1.100 1.0000 1.2000 1.6000 1.5000
conditions: {1: np.True_, 2: True, 3: np.True_, 4: False, 5: True}
steps of the maximum for beta above the best: [-0.0031 -0.0013 -0.0021 -0.0013 0.0008 -0.0062 0.0014 -0.0168 -0.02 -0.0203] ; rises: 2 of 10
variant, not judged: Spearman rho of beta and the maximum, beta >= best: -0.982 (p = 8.4e-08); drop from beta 0.3 to 8: 0.0689
performance_scan_figure3d.py, lines that set the grid cell or a seed: ['kappa_star = float(sys.argv[1])', 'beta = float(sys.argv[2])', 'nIter = 10000']
SE of one cell's normalized performance 0.0018 (distribution-free bound 0.0032); SE of a difference of two cells 0.0026
rises of the maximum over kappa*, in SEs of a difference: [0.32 0.52]
Result: partially reproduced. The best network (kappa* = 0.7, beta = 0.3) decays slowly, differs from the Bayesian attractor's parameters, and the Bayesian attractor is 0.57% below it. Above the best beta the kappa* that maximizes performance is always at least 0.7, rising to 1.5 or 1.6. The maximum over kappa* falls overall from 1.0012 to 0.9322 (Spearman rho printed above), but not at every step: it rises twice, by 0.0008 (beta 0.7 to 0.8) and 0.0014 (0.9 to 1), so condition (4) fails as written. The two rises are 0.32 and 0.52 standard errors of a difference between two grid cells (0.0026; the scan script sets no seed and runs each cell as its own process, so the cells' trajectories are independent), so they are within the scan's sampling error; the criterion gives no tolerance, so the outcome stays partial.
fig3f-noise-drop-below-70¶
In a network of 64 neurons with additive neural noise, inference accuracy drops below 70% of its noise-free value once the noise is large enough.
grids = {s: code.npz(f"figure3f/sigma={s}.npz") for s in bra.SIGMAS_3F}
M0 = grids[0.0]["MSE_array"] / bra.REF_PF_3F
idx_ref = np.unravel_index(np.nanargmax(M0), M0.shape)
print("noise-free optimum: kappa* = %g, beta = %g" % (bra.KSTAR_3F[idx_ref[0]], bra.BETA_GRID[idx_ref[1]]))
rows = []
for s in bra.SIGMAS_3F:
M, act = grids[s]["MSE_array"] / bra.REF_PF_3F, grids[s]["avg_network_activity"]
idx = np.unravel_index(np.nanargmax(M), M.shape)
rows.append(dict(sigma=s, retuned_pf=100 * M[idx], untuned_pf=100 * M[idx_ref],
retuned_self=100 * M[idx] / M0[idx_ref], untuned_self=100 * M[idx_ref] / M0[idx_ref],
snr_retuned=act[idx] / s * np.sqrt(bra.N_3F / 2) if s else np.nan,
snr_untuned=act[idx_ref] / s * np.sqrt(bra.N_3F / 2) if s else np.nan,
best_kstar=bra.KSTAR_3F[idx[0]], best_beta=bra.BETA_GRID[idx[1]], nan_cells=int(np.isnan(M).sum())))
f3f = pd.DataFrame(rows).set_index("sigma")
print(f3f.round(3).to_string())
for norm_ in ("pf", "self"):
for smax in (10, 20):
v = f3f.loc[float(smax), "untuned_" + norm_]
print("normalization %-4s, largest SD %d: untuned performance %.2f%% -> below 70%%: %s" % (norm_, smax, v, v < 70))
noise-free optimum: kappa* = 0.5, beta = 0.2
retuned_pf untuned_pf retuned_self untuned_self snr_retuned snr_untuned best_kstar best_beta nan_cells
sigma
0.0 100.081 100.081 100.000 100.000 NaN NaN 0.5 0.2 15
1.0 99.445 99.155 99.364 99.074 14.497 13.860 0.3 0.1 15
2.0 97.784 97.349 97.704 97.270 7.849 6.758 0.9 0.3 15
3.0 95.378 93.362 95.300 93.286 5.352 4.356 1.2 0.4 15
4.0 92.669 89.050 92.593 88.978 4.423 3.152 1.7 0.5 15
5.0 89.653 84.511 89.580 84.443 3.620 2.450 2.2 0.7 15
6.0 86.513 79.834 86.443 79.769 3.036 1.997 2.7 1.0 15
7.0 83.686 75.333 83.618 75.272 3.081 1.678 3.2 0.8 15
8.0 81.073 71.199 81.007 71.141 2.719 1.447 3.2 0.7 15
9.0 78.335 67.440 78.271 67.385 2.658 1.271 4.2 1.0 15
10.0 75.897 63.123 75.835 63.072 2.325 1.136 4.8 2.0 15
15.0 65.030 48.385 64.977 48.345 2.053 0.773 4.6 0.4 15
20.0 56.541 38.094 56.495 38.063 1.850 0.616 2.4 0.1 15
normalization pf , largest SD 10: untuned performance 63.12% -> below 70%: True
normalization pf , largest SD 20: untuned performance 38.09% -> below 70%: True
normalization self, largest SD 10: untuned performance 63.07% -> below 70%: True
normalization self, largest SD 20: untuned performance 38.06% -> below 70%: True
Result: reproduced. Without retuning, the network's relative performance at noise SD 10 is 63.1% under both normalizations and at SD 20 it is 38.1%, below 70% under all four combinations.
fig3f-retuning-reduces-drop¶
Retuning the bump decay speed and baseline amplitude for each noise level reduces the drop in inference accuracy by about 35% compared with the network tuned for noise-free performance, and raises the signal-to-noise ratio.
Expected: reduction of the drop 35 +/- 3.5%.
red = {}
for norm_ in ("pf", "self"):
drop_r = 100 - f3f["retuned_" + norm_]
drop_u = 100 - f3f["untuned_" + norm_]
r_ = 100 * (1 - drop_r / drop_u)
red[(norm_, "a: SD 10")] = r_.loc[10.0]
red[(norm_, "b: SD 20")] = r_.loc[20.0]
red[(norm_, "c: mean SD 1-10")] = r_.loc[1.0:10.0].mean()
print("normalization", norm_, "reduction by SD:", r_.round(2).to_dict())
red = pd.Series(red)
print("\nreadings of the 35%:\n", red.round(2).to_string())
print("readings within 31.5 to 38.5%:", list(red[(red >= 31.5) & (red <= 38.5)].index))
snr_ok = bool(np.all(f3f.snr_retuned.iloc[1:] >= f3f.snr_untuned.iloc[1:]))
print("retuned S/N >= untuned at every SD > 0:", snr_ok, "; ratio range %.2f to %.2f"
% ((f3f.snr_retuned / f3f.snr_untuned).iloc[1:].min(), (f3f.snr_retuned / f3f.snr_untuned).iloc[1:].max()))
print("Fig. S4, reported only: best kappa* by SD", f3f.best_kstar.to_dict(), "\n best beta by SD", f3f.best_beta.to_dict())
normalization pf reduction by SD: {0.0: 0.0, 1.0: 34.35, 2.0: 16.4, 3.0: 30.37, 4.0: 33.05, 5.0: 33.2, 6.0: 33.12, 7.0: 33.86, 8.0: 34.28, 9.0: 33.46, 10.0: 34.64, 15.0: 32.25, 20.0: 29.8}
normalization self reduction by SD: {0.0: nan, 1.0: 31.33, 2.0: 15.91, 3.0: 30.0, 4.0: 32.8, 5.0: 33.02, 6.0: 32.99, 7.0: 33.75, 8.0: 34.19, 9.0: 33.38, 10.0: 34.56, 15.0: 32.2, 20.0: 29.76}
readings of the 35%:
pf a: SD 10 34.64
b: SD 20 29.80
c: mean SD 1-10 31.67
self a: SD 10 34.56
b: SD 20 29.76
c: mean SD 1-10 31.19
readings within 31.5 to 38.5%: [('pf', 'a: SD 10'), ('pf', 'c: mean SD 1-10'), ('self', 'a: SD 10')]
retuned S/N >= untuned at every SD > 0: True ; ratio range 1.05 to 3.00
Fig. S4, reported only: best kappa* by SD {0.0: 0.5, 1.0: 0.3, 2.0: 0.9, 3.0: 1.2, 4.0: 1.7, 5.0: 2.2, 6.0: 2.7, 7.0: 3.2, 8.0: 3.2, 9.0: 4.2, 10.0: 4.8, 15.0: 4.6, 20.0: 2.4}
best beta by SD {0.0: 0.2, 1.0: 0.1, 2.0: 0.3, 3.0: 0.4, 4.0: 0.5, 5.0: 0.7, 6.0: 1.0, 7.0: 0.8, 8.0: 0.7, 9.0: 1.0, 10.0: 2.0, 15.0: 0.4, 20.0: 0.1}
Result: reproduced. At noise SD 10 retuning reduces the drop by 34.6% under both normalizations, and the mean over SD 1 to 10 is 31.7% with the particle filter normalization, all inside 31.5 to 38.5%; the other readings give 31.2% (mean, legend normalization) and 29.8% (SD 20), outside it. The retuned signal-to-noise ratio is above the untuned one at every noise level (ratio 1.05 to 3.00). For Fig. S4, reported only: the best kappa* grows with the noise from SD 1 to 10 (0.3 to 4.8), while the best beta does not grow steadily (0.1 to 2.0, then 0.4 and 0.1 at SD 15 and 20).
noise-effects-on-mu-and-kappa¶
Neural noise makes the bump position diffuse without bias and makes the bump amplitude drift upward and diffuse, which leads to overconfidence; for additive noise, a four-fold larger network halves the noise's effective SD.
RNN_filter_Run at kappa* = 1, beta = 0.5, tau = 1 (w_even = 1.5, w_quad = 0.5), no inputs,
stoch_corr = 0, from a bump of amplitude 1, 2,000 runs of T = 5, noise SD 1 and 2, N = 64 and 256.
The runs use a batched copy of the function (bra.rnn_batch), first checked against the released
function on the same noise.
# the batched copy against the released RNN_filter_Run, same noise draws
np.random.seed(3)
_, mu_ref, k_ref = flt.RNN_filter_Run(0.5, N=16, dt=0.01, w_even=1.5, w_quad=0.5, tau=1, I_ext=0,
stoch_corr=0, phi_0=0, kappa_0=1, sigma_N=1.0)
np.random.seed(3)
noise = np.random.randn(50, 16)[:, None, :]
mu_b, k_b = bra.rnn_batch(None, 1, 0.5, 16, 1.0, 1.5, 0.5, noise=noise)
print("batched vs released, max abs difference: position %.2e, amplitude %.2e"
% (np.abs(mu_b[:, 0] - mu_ref).max(), np.abs(k_b[:, 0] - k_ref).max()))
rng = np.random.default_rng(289)
noise_res, kmean_t = {}, {}
for N in (64, 256):
_, k0 = bra.rnn_batch(rng, 1, 5, N, 0.0, 1.5, 0.5)
for sig in (1.0, 2.0):
mu, k = bra.rnn_batch(rng, 2000, 5, N, sig, 1.5, 0.5)
dmu = bra.circ_diff(mu[-1], mu[0])
var_t = [np.var(bra.circ_diff(mu[int(t / 0.01) - 1], mu[0])) for t in (1, 2, 3, 4, 5)]
kmean_t[(N, sig)] = [k[int(t / 0.01) - 1].mean() for t in (0.1, 0.5, 1, 2, 3, 4, 5)]
noise_res[(N, sig)] = dict(mean_dmu=dmu.mean(), se=dmu.std(ddof=1) / np.sqrt(len(dmu)),
z=dmu.mean() / (dmu.std(ddof=1) / np.sqrt(len(dmu))),
var_t1_to_t5=np.round(var_t, 4), var_grows=bool(np.all(np.diff(var_t) > 0)),
amp_end=k[-1].mean(), amp_end_sd=k[-1].std(), amp_noise_free=k0[-1, 0],
inc_sd=np.diff(k, axis=0).std())
nr = pd.DataFrame(noise_res).T
print(nr.to_string())
print("\nmean amplitude at t = 0.1, 0.5, 1, 2, 3, 4, 5 s:")
for key, v in kmean_t.items():
print(key, np.round(v, 3))
ratio = {sig: noise_res[(256, sig)]["inc_sd"] / noise_res[(64, sig)]["inc_sd"] for sig in (1.0, 2.0)}
c1 = all(abs(v["z"]) <= 2 and v["var_grows"] for v in noise_res.values())
c2 = all(v["amp_end"] > v["amp_noise_free"] for v in noise_res.values())
c3 = all(0.45 <= r_ <= 0.55 for r_ in ratio.values())
print("\n(1) unbiased and growing variance in all four:", c1, [round(v["z"], 2) for v in noise_res.values()])
print("(2) mean end amplitude above noise-free in all four:", c2)
print("(3) increment SD ratio N=256/N=64:", {k_: round(v, 4) for k_, v in ratio.items()}, c3)
# Diagnostics (the verifier's request). (i) In the full network, the run-mean global inhibition
# w_quad pi/N sum max(r, 0) and the zero mode (mean rate) over time. (ii) The same runs with the
# global inhibition held at its noise-free level, so the noise acts only as one more input to
# a constant. This changes two things at once: it removes the noise's effect on the inhibition, and
# the inhibition's dependence on the amplitude, the restoring force toward kappa* = 1, so the first
# mode is neutral (w_even - 1/tau - inhibition = 0). SI Eq. S124 gives each coordinate of (x1, x2)
# the variance 2 sigma^2 / N per unit time, so the amplitude after T is Rice distributed.
# (iii) The verifier's cleaner control: the inhibition set to the value a noise-free cosine bump of
# the run's decoded amplitude gives (inh_per_amp times kappa), so the restoring force is kept and only
# the noise's own contribution to the rectified sum is removed. Its reference is the reduced
# two-dimensional SDE of the first mode, dx = (w_even - 1/tau - w_quad kappa) x dt + noise of
# variance 2 sigma^2 / N per coordinate and unit time, which has the SI's Ito drift sigma^2/(kappa N).
rows = []
for N in (64, 256):
tr0 = {}
_, _k0 = bra.rnn_batch(rng, 1, 5, N, 0.0, 1.5, 0.5, traces=tr0)
inh0, k_nf = tr0["inh"][-1], _k0[-1, 0]
for sig in (1.0, 2.0):
tr = {}
_, k = bra.rnn_batch(rng, 2000, 5, N, sig, 1.5, 0.5, traces=tr)
_, kf = bra.rnn_batch(rng, 2000, 5, N, sig, 1.5, 0.5, inh_fixed=inh0)
_, ka = bra.rnn_batch(rng, 2000, 5, N, sig, 1.5, 0.5, inh_per_amp=inh0 / k_nf)
s_T = np.sqrt(2 * sig ** 2 / N * 5)
x = np.tile([1.0, 0.0], (20000, 1))
for _ in range(499):
kx = np.hypot(x[:, 0], x[:, 1])[:, None]
x = x + (0.5 - 0.5 * kx) * x * 0.01 + np.sqrt(2 * sig ** 2 / N * 0.01) * rng.standard_normal(x.shape)
rows.append(dict(N=N, sigma=sig, inh_noise_free=inh0,
inh_t01=tr["inh"][9], inh_t1=tr["inh"][99], inh_t5=tr["inh"][-1],
zero_mode_t5=tr["zero_mode"][-1], amp_end_full=k[-1].mean(),
amp_end_inh_fixed=kf[-1].mean(), se_inh_fixed=kf[-1].std(ddof=1) / np.sqrt(len(kf[-1])),
amp_end_SI_rice=stats.rice(b=1 / s_T, scale=s_T).mean(),
amp_noise_free=k_nf, amp_end_inh_of_amp=ka[-1].mean(),
se_inh_of_amp=ka[-1].std(ddof=1) / np.sqrt(len(ka[-1])),
amp_end_2d_sde=np.hypot(x[:, 0], x[:, 1]).mean()))
dg = pd.DataFrame(rows)
print("\nglobal inhibition and zero mode in the full network, and the end amplitude with the inhibition held:")
print(dg.iloc[:, :11].round(4).to_string(index=False))
print("\nend amplitude with the inhibition set by the decoded amplitude (iii), against the noise-free run and the 2-D SDE:")
dg["z_vs_noise_free"] = (dg.amp_end_inh_of_amp - dg.amp_noise_free) / dg.se_inh_of_amp
dg["z_vs_2d_sde"] = (dg.amp_end_inh_of_amp - dg.amp_end_2d_sde) / dg.se_inh_of_amp
print(dg[["N", "sigma", "amp_noise_free", "amp_end_inh_of_amp", "se_inh_of_amp", "amp_end_2d_sde",
"z_vs_noise_free", "z_vs_2d_sde", "amp_end_full"]].round(4).to_string(index=False))
batched vs released, max abs difference: position 1.11e-16, amplitude 4.44e-16
mean_dmu se z var_t1_to_t5 var_grows amp_end amp_end_sd amp_noise_free inc_sd
64 1.0 -0.002646 0.010771 -0.245626 [0.0338, 0.0748, 0.1244, 0.1782, 0.2319] True 0.832435 0.209167 1.000738 0.017703
2.0 0.025019 0.030527 0.819556 [0.1765, 0.626, 1.0294, 1.47, 1.8629] True 0.646362 0.309351 1.000738 0.035337
256 1.0 -0.003728 0.005008 -0.744347 [0.0082, 0.0174, 0.0264, 0.0369, 0.0501] True 0.815035 0.115057 1.000046 0.008851
2.0 -0.008879 0.021326 -0.416355 [0.0394, 0.1221, 0.2852, 0.541, 0.9091] True 0.421835 0.192841 1.000046 0.017703
mean amplitude at t = 0.1, 0.5, 1, 2, 3, 4, 5 s:
(64, 1.0) [1. 0.986 0.955 0.909 0.873 0.851 0.832]
(64, 2.0) [1. 0.945 0.854 0.741 0.683 0.656 0.646]
(256, 1.0) [0.999 0.982 0.953 0.898 0.863 0.835 0.815]
(256, 2.0) [0.998 0.928 0.826 0.66 0.55 0.47 0.422]
(1) unbiased and growing variance in all four: True [np.float64(-0.25), np.float64(0.82), np.float64(-0.74), np.float64(-0.42)]
(2) mean end amplitude above noise-free in all four: False
(3) increment SD ratio N=256/N=64: {1.0: np.float64(0.5), 2.0: np.float64(0.501)} True
global inhibition and zero mode in the full network, and the end amplitude with the inhibition held: N sigma inh_noise_free inh_t01 inh_t1 inh_t5 zero_mode_t5 amp_end_full amp_end_inh_fixed se_inh_fixed amp_end_SI_rice 64 1.0 0.5 0.5196 0.5695 0.5281 0.0007 0.8249 1.0756 0.0085 1.0824 64 2.0 0.5 0.5889 0.7680 0.7279 -0.0096 0.6394 1.3802 0.0148 1.3522 256 1.0 0.5 0.5202 0.5662 0.5231 -0.0022 0.8214 1.0114 0.0044 1.0197 256 2.0 0.5 0.5877 0.7587 0.7071 -0.0019 0.4297 1.0820 0.0084 1.0824 end amplitude with the inhibition set by the decoded amplitude (iii), against the noise-free run and the 2-D SDE: N sigma amp_noise_free amp_end_inh_of_amp se_inh_of_amp amp_end_2d_sde z_vs_noise_free z_vs_2d_sde amp_end_full 64 1.0 1.0007 1.0021 0.0039 1.0037 0.3501 -0.3970 0.8249 64 2.0 1.0007 1.0453 0.0071 1.0373 6.2475 1.1210 0.6394 256 1.0 1.0000 1.0026 0.0020 1.0011 1.2648 0.7300 0.8214 256 2.0 1.0000 1.0046 0.0040 1.0018 1.1399 0.6857 0.4297
Result: partially reproduced. The batched network matches the released function to rounding error. The position diffuses without bias: its mean change is within 2 standard errors of 0 in all four settings (z = -0.74 to 0.82) and its variance grows at every second. The SD of the one-step amplitude increments at N = 256 is half that at N = 64 (ratios 0.500 and 0.501). But in the full network the mean amplitude does not drift upward: it falls throughout the run and ends below the noise-free value of 1.0 in all four settings (0.42 to 0.83), so condition (2) fails as the criterion runs it.
The diagnostics show why. The noise raises the rectified sum that drives the global inhibition, from 0.5 without noise to 0.52 to 0.59 at 0.1 s and 0.57 to 0.77 at 1 s; the zero mode stays near 0 (at most 0.01 in size). Holding the inhibition at its noise-free constant (ii) makes the mean amplitude rise to 1.08, 1.38, 1.01 and 1.08 after 5 s, against 1.08, 1.35, 1.02 and 1.08 from the Rice mean that SI Eq. S124 implies, each within 1.9 standard errors. That shows the SI's noise term is right, but this control changes two things at once: it removes the noise's effect on the inhibition and also the restoring force toward kappa* = 1, which leaves a free random walk of the first mode. The control (iii) removes only the first: the inhibition follows the decoded amplitude as it would for a noise-free bump. Then the mean amplitude after 5 s is 1.002, 1.045, 1.003 and 1.005, at or above the noise-free 1.000 to 1.001 in all four settings (z = 0.35, 6.2, 1.3 and 1.1) and within 1.2 standard errors of the reduced two-dimensional SDE of the first mode. So with the restoring force kept, the noise's direct effect gives a small upward drift, clearly resolved only at N = 64 and SD 2, and the fall of the amplitude in the released network comes from the noise's own contribution to the rectified sum that drives the global inhibition. The criterion tests the full network; whether that is what the paper claims is asked of a reader (questions/000289.yaml).
fig4b-connectome-motifs¶
The connectivity among EPG, Delta7 and PEN1 neurons in the hemibrain connectome resembles the connectivity of a network that implements the approximate circular Kalman filter, with shifted PEN1 to EPG projections and broad EPG to Delta7 to EPG inhibition.
Total synapse count per pair, from the public hemibrain v1.2 export (the paper used v1.2.1 in neuPrint, which now needs a personal token). Neurons are sorted by the PB glomerulus in their name, L9 to R9 (reading (a)); reading (b), by ellipsoid body or gall position in the name, cannot be formed because no name holds that position. (1) For each PEN1 neuron, the synapse-weighted circular mean position (period 8 glomeruli) of its EPG targets minus that of its EPG inputs, among EPG neurons on the same side of the bridge; the claim holds when the offsets of left and right PEN1 neurons have opposite signs. (2) The number of the 16 EPG glomeruli from which each Delta7 receives, and to which it sends, at least one synapse.
pos = {g: k for k, g in enumerate(bra.PB_ORDER)}
cx["x"] = cx.g.map(pos)
pre = cx[["bodyId", "t", "g", "x"]].add_suffix("_pre")
post = cx[["bodyId", "t", "g", "x"]].add_suffix("_post")
cc = conns.merge(pre, on="bodyId_pre").merge(post, on="bodyId_post")
def cmean(x, w):
return np.angle(np.sum(w * np.exp(2j * np.pi * x / 8)) / np.sum(w)) * 8 / (2 * np.pi)
rows = []
for b, p in cx[cx.t == "PEN1"].set_index("bodyId").iterrows():
side = p.g[0]
inp = cc[(cc.bodyId_post == b) & (cc.t_pre == "EPG")]
out = cc[(cc.bodyId_pre == b) & (cc.t_post == "EPG")]
for sel in ("same side", "both sides"):
i_ = inp[inp.g_pre.str[0] == side] if sel == "same side" else inp
o_ = out[out.g_post.str[0] == side] if sel == "same side" else out
off = (cmean(o_.x_post.values, o_.weight.values) - cmean(i_.x_pre.values, i_.weight.values) + 4) % 8 - 4
rows.append(dict(glomerulus=p.g, side=side, partners=sel, offset=off,
syn_in=i_.weight.sum(), syn_out=o_.weight.sum()))
pen = pd.DataFrame(rows)
print(pen[pen.partners == "same side"].sort_values("glomerulus", key=lambda s: s.map(pos)).round(3).to_string(index=False))
summ = pen.groupby(["partners", "side"]).offset.agg(["mean", "median", "min", "max", lambda s: int((s > 0).sum()), "size"])
print(summ.round(3).to_string())
for sel in ("same side", "both sides"):
L = pen[(pen.partners == sel) & (pen.side == "L")].offset
Rr = pen[(pen.partners == sel) & (pen.side == "R")].offset
print("%s: left mean %+.3f (Wilcoxon vs 0 p = %.3g), right mean %+.3f (p = %.3g); opposite signs: %s"
% (sel, L.mean(), stats.wilcoxon(L).pvalue, Rr.mean(), stats.wilcoxon(Rr).pvalue, np.sign(L.mean()) != np.sign(Rr.mean())))
rows = []
for b, p in cx[cx.t == "Delta7"].set_index("bodyId").iterrows():
i_ = cc[(cc.bodyId_post == b) & (cc.t_pre == "EPG")].groupby("g_pre").weight.sum()
o_ = cc[(cc.bodyId_pre == b) & (cc.t_post == "EPG")].groupby("g_post").weight.sum()
rows.append({"instance": p.instance, **{f"in>={th}": int((i_ >= th).sum()) for th in (1, 3, 5)},
**{f"out>={th}": int((o_ >= th).sum()) for th in (1, 3, 5)}, "syn_in": i_.sum(), "syn_out": o_.sum()})
d7 = pd.DataFrame(rows)
print("\nEPG glomeruli (of 16) connected to each Delta7, by minimum synapse count per glomerulus:")
print(d7.describe().loc[["min", "50%", "max"]].round(1).to_string())
for th in (1, 3, 5):
print("threshold %d synapse(s): Delta7 with more than 8 input glomeruli %d of %d, more than 8 output glomeruli %d of %d"
% (th, (d7[f"in>={th}"] > 8).sum(), len(d7), (d7[f"out>={th}"] > 8).sum(), len(d7)))
glomerulus side partners offset syn_in syn_out
L9 L same side 0.046 150 529
L8 L same side 0.520 154 406
L7 L same side 0.537 124 351
L6 L same side 0.530 108 178
L6 L same side 0.684 81 157
L5 L same side 0.598 181 328
L4 L same side 0.702 146 256
L4 L same side 0.776 157 269
L3 L same side 0.706 180 349
L2 L same side 0.543 209 378
R2 R same side -0.509 163 363
R3 R same side -0.551 203 413
R4 R same side -0.320 141 353
R4 R same side -0.453 124 298
R5 R same side -0.611 203 480
R6 R same side -0.575 143 231
R6 R same side -0.482 139 250
R7 R same side -0.690 219 363
R8 R same side -0.601 249 420
R9 R same side 0.019 111 378
mean median min max <lambda_0> size
partners side
both sides L 0.374 0.360 0.083 0.670 10 10
R -0.346 -0.351 -0.489 0.004 1 10
same side L 0.564 0.570 0.046 0.776 10 10
R -0.477 -0.530 -0.690 0.019 1 10
same side: left mean +0.564 (Wilcoxon vs 0 p = 0.00195), right mean -0.477 (p = 0.00391); opposite signs: True
both sides: left mean +0.374 (Wilcoxon vs 0 p = 0.00195), right mean -0.346 (p = 0.00391); opposite signs: True
EPG glomeruli (of 16) connected to each Delta7, by minimum synapse count per glomerulus:
in>=1 in>=3 in>=5 out>=1 out>=3 out>=5 syn_in syn_out
min 13.0 11.0 10.0 9.0 4.0 2.0 332.0 115.0
50% 16.0 14.0 13.0 12.0 7.0 5.0 404.5 155.5
max 16.0 16.0 15.0 15.0 11.0 7.0 566.0 212.0
threshold 1 synapse(s): Delta7 with more than 8 input glomeruli 42 of 42, more than 8 output glomeruli 42 of 42
threshold 3 synapse(s): Delta7 with more than 8 input glomeruli 42 of 42, more than 8 output glomeruli 7 of 42
threshold 5 synapse(s): Delta7 with more than 8 input glomeruli 42 of 42, more than 8 output glomeruli 0 of 42
Result: partially reproduced. The criterion asks for both sorts, and only reading (a), by the PB glomerulus in the name, can be formed; reading (b) cannot (data check). Under reading (a), and counting any synapse, both parts hold. Left PEN1 neurons send their EPG output a mean of +0.56 glomeruli from their EPG input (10 of 10 positive) and right ones -0.48 (9 of 10 negative), opposite directions, both different from 0 (Wilcoxon, printed above); with partners from both sides, +0.37 and -0.35. Every Delta7 receives EPG input from at least 13 of the 16 glomeruli and sends output to EPG neurons in at least 9. The output side depends on the threshold: counting only glomeruli with at least 3 or 5 synapses, Delta7 outputs reach more than half the glomeruli in 7 and 0 of 42 neurons (printed above), while the inputs still do in 42 of 42.
fig4gh-drosophila-network-matches¶
A Drosophila-like multipopulation network, constrained by the connectome, tracks HD with an accuracy indistinguishable from the idealized Bayesian ring attractor across observation reliabilities.
g4 = code.npz("figure4g.npz")
cxa, bqa = g4["CX_stats"][:, 1], g4["vonMises_q_stats"][:, 1]
d4 = np.abs(cxa - bqa)
print(pd.DataFrame(dict(rate=R, drosophila=cxa, bayes=bqa, abs_diff=d4,
cert_drosophila=g4["CX_stats"][:, 2], cert_bayes=g4["vonMises_q_stats"][:, 2])).round(4).to_string(index=False))
for bound in (0.01, 0.02):
print("max |diff| %.4f; within %.2f at every rate: %s; rates outside: %s"
% (d4.max(), bound, bool(np.all(d4 <= bound)), R[d4 > bound].tolist()))
# Monte Carlo error (10,000 runs per rate, performance_scan_figure4g.py). The two networks share
# their runs, so the SE of their difference is below that of two independent estimates, sqrt(2) SE.
se4 = bra.resultant_se(bqa)[0]
print("\nSE of the attractor's accuracy at 0.1, 0.2, 0.3/s: %s; differences in units of sqrt(2) SE: %s"
% (se4[2:5].round(4), (d4[2:5] / (np.sqrt(2) * se4[2:5])).round(2)))
print("the same Bayesian attractor in figure4g.npz (dt 0.001) and figure2i.npz (dt 0.01): largest |diff| %.4f at %g/s, %.2f sqrt(2) SE"
% (np.abs(bqa - accB).max(), R[np.argmax(np.abs(bqa - accB))], np.abs(bqa - accB).max() / (np.sqrt(2) * se4[np.argmax(np.abs(bqa - accB))])))
rate drosophila bayes abs_diff cert_drosophila cert_bayes 0.01 0.0285 0.0228 0.0057 0.8848 1.0006 0.03 0.0784 0.0719 0.0065 0.8881 1.0030 0.10 0.2133 0.1943 0.0190 0.9511 1.0568 0.20 0.3300 0.3160 0.0140 1.0634 1.1644 0.30 0.4218 0.4117 0.0101 1.1728 1.2758 0.40 0.4546 0.4513 0.0033 1.2617 1.3701 0.50 0.4958 0.4920 0.0038 1.3558 1.4692 0.60 0.5259 0.5239 0.0020 1.4229 1.5390 0.70 0.5576 0.5558 0.0018 1.4956 1.6179 0.80 0.5885 0.5878 0.0007 1.5586 1.6857 0.90 0.5962 0.5980 0.0018 1.6323 1.7637 1.00 0.6405 0.6414 0.0009 1.7006 1.8354 2.00 0.7324 0.7343 0.0018 2.2427 2.4167 3.00 0.7874 0.7892 0.0017 2.6353 2.8399 4.00 0.8112 0.8140 0.0028 2.9912 3.2202 5.00 0.8280 0.8306 0.0026 3.3128 3.5652 6.00 0.8451 0.8482 0.0031 3.5764 3.8460 7.00 0.8581 0.8598 0.0018 3.8184 4.1027 8.00 0.8670 0.8687 0.0017 4.0690 4.3704 9.00 0.8750 0.8767 0.0018 4.3071 4.6225 10.00 0.8827 0.8845 0.0018 4.5206 4.8494 30.00 0.9317 0.9337 0.0020 7.6173 8.0913 100.00 0.9617 0.9632 0.0015 13.8363 14.4548 max |diff| 0.0190; within 0.01 at every rate: False; rates outside: [0.1, 0.2, 0.3] max |diff| 0.0190; within 0.02 at every rate: True; rates outside: [] SE of the attractor's accuracy at 0.1, 0.2, 0.3/s: [0.0069 0.0065 0.0061]; differences in units of sqrt(2) SE: [1.95 1.51 1.16] the same Bayesian attractor in figure4g.npz (dt 0.001) and figure2i.npz (dt 0.01): largest |diff| 0.0216 at 0.1/s, 2.22 sqrt(2) SE
Result: partially reproduced. The Drosophila-like network's accuracy is within 0.02 of the Bayesian ring attractor's at every rate (largest 0.019, at 0.1/s), but not within 0.01: it differs by 0.019, 0.014 and 0.010 at 0.1, 0.2 and 0.3/s. Since the criterion requires both readings, the claim holds in part. The three differences above 0.01 are 1.95, 1.51 and 1.16 times sqrt(2) standard errors (the SE of a difference of two independent scans; the two networks share their runs, so the actual SE of their difference is smaller). The same Bayesian attractor differs by 0.022 between figure4g.npz (dt 0.001) and figure2i.npz (dt 0.01) at 0.1/s, 2.2 such units, so the 0.01 bound is near the scans' sampling error at the low rates.
figs2-performance-independent-of-n¶
The Bayesian ring attractor's inference accuracy depends little on its number of neurons.
s2 = code.npz("figureS2.npz")["MSE_array"]
ref2 = code.npz("vonMises_q_stats.npz")["stats"][:, 1]
print("vonMises_q_stats.npz equals the attractor's accuracy in figure2i.npz:", np.allclose(ref2, accB))
dev = pd.DataFrame(np.abs(s2 - ref2), index=bra.N_S2, columns=R)
print("missing values (NaN) per N:", dict(zip(bra.N_S2.tolist(), np.isnan(s2).sum(1).tolist())))
print("max |diff| per N:", dev.max(axis=1).round(4).to_dict())
print("cells above 0.02:", [(int(n), float(r)) for n in dev.index for r in dev.columns if dev.loc[n, r] > 0.02])
for nmin in (8, 16):
sub = dev.loc[dev.index >= nmin]
print("N >= %d: max |diff| %.4f, within 0.02 at every rate: %s" % (nmin, np.nanmax(sub.values), bool(np.nanmax(sub.values) <= 0.02)))
print("N = 3 and 4, reported: max |diff| %.4f and %.4f" % (dev.loc[3].max(), dev.loc[4].max()))
print("NaN cells (N, rate):", [(int(bra.N_S2[a]), float(R[b])) for a, b in np.argwhere(np.isnan(s2))])
# Monte Carlo error: each N is its own scan (performance_scan_figureS2.py, 10,000 runs per rate),
# independent of the reference runs (figure2i scan), so a difference has SE about sqrt(2) SE.
se2 = bra.resultant_se(ref2)
zs2 = (s2 - ref2) / (np.sqrt(2) * se2[0])
print("SE of the reference at 0.1/s: %.4f (bound %.4f); a difference of two independent values: %.4f"
% (se2[0][2], se2[1][2], np.sqrt(2) * se2[0][2]))
print("N = 16 at 0.1/s: %.4f, %.1f SEs of a difference; |z| above 3 elsewhere: %d of %d"
% ((s2 - ref2)[3, 2], zs2[3, 2], int(np.sum(np.abs(zs2) > 3)) - 1, int(np.isfinite(zs2).sum()) - 1))
print("at 0.1/s, every N minus the reference:", dict(zip(bra.N_S2.tolist(), (s2[:, 2] - ref2[2]).round(4))))
print("N = 16 at 0.1/s minus the mean of the other six N: %.4f" % (s2[3, 2] - np.delete(s2[:, 2], 3).mean()))
expected = np.nansum(2 * stats.norm.sf(0.02 / (np.sqrt(2) * se2[0]))) * len(bra.N_S2[bra.N_S2 >= 8])
print("cells with N >= 8 expected beyond 0.02 by sampling alone, if every network equalled the reference: %.2f" % expected)
vonMises_q_stats.npz equals the attractor's accuracy in figure2i.npz: True
missing values (NaN) per N: {3: 0, 4: 0, 8: 0, 16: 0, 32: 0, 64: 0, 128: 1}
max |diff| per N: {3: 0.0105, 4: 0.0186, 8: 0.0183, 16: 0.0414, 32: 0.009, 64: 0.0142, 128: 0.0129}
cells above 0.02: [(16, 0.1)]
N >= 8: max |diff| 0.0414, within 0.02 at every rate: False
N >= 16: max |diff| 0.0414, within 0.02 at every rate: False
N = 3 and 4, reported: max |diff| 0.0105 and 0.0186
NaN cells (N, rate): [(128, 100.0)]
SE of the reference at 0.1/s: 0.0068 (bound 0.0098); a difference of two independent values: 0.0096
N = 16 at 0.1/s: -0.0414, -4.3 SEs of a difference; |z| above 3 elsewhere: 0 of 159
at 0.1/s, every N minus the reference: {3: np.float64(-0.0105), 4: np.float64(-0.0186), 8: np.float64(-0.0183), 16: np.float64(-0.0414), 32: np.float64(-0.0017), 64: np.float64(-0.0142), 128: np.float64(-0.0129)}
N = 16 at 0.1/s minus the mean of the other six N: -0.0287
cells with N >= 8 expected beyond 0.02 by sampling alone, if every network equalled the reference: 1.19
Result: not reproduced. Of the 7 x 23 values, one is more than 0.02 from the reference: N = 16 at 0.1/s, 0.041 below it. All others, N = 3 and 4 included, are within 0.019. Both readings (N of 8 and more, N of 16 and more) include N = 16, so neither holds, and the criterion requires both. The one deviation is not sampling error: it is 4.3 standard errors of a difference of two independent scans (0.0096), no other value is beyond 3, and it is 0.029 below the mean of the other six sizes at the same rate. Sampling alone would put about 1.2 of the values with N of 8 and more beyond 0.02 if every network equalled the reference, so the bound is near the scans' sampling error, but this value is far beyond it. At every other network size and rate the accuracy is close to the reference. One value (N = 128 at 100/s) is missing in the stored array.
prediction-bump-amplitude-cue-reliability¶
The paper predicts that the bump amplitude in the HD system is higher in the presence of reliable external HD cues, such as salient visual landmarks. The paper does not test this.
Bump amplitude is the length of the population vector at the bump's spatial frequency across the PB midline ROIs (two bump copies across the bridge, frequency 2, the peak of every fly's spatial spectrum), scaled by 2/n as in the paper's code. More reliable cue = high contrast, less reliable = low contrast, as the trials table labels them. Per fly, the mean amplitude over all frames of each contrast. Signal (i) is the stored dF/F. Signal (ii) is raw fluorescence minus baseline: the TwoPhotonSeries averaged over planes, smoothed in the plane (Gaussian SD 4 pixels) and averaged in each ROI mask, minus each ROI's 10th percentile over the session, read in blocks of 25 volumes every 250 (10% of the frames, about 0.9 GB per fly). The projection and smoothing were chosen by their agreement with the stored dF/F, shown first.
print("spatial frequency of the peak of the dF/F spectrum, per fly:",
[int(bra.spectrum(s["dff"]).argmax()) for s in sessions])
# which projection reproduces the stored dF/F? sub-1, 10 blocks of 25 volumes
s1 = sessions[0]
f1 = bra.open_nwb(contrast_assets[0])
blocks = [(a, a + 25) for a in range(0, 2500, 250)]
idx = np.concatenate([np.arange(a, b) for a, b in blocks])
for proj, sm in [("mean", 0), ("max", 0), ("mean", 2), ("mean", 4), ("max", 4)]:
F = bra.raw_roi_fluorescence(f1, blocks, s1["masks"], proj, sm)
r_ = [np.corrcoef(F[:, j], s1["dff"][idx, j])[0, 1] for j in range(F.shape[1])]
print("projection %-4s smoothing %d px: median correlation with stored dF/F over ROIs %.3f (min %.3f)"
% (proj, sm, np.median(r_), np.min(r_)))
spatial frequency of the peak of the dF/F spectrum, per fly: [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
projection mean smoothing 0 px: median correlation with stored dF/F over ROIs 0.659 (min 0.135) projection max smoothing 0 px: median correlation with stored dF/F over ROIs 0.595 (min 0.076)
projection mean smoothing 2 px: median correlation with stored dF/F over ROIs 0.894 (min 0.233)
projection mean smoothing 4 px: median correlation with stored dF/F over ROIs 0.925 (min 0.248)
projection max smoothing 4 px: median correlation with stored dF/F over ROIs 0.910 (min 0.216)
def amplitude_by_contrast(s, act, t):
_, amp = bra.fourier(act)
out = {}
for c in (1, 2, 3):
m = np.zeros(len(t), bool)
for a, b, l in zip(s["trials"]["start_time"], s["trials"]["stop_time"], s["trials"]["cue_contrast"]):
if l == c:
m |= (t >= a) & (t < b)
out[c] = amp[m].mean()
return out
def raw_signal(a, s):
f = bra.open_nwb(a)
n = s["tps_shape"][0]
blocks = [(i, min(i + 25, n)) for i in range(0, n, 250)]
F = bra.raw_roi_fluorescence(f, blocks, s["masks"], "mean", (0, 2, 4))
t = s["t"][np.concatenate([np.arange(i, j) for i, j in blocks])]
return F, t
t1 = time.time()
with ThreadPoolExecutor(4) as ex:
raws = list(ex.map(raw_signal, contrast_assets, sessions))
print("raw fluorescence read in %.0f s" % (time.time() - t1))
rows = []
for s, (Fs, tF) in zip(sessions, raws):
F = Fs[4] - np.percentile(Fs[4], 10, axis=0) # the main raw signal: smoothing 4 px, 10th percentile
a_d = amplitude_by_contrast(s, s["dff"], s["t"])
a_r = amplitude_by_contrast(s, F, tF)
rows.append(dict(fly=s["subject"], dff_zero=a_d[1], dff_low=a_d[2], dff_high=a_d[3],
raw_zero=a_r[1], raw_low=a_r[2], raw_high=a_r[3]))
amp = pd.DataFrame(rows).set_index("fly")
print(amp.round(3).to_string())
pred = {}
for sig in ("dff", "raw"):
hi, lo, ze = amp[sig + "_high"], amp[sig + "_low"], amp[sig + "_zero"]
n_hi = int((hi > lo).sum())
w = stats.wilcoxon(hi, lo, alternative="greater")
pred[(sig, "a")] = n_hi > len(amp) / 2
pred[(sig, "b")] = w.pvalue < 0.05
print("\n%s: high > low in %d of %d flies; median high/low ratio %.3f; median difference %.3f; one-sided Wilcoxon W = %g, p = %.4f"
% (sig, n_hi, len(amp), np.median(hi / lo), np.median(hi - lo), w.statistic, w.pvalue))
w0 = stats.wilcoxon(hi, ze, alternative="greater")
print(" variant, high against zero contrast: %d of %d flies higher, median ratio %.3f, p = %.4f"
% (int((hi > ze).sum()), len(amp), np.median(hi / ze), w0.pvalue))
w1 = stats.wilcoxon(lo, ze, alternative="greater")
print(" variant, low against zero contrast: %d of %d flies higher, p = %.4f" % (int((lo > ze).sum()), len(amp), w1.pvalue))
print(" without sub-2 (one high block): high > low in %d of 14, p = %.4f"
% (int((hi.drop("2") > lo.drop("2")).sum()), stats.wilcoxon(hi.drop("2"), lo.drop("2"), alternative="greater").pvalue))
print("\nreadings:", pred)
# raw-signal variants (the verifier's request): in-plane smoothing 0, 2, 4 px and baseline at the
# 5th, 10th or 20th percentile of each ROI, plus the baseline-normalized signal (F - F0) / F0
rows = []
for sm in (0, 2, 4):
for pct in (5, 10, 20):
for norm_ in ("F - F0", "(F - F0) / F0"):
hl = []
for s, (Fs, tF) in zip(sessions, raws):
F0 = np.percentile(Fs[sm], pct, axis=0)
X = Fs[sm] - F0 if norm_ == "F - F0" else (Fs[sm] - F0) / F0
a_ = amplitude_by_contrast(s, X, tF)
hl.append((a_[3], a_[2]))
hl = np.array(hl)
rows.append(dict(smooth_px=sm, baseline_pct=pct, signal=norm_, high_gt_low=int((hl[:, 0] > hl[:, 1]).sum()),
median_ratio=np.median(hl[:, 0] / hl[:, 1]),
p_wilcoxon=stats.wilcoxon(hl[:, 0], hl[:, 1], alternative="greater").pvalue))
rv = pd.DataFrame(rows)
print("\nraw-signal variants, high against low contrast over 15 flies:")
print(rv.round(4).to_string(index=False))
print("variants with p < 0.05: %d of %d; with high > low in more than half the flies: %d of %d"
% ((rv.p_wilcoxon < 0.05).sum(), len(rv), (rv.high_gt_low > 7).sum(), len(rv)))
raw fluorescence read in 402 s
dff_zero dff_low dff_high raw_zero raw_low raw_high
fly
1 0.471 0.461 0.538 45.047 40.951 45.534
2 0.572 0.567 0.722 54.612 55.473 66.323
3 0.631 0.667 0.680 40.693 48.959 48.962
4 0.802 0.823 0.812 143.750 135.246 149.707
5 0.301 0.292 0.327 25.332 17.785 26.943
6 0.838 0.748 0.678 60.155 56.581 51.814
7 0.393 0.456 0.552 27.050 34.946 43.872
8 0.536 0.536 0.600 36.313 47.755 51.333
9 0.486 0.433 0.553 46.099 44.369 59.152
10 0.646 0.674 0.473 63.179 74.403 50.905
11 0.516 0.503 0.570 70.101 65.238 86.060
12 0.539 0.616 0.703 49.545 60.033 69.591
13 0.731 0.890 0.779 74.765 74.568 84.380
14 0.409 0.597 0.545 35.524 48.836 41.587
15 0.721 0.655 0.710 88.856 82.788 90.303
dff: high > low in 10 of 15 flies; median high/low ratio 1.120; median difference 0.055; one-sided Wilcoxon W = 80, p = 0.1384
variant, high against zero contrast: 12 of 15 flies higher, median ratio 1.105, p = 0.0473
variant, low against zero contrast: 7 of 15 flies higher, p = 0.2622
without sub-2 (one high block): high > low in 9 of 14, p = 0.2131
raw: high > low in 12 of 15 flies; median high/low ratio 1.112; median difference 8.926; one-sided Wilcoxon W = 96, p = 0.0206
variant, high against zero contrast: 13 of 15 flies higher, median ratio 1.171, p = 0.0062
variant, low against zero contrast: 7 of 15 flies higher, p = 0.2106
without sub-2 (one high block): high > low in 11 of 14, p = 0.0338
readings: {('dff', 'a'): True, ('dff', 'b'): np.False_, ('raw', 'a'): True, ('raw', 'b'): np.True_}
raw-signal variants, high against low contrast over 15 flies:
smooth_px baseline_pct signal high_gt_low median_ratio p_wilcoxon
0 5 F - F0 11 1.0894 0.0416
0 5 (F - F0) / F0 7 0.9976 0.6807
0 10 F - F0 11 1.0991 0.0319
0 10 (F - F0) / F0 7 0.9551 0.8853
0 20 F - F0 11 1.1029 0.0365
0 20 (F - F0) / F0 9 1.1151 0.1147
2 5 F - F0 11 1.1248 0.0277
2 5 (F - F0) / F0 9 1.0736 0.2444
2 10 F - F0 11 1.1175 0.0277
2 10 (F - F0) / F0 9 1.0780 0.2271
2 20 F - F0 12 1.1133 0.0240
2 20 (F - F0) / F0 10 1.1017 0.1262
4 5 F - F0 11 1.1117 0.0240
4 5 (F - F0) / F0 10 1.0638 0.1651
4 10 F - F0 12 1.1119 0.0206
4 10 (F - F0) / F0 10 1.0596 0.1514
4 20 F - F0 12 1.1264 0.0206
4 20 (F - F0) / F0 10 1.0932 0.1147
variants with p < 0.05: 9 of 18; with high > low in more than half the flies: 16 of 18
Result: blocked (methods_underspecified). The readings disagree. Bump amplitude is higher at high than at low contrast in 10 of 15 flies with dF/F and in 12 of 15 with raw fluorescence minus baseline, so reading (a) holds for both signals. The one-sided Wilcoxon test over flies (reading (b)) gives p = 0.021 with raw fluorescence but p = 0.138 with dF/F (median high/low ratio 1.11 and 1.12). The criterion blocks the claim when its readings differ. Variants: against zero contrast, the high-contrast amplitude is higher in 12 (dF/F, p = 0.047) and 13 of 15 flies (raw, p = 0.006); low contrast is not above zero contrast (7 of 15 flies with either signal). Across 18 variants of the raw signal (smoothing 0, 2 or 4 pixels; baseline at the 5th, 10th or 20th percentile; F - F0 or (F - F0) / F0), high exceeds low contrast in more than half the flies in 16, and p < 0.05 in 9: in all nine F - F0 variants (p 0.021 to 0.042) and in none of the nine (F - F0) / F0 variants (p 0.11 to 0.89). The split between the readings follows the division by the baseline, not the smoothing or the percentile.
Summary¶
| Claim | Outcome |
|---|---|
| fig1-conflict-decreases-certainty | reproduced |
| fig2h-bump-tracks-circkf | reproduced |
| fig2i-accuracy-and-certainty-vs-reliability | reproduced |
| fig3b-bayesian-beats-conventional | reproduced |
| fig3c-update-weights | reproduced |
| fig3d-best-network-slow-decay | partially reproduced |
| fig3f-noise-drop-below-70 | reproduced |
| fig3f-retuning-reduces-drop | reproduced (34.6%) |
| noise-effects-on-mu-and-kappa | partially reproduced |
| fig4b-connectome-motifs | partially reproduced, assumed (one sort of two, any synapse) |
| fig4gh-drosophila-network-matches | partially reproduced |
| figs2-performance-independent-of-n | not reproduced |
| prediction-bump-amplitude-cue-reliability | blocked, methods underspecified |
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.
How the analysis differs from the paper's. Twelve claims are tested on the paper's own released
outputs and code, not on the dandiset, as the claims file sets out. The scan scripts were not
re-run (10,000 trajectories per setting is beyond the time limit); their Monte Carlo error is
estimated from the stored accuracies. The noise test uses a batched
copy of RNN_filter_Run, checked against it, with a variant that holds the global inhibition at its noise-free level. The connectome is the public v1.2 export, not
v1.2.1 in neuPrint, and only the PB-glomerulus sort can be formed from its names. The prediction is
tested on the dandiset's cue-contrast experiment, from a different study; the raw-fluorescence signal
is read on 10% of the frames.
What was found about the dataset. The dandiset holds no data of this paper, which is a theory
paper; its related resources point to a bioRxiv preprint of this paper (version 4) for a study whose
data it holds (Basnak et al., in preparation). In the cue-contrast experiment, imaging_rate says
9.18 Hz for every file, while sub-4 and sub-10 are imaged at 7.87 Hz (0.127 s frames, 14 planes);
sub-2 has one high-contrast block, not two, and a last block of 405 s labelled low contrast. The
ROI masks are 3 to 11 pixels, and the stored dF/F agrees with the raw imaging
only after spatial smoothing, which the files do not describe.
print("notebook ran in %.1f min" % ((time.time() - started) / 60))
notebook ran in 8.4 min