Reproducing Fan et al. (2024) from dandiset 000959¶
Paper. Fan, Dal Monte, Nair, Fagan & Chang (2024), "Closed-loop microstimulations of the orbitofrontal cortex during real-life gaze interaction enhance dynamic social attention", Neuron 112, 2631-2644, doi:10.1016/j.neuron.2024.05.004.
Claims. claims/000959/10.1016-j.neuron.2024.05.004.yaml, written from the accepted author
manuscript (PMC11309918) with its supplement, and the paper's code
(github.com/changlabneuro/microstim-ofc-social-attention at cc0bf994).
Dataset. DANDI:000959, version 0.261020.0242. The paper cites the dandiset without a version. Two versions are published, 0.260531.2104 and 0.261020.0242; the most recent is pinned, as the protocol asks when the paper cites none. The data check below compares the two asset lists: both hold the same 1705 assets with the same identifiers, paths and sizes, so the choice does not change the data.
Each file is one 300-s run (a "session" in the paper) and holds M1's gaze and pupil at 1 kHz,
M2's gaze and pupil in social runs, the stimulation and sham times, the day's ROIs, and
fixation events as interval tables (M1-only, M2-only and mutual events by ROI). Files are read
over HTTP into memory at the pinned version (nothing is written to disk) and reduced to small
tables by microstim.py.
| Claim | Tested |
|---|---|
| design-days-per-area, design-sessions-per-day | tested |
| fig1e-stim-counts-comparable, figs1d-stim-counts-nonsocial | tested |
| figs1e-eyes-more-than-face, figs1e-eyes-vs-rdm-comparable | tested |
| fig2a-density-map-ofc | tested |
| fig2c-ofc-social-gaze-distance, fig2c-dmpfc-accg-no-distance-effect | tested |
| fig2e-ofc-contralateral, fig2e-dmpfc-accg-no-hemifield-effect | tested |
| fig2f-nonsocial-no-distance-effect, fig2gh-nonsocial-no-hemifield-effect | tested |
| methods-ofc-distance-longer-windows, methods-ofc-no-change-in-fixation-number, methods-ofc-longer-face-looking | tested |
| methods-mouth-stimulation-no-effect, methods-eyes-vs-mouth-bootstrap | tested |
| fig3b-ofc-inter-looking-interval, ili-nonsocial-no-effect, ili-missing-trials-41, figs3b-ili-trial-level | tested |
| figs3c-reciprocation-day-level-none, fig3d-ofc-reciprocation-latency | tested |
| partner-monkey-unchanged | tested |
| epoch-days-76-of-81 | tested |
| figs4a-causal-strength-no-stim-effect, figs4b-causal-strength-no-epoch-effect, ofc-slope-unchanged-by-epoch, fig4c-dmpfc-slope-late-greater, fig4d-dmpfc-slope-permutation, dmpfc-slope-positive | blocked, too_large (the causal decomposition; timed below) |
| fig5a-trigger-duration-unchanged, fig5b-pupil-unchanged | tested |
| fig5cd-saccade-counts-unchanged, fig5e-h-saccade-kinematics-unchanged | tested, with outside software (a Python port of the Otero-Millan et al. 2014 saccade detector, pinned; the dataset holds no saccades) |
import time
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
from scipy.stats import mannwhitneyu, wilcoxon
from statsmodels.stats.multitest import multipletests
import microstim as ms
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
AREAS = ["OFC", "dmPFC", "ACCg"]
def signrank(a, b=None):
"""Two-sided Wilcoxon signed-rank test (scipy's default method: exact for small n
without ties). Returns (n, p)."""
a = np.asarray(a, float)
d = a if b is None else a - np.asarray(b, float)
d = d[~np.isnan(d)]
if len(d) == 0 or np.all(d == 0):
return len(d), np.nan
return len(d), wilcoxon(d).pvalue
def ranksum(a, b):
a = np.asarray(a, float)
b = np.asarray(b, float)
a, b = a[~np.isnan(a)], b[~np.isnan(b)]
return mannwhitneyu(a, b, alternative="two-sided").pvalue
def day_test(dm, region, subject=None):
"""Signed-rank of stim against sham day means; returns n, median stim - sham, p."""
x = dm[(dm.region == region) & ((dm.subject == subject) if subject else True)]
x = x.dropna(subset=["stim", "sham"])
n, p = signrank(x.stim, x.sham)
return dict(region=region, subject=subject or "both", n=n,
median_diff=np.median(x.stim - x.sham) if n else np.nan,
mean_diff=np.mean(x.stim - x.sham) if n else np.nan, p=p)
def day_table(dm, subjects=(None, "L", "T"), areas=AREAS):
return pd.DataFrame([day_test(dm, a, s) for a in areas for s in subjects])
Data check¶
t_start = time.time()
assets = ms.list_assets()
other = ms.list_assets("0.260531.2104")
same = assets.merge(other, on=["path", "asset_id", "size"], how="outer", indicator=True)
print(f"version {ms.DANDISET_VERSION}: {len(assets)} assets, {assets['size'].sum() / 1e9:.1f} GB; "
f"version 0.260531.2104: {len(other)} assets; rows not shared: {(same._merge != 'both').sum()}")
runs = ms.load_all(assets)
meta, trials, events = ms.prepare(runs)
print(f"loaded {len(runs)} runs in {(time.time() - t_start) / 60:.1f} min")
lay = ms.file_layout(assets.asset_id.iloc[0])
print("groups of one file,", assets.path.iloc[0])
for g, names in lay.items():
print(f" {g}: {names}")
print("tables named like saccade in any group of that file:",
[n for names in lay.values() if names for n in names if "sacc" in n.lower()])
print("\nThe description, the scratch fields and the file name agree:",
bool((meta.region == meta.region_scratch).all() and (meta.roi == meta.roi_scratch).all()
and (meta.cond == meta.cond_scratch).all()
and (meta.path.str.slice(4, 5) == meta.subject).all()))
print("runs by subject, triggering ROI, region and condition:")
print(meta.groupby(["subject", "roi", "region", "cond"]).size().unstack("cond"))
print("\ndays by subject, triggering ROI and region:")
print(meta.groupby(["subject", "roi", "region"]).date.nunique().unstack("region"))
print("days with more than one region:", (meta.groupby("date").region.nunique() > 1).sum())
dup_key = meta.apply(lambda r: (r.n_samples, round(r.t_first, 3), round(r.t_last, 3)), axis=1)
dup_stim = trials.groupby("run_id").time.apply(lambda s: tuple(np.round(s, 3)))
print("runs sharing sample count and time span:", dup_key.duplicated().sum(),
"; runs sharing all trial times:", dup_stim.duplicated().sum(),
"; identifiers shared by a Social and a Dot run (run numbers restart per condition):",
meta.identifier.duplicated().sum(), "; (date, condition, run) repeated:",
meta.duplicated(["date", "cond", "run"]).sum())
print("\nunits:", meta.gaze_unit.unique(), "; sampling rate (Hz):", meta.rate.unique(),
"; median and largest step of t (s):", meta.dt_median.median(), meta.dt_max.max())
print("valid run length (s) by condition:")
print((meta.t_last - meta.t_first).groupby(meta.cond).describe()[["count", "min", "50%", "max"]])
print("M1 gaze x range (pixels):", np.min([r[0] for r in meta.gaze_x_range]), np.max([r[1] for r in meta.gaze_x_range]),
"; y range:", np.min([r[0] for r in meta.gaze_y_range]), np.max([r[1] for r in meta.gaze_y_range]),
"; median fraction of NaN gaze in the valid run:", round(meta.gaze_nan.median(), 3))
print("M1 pupil range (a.u.):", np.min([r[0] for r in meta.pupil_range]), np.max([r[1] for r in meta.pupil_range]))
print("times (s): first valid sample", round(meta.t_first.min(), 3), "to", round(meta.t_first.max(), 3),
"; last valid sample", round(meta.t_last.min(), 3), "to", round(meta.t_last.max(), 3),
"; event starts", round(events.start.min(), 3), "to", round(events.start.max(), 3),
"; trial times", round(trials.time.min(), 3), "to", round(trials.time.max(), 3),
"(all below the 1e6 s run offset used by window_pairs)")
print("event start indices are 1-based: largest |t[start_index - 1] - start time| =",
meta.index_time_max_err.max(), "s")
print("\ntrials by condition and label (all events in the files):")
print(trials.groupby(["cond", "roi", "type"]).size().unstack("type"))
gap = trials.sort_values(["run_id", "time"]).groupby("run_id").time.diff()
print("smallest interval between consecutive trials (s):", round(gap.min(), 3),
"; intervals under 4.99 s:", int((gap < 4.99).sum()))
print("kept for analysis:", int(trials.kept.sum()), "of", len(trials), "; dropped as first of run:",
int(trials["first"].sum()), "; outside the run or in its last 5 s:", int((~trials.in_run & ~trials["first"]).sum()))
blk = trials.sort_values(["run_id", "time"]).copy()
blk["k"] = blk.groupby("run_id").cumcount() // 4
full = blk.groupby(["run_id", "k"]).filter(lambda g: len(g) == 4)
print("blocks of four consecutive trials with two stimulations:",
round((full.groupby(["run_id", "k"]).type.apply(lambda s: (s == "stim").sum()) == 2).mean(), 3))
m1 = events[events.who == "m1"]
print("\nevent tables: m1_all = m1_everywhere + m1_wholeface and m1_wholeface = m1_eye + m1_neface in every run:",
bool((m1.groupby(["run_id", "table"]).size().unstack(fill_value=0).eval(
"m1_all - m1_everywhere - m1_wholeface").abs().max() == 0)))
mw = set(zip(m1[m1.table == "mutual_wholeface"].run_id, m1[m1.table == "mutual_wholeface"].start.round(4)))
for t in ["mutual_eye", "mutual_face", "mutual_mouth"]:
s = set(zip(m1[m1.table == t].run_id, m1[m1.table == t].start.round(4)))
print(f" {t}: {len(s)} events, in mutual_wholeface: {len(s & mw)}")
print(" mutual_wholeface:", int((m1.table == "mutual_wholeface").sum()), "events; mutual_eye + mutual_face + mutual_mouth:",
sum(int((m1.table == t).sum()) for t in ["mutual_eye", "mutual_face", "mutual_mouth"]),
"(mutual_face is the face outside Eyes and Mouth)")
mo = set(zip(m1[m1.table == "m1_mouth"].run_id, m1[m1.table == "m1_mouth"].start.round(4)))
ne = set(zip(m1[m1.table == "m1_neface"].run_id, m1[m1.table == "m1_neface"].start.round(4)))
print(" m1_mouth events:", len(mo), "of which in m1_neface:", len(mo & ne), "(m1_neface is the face outside Eyes, Mouth included)")
alls = set(zip(m1[m1.table == "m1_all"].run_id, m1[m1.table == "m1_all"].start.round(4)))
print(" M1 sides of mutual events also in m1_all:", len(mw & alls))
print("M1 events with no gaze position (all samples NaN; left out of distances):",
int(m1.x.isna().sum()), "of", len(m1), "; M2 events:", int(events[events.who == "m2"].x.isna().sum()))
lab = m1.dropna(subset=["hemi"])
print("hemifield from the centre of roi_m2_face agrees with the files' laterality_label for",
f"{(lab.hemi == lab.laterality_label).mean():.4f} of M1 events "
f"(social {(lab[lab.cond == 'Social'].hemi == lab[lab.cond == 'Social'].laterality_label).mean():.4f}, "
f"Dot {(lab[lab.cond == 'Dot'].hemi == lab[lab.cond == 'Dot'].laterality_label).mean():.4f})")
# fixation positions against their ROI labels
rows = []
for (cond, table), g in m1[m1.table.isin(["m1_eye", "m1_wholeface", "m1_everywhere"])].groupby(["cond", "table"]):
r = []
for rid, gg in g.groupby("run_id"):
ro = runs[rid]["rois"]
rect = ro["roi_m1_eyes"] if table == "m1_eye" else ro["roi_m1_face"]
r.append(ms.in_rect(gg[["x", "y"]].to_numpy(), rect))
rows.append(dict(cond=cond, table=table, n=len(g), mean_position_in_rectangle=np.concatenate(r).mean()))
print("\nfraction of events whose mean position lies in the ROI rectangle of the file "
"(Eyes for m1_eye, whole Face for the others):")
print(pd.DataFrame(rows))
nvar = [d for d, g in meta.groupby("date") if len({tuple(np.round(runs[i]["rois"]["roi_m1_eyes"], 2)) for i in g.run_id}) > 1]
print("days whose Eyes ROI differs between runs (recalibrated; each run's own ROIs are used):", nvar)
version 0.261020.0242: 1705 assets, 26.0 GB; version 0.260531.2104: 1705 assets; rows not shared: 0
loaded 1705 runs in 18.0 min
groups of one file, sub-L/sub-L_ses-20190718_obj-1wqjp5d_behavior.nwb
/: ['acquisition', 'analysis', 'file_create_date', 'general', 'identifier', 'intervals', 'processing', 'scratch', 'session_description', 'session_start_time', 'specifications', 'stimulus', 'timestamps_reference_time']
intervals: ['m1_all', 'm1_everywhere', 'm1_eye', 'm1_mouth', 'm1_neface', 'm1_nenmface', 'm1_wholeface']
processing: []
acquisition: ['EyeTracking_m1', 'pupil_m1', 'roi_definitions', 'stimulation_events', 't']
scratch: ['stimulated_region', 'stimulation_condition', 'stimulation_roi']
analysis: []
tables named like saccade in any group of that file: []
The description, the scratch fields and the file name agree: True
runs by subject, triggering ROI, region and condition:
cond Dot Social
subject roi region
L Eyes ACCg 91 150
OFC 92 150
dmPFC 88 150
T Eyes ACCg 60 168
OFC 60 168
dmPFC 60 168
Mouth ACCg 25 75
OFC 25 75
dmPFC 25 75
days by subject, triggering ROI and region:
region ACCg OFC dmPFC
subject roi
L Eyes 15 15 15
T Eyes 12 12 12
Mouth 5 5 5
days with more than one region: 0
runs sharing sample count and time span: 0 ; runs sharing all trial times: 0 ; identifiers shared by a Social and a Dot run (run numbers restart per condition): 526 ; (date, condition, run) repeated: 0
units: <StringArray>
['pixels']
Length: 1, dtype: str ; sampling rate (Hz): [1000.] ; median and largest step of t (s): 0.0010000000002037268 0.0010049950050188272
valid run length (s) by condition:
count min 50% max
cond
Dot 526.0 295.839 295.841 297.290
Social 1179.0 295.872 296.497 296.701
M1 gaze x range (pixels): 0.20001220703125 3072.10009765625 ; y range: -768.0999755859375 1536.0999755859375 ; median fraction of NaN gaze in the valid run: 0.062
M1 pupil range (a.u.): 0.0 6258.0
times (s): first valid sample 13.757 to 9659.214 ; last valid sample 310.239 to 9955.055 ; event starts 13.757 to 9953.872 ; trial times 15.554 to 9945.824 (all below the 1e6 s run offset used by window_pairs)
event start indices are 1-based: largest |t[start_index - 1] - start time| = 0.0 s
trials by condition and label (all events in the files):
type sham stim
cond roi
Dot Eyes 5408 5403
Mouth 1033 1034
Social Eyes 10443 10395
Mouth 1355 1352
smallest interval between consecutive trials (s): 4.999 ; intervals under 4.99 s: 0
kept for analysis: 34019 of 36423 ; dropped as first of run: 1705 ; outside the run or in its last 5 s: 699
blocks of four consecutive trials with two stimulations: 0.995
event tables: m1_all = m1_everywhere + m1_wholeface and m1_wholeface = m1_eye + m1_neface in every run: True
mutual_eye: 2912 events, in mutual_wholeface: 2912 mutual_face: 2384 events, in mutual_wholeface: 2384
mutual_mouth: 263 events, in mutual_wholeface: 263
mutual_wholeface: 5896 events; mutual_eye + mutual_face + mutual_mouth: 5896 (mutual_face is the face outside Eyes and Mouth)
m1_mouth events: 23641 of which in m1_neface: 23641 (m1_neface is the face outside Eyes, Mouth included)
M1 sides of mutual events also in m1_all: 0
M1 events with no gaze position (all samples NaN; left out of distances): 0 of 1929396 ; M2 events: 188
hemifield from the centre of roi_m2_face agrees with the files' laterality_label for 0.9998 of M1 events (social 0.9997, Dot 1.0000)
fraction of events whose mean position lies in the ROI rectangle of the file (Eyes for m1_eye, whole Face for the others):
cond table n mean_position_in_rectangle
0 Dot m1_everywhere 247506 0.000000
1 Dot m1_eye 24703 0.893009
2 Dot m1_wholeface 46318 0.900751
3 Social m1_everywhere 546535 0.000000
4 Social m1_eye 38087 0.937065
5 Social m1_wholeface 51760 0.969494
days whose Eyes ROI differs between runs (recalibrated; each run's own ROIs are used): ['08142019', '12102019', '12132019']
Clock check. Each trial was triggered online by an M1 fixation in the targeted ROI, so the trial times should fall inside the M1 events labelled with that ROI far more often than the same times shifted by a few hundred milliseconds or seconds. The contrast between the true and shifted times tests that the trial times and the event tables share a clock.
def inside_rate(tr, ev, shift):
p = ms.window_pairs(tr.assign(time=tr.time + shift), ev, -10.0, 0.0005)
p = p[p.stop >= p.start - p.lag]
return tr.index.isin(p.trial).mean()
clock = []
for (cond, roi), tr in trials[trials.kept].groupby(["cond", "roi"]):
tab = "eye" if roi == "Eyes" else "mouth"
ev = events[(events.who == "m1") & (events.target == tab)]
for shift in [-2.0, -0.5, 0.0, 0.5, 2.0]:
clock.append(dict(cond=cond, roi=roi, shift_s=shift, fraction_inside=inside_rate(tr, ev, shift)))
clock = pd.DataFrame(clock).pivot(index=["cond", "roi"], columns="shift_s", values="fraction_inside")
print("fraction of trials whose time lies inside an M1 event in the triggering ROI, by time shift:")
print(clock.round(3))
by_run = trials[trials.kept & (trials.cond == "Social") & (trials.roi == "Eyes")]
ev_e = events[(events.who == "m1") & (events.target == "eye")]
pr = ms.window_pairs(by_run, ev_e, -10.0, 0.0005)
pr = pr[pr.stop >= pr.start - pr.lag]
r_in = by_run.assign(inside=by_run.index.isin(pr.trial)).groupby("run_id").inside.mean()
print("social Eyes runs: fraction inside per run, quantiles:", r_in.quantile([0, 0.05, 0.5]).round(3).to_dict(),
"; runs below 0.5:", int((r_in < 0.5).sum()), "of", len(r_in))
fraction of trials whose time lies inside an M1 event in the triggering ROI, by time shift:
shift_s -2.0 -0.5 0.0 0.5 2.0
cond roi
Dot Eyes 0.147 0.152 0.619 0.264 0.192
Mouth 0.018 0.031 0.333 0.125 0.058
Social Eyes 0.054 0.055 0.708 0.200 0.100
Mouth 0.011 0.023 0.280 0.130 0.048
social Eyes runs: fraction inside per run, quantiles: {0.0: 0.261, 0.05: 0.467, 0.5: 0.722} ; runs below 0.5: 65 of 954
design-days-per-area¶
Eyes stimulation was applied on 27 days per stimulated area (OFC, dmPFC and ACCg), 15 per area from monkey L and 12 per area from monkey T, and Mouth stimulation on 5 days per area from monkey T.
kept = trials[trials.kept]
dd = kept.groupby(["subject", "roi", "region", "date"]).type.agg(lambda s: set(s) >= {"stim", "sham"})
days = dd[dd].reset_index().groupby(["subject", "roi", "region"]).date.nunique().unstack("region")
print("days holding kept stimulation and sham trials:")
print(days)
days holding kept stimulation and sham trials:
region ACCg OFC dmPFC
subject roi
L Eyes 15 15 15
T Eyes 12 12 12
Mouth 5 5 5
Result: reproduced. Every day holds kept stimulation and sham trials: 15 Eyes days per area for monkey L, 12 per area for monkey T, and 5 Mouth days per area for monkey T.
design-sessions-per-day¶
Each day held on average 10 live social gaze sessions (9 to 11) and 5 non-social control sessions (4 to 8) for monkey L, and 15 social (14 to 15) and 5 control sessions for monkey T, each 300 s long.
eyes_days = meta[meta.roi == "Eyes"]
spd = eyes_days.groupby(["subject", "date", "cond"]).size().unstack("cond")
print(spd.groupby("subject").agg(["mean", "min", "max"]).round(2))
print("\nnumber of control runs per day, monkey L:", spd.loc["L", "Dot"].value_counts().sort_index().to_dict())
d8 = spd.loc["L"][spd.loc["L", "Dot"] == 8].index
l8 = eyes_days[(eyes_days.subject == "L") & eyes_days.date.isin(d8) & (eyes_days.cond == "Dot")]
print("monkey L days with 8 control runs:", len(d8), "; their control run numbers:",
sorted({tuple(sorted(g.run)) for _, g in l8.groupby("date")}),
"; run length (s)", round((l8.t_last - l8.t_first).min(), 1), "to", round((l8.t_last - l8.t_first).max(), 1),
"; kept trials per run", trials[trials.kept & trials.run_id.isin(l8.run_id)].groupby("run_id").size().agg(["min", "median", "max"]).to_dict())
print("number of social runs per day, monkey T:", spd.loc["T", "Social"].value_counts().sort_index().to_dict())
dur = eyes_days.t_last - eyes_days.t_first
print(f"\nrun length (s), valid samples: mean {dur.mean():.1f}, min {dur.min():.1f}, max {dur.max():.1f}; "
f"whole file: mean {(eyes_days.n_samples / 1000).mean():.1f} s")
spd_m = meta[meta.roi == "Mouth"].groupby(["subject", "date", "cond"]).size().unstack("cond")
print("Mouth days (not part of the claim): runs per day", spd_m.agg(["mean", "min", "max"]).to_dict())
cond Dot Social
mean min max mean min max
subject
L 6.02 4 8 10.0 9 11
T 5.00 5 5 14.0 14 14
number of control runs per day, monkey L: {4: 2, 5: 24, 6: 4, 7: 1, 8: 14}
monkey L days with 8 control runs: 14 ; their control run numbers: [(1, 2, 3, 4, 5, 6, 7, 8)] ; run length (s) 295.8 to 296.9 ; kept trials per run {'min': 9.0, 'median': 22.0, 'max': 32.0}
number of social runs per day, monkey T: {14: 36}
run length (s), valid samples: mean 296.4, min 295.8, max 297.3; whole file: mean 301.6 s
Mouth days (not part of the claim): runs per day {'Dot': {'mean': 5.0, 'min': 5.0, 'max': 5.0}, 'Social': {'mean': 15.0, 'min': 15.0, 'max': 15.0}}
Result: partially reproduced. Monkey L has 10.0 social runs per day (9 to 11) and 6.02 control runs (4 to 8); monkey T has 14 social runs on every Eyes day and 5 control runs. Runs hold 295.8 to 297.3 s of valid samples. The ranges and three of the four means hold; monkey L's mean of 6.02 control runs lies outside 5 ± 0.5 (14 of its 45 days have 8 control runs). Monkey T's 14 social runs on every day are within the tolerance of the reported mean of 15 but never reach the reported 15 (the Mouth days have 15).
fig1e-stim-counts-comparable and figs1d-stim-counts-nonsocial¶
In the live social gaze condition, the total numbers of microstimulations and of shams per day were comparable across the three stimulated areas and between the two monkeys.
In the non-social control condition the numbers of microstimulations and shams per day were also comparable across areas and between monkeys.
def count_tests(cond):
k = trials[trials.kept & (trials.roi == "Eyes") & (trials.cond == cond)]
c = k.groupby(["date", "subject", "region", "type"]).size().unstack("type", fill_value=0).reset_index()
out = []
for typ in ["stim", "sham"]:
tests = []
for subj in [None, "L", "T"]:
x = c if subj is None else c[c.subject == subj]
for a, b in [("OFC", "dmPFC"), ("OFC", "ACCg"), ("dmPFC", "ACCg")]:
tests.append(dict(type=typ, comparison=f"{a} vs {b}, {subj or 'both'}",
p=ranksum(x[x.region == a][typ], x[x.region == b][typ])))
for a in AREAS:
x = c[c.region == a]
tests.append(dict(type=typ, comparison=f"L vs T, {a}",
p=ranksum(x[x.subject == "L"][typ], x[x.subject == "T"][typ])))
t = pd.DataFrame(tests)
t["p_fdr"] = multipletests(t.p, method="fdr_bh")[1]
out.append(t)
summ = c.groupby(["subject", "region"])[["stim", "sham"]].agg(["mean", "min", "max"]).round(1)
return summ, pd.concat(out, ignore_index=True)
summ, t = count_tests("Social")
print("kept trials per day, social runs:")
print(summ)
print(t.round(4).to_string())
for typ in ["stim", "sham"]:
x = t[t.type == typ]
print(f"{typ}: smallest corrected p {x.p_fdr.min():.3f}; all corrected p > 0.90: {bool((x.p_fdr > 0.9).all())}; "
f"smallest uncorrected p {x.p.min():.3f}")
kept trials per day, social runs:
type stim sham
mean min max mean min max
subject region
L ACCg 121.7 93 141 124.5 88 145
OFC 121.0 87 143 123.1 88 146
dmPFC 118.6 94 141 121.1 98 145
T ACCg 115.9 50 161 117.8 54 157
OFC 118.9 80 146 120.7 76 147
dmPFC 118.3 64 170 118.6 68 170
type comparison p p_fdr
0 stim OFC vs dmPFC, both 0.8693 0.9080
1 stim OFC vs ACCg, both 0.7097 0.9080
2 stim dmPFC vs ACCg, both 0.8966 0.9080
3 stim OFC vs dmPFC, L 0.6329 0.9080
4 stim OFC vs ACCg, L 0.6475 0.9080
5 stim dmPFC vs ACCg, L 0.5606 0.9080
6 stim OFC vs dmPFC, T 0.9080 0.9080
7 stim OFC vs ACCg, T 0.7073 0.9080
8 stim dmPFC vs ACCg, T 0.7506 0.9080
9 stim L vs T, OFC 0.7325 0.9080
10 stim L vs T, dmPFC 0.7509 0.9080
11 stim L vs T, ACCg 0.5092 0.9080
12 sham OFC vs dmPFC, both 0.7034 0.9416
13 sham OFC vs ACCg, both 0.7291 0.9416
14 sham dmPFC vs ACCg, both 0.6841 0.9416
15 sham OFC vs dmPFC, L 0.6184 0.9416
16 sham OFC vs ACCg, L 0.8680 0.9416
17 sham dmPFC vs ACCg, L 0.3502 0.9416
18 sham OFC vs dmPFC, T 0.7948 0.9416
19 sham OFC vs ACCg, T 0.8397 0.9416
20 sham dmPFC vs ACCg, T 0.7948 0.9416
21 sham L vs T, OFC 0.8643 0.9416
22 sham L vs T, dmPFC 0.9416 0.9416
23 sham L vs T, ACCg 0.1873 0.9416
stim: smallest corrected p 0.908; all corrected p > 0.90: True; smallest uncorrected p 0.509
sham: smallest corrected p 0.942; all corrected p > 0.90: True; smallest uncorrected p 0.187
Result: reproduced. Of the 24 rank-sum tests on the social runs, the smallest FDR-corrected p is 0.908 for stimulations and 0.942 for shams, all above 0.90 as reported.
Result: reproduced. No corrected p is below 0.05 (smallest 0.108 for stimulations, 0.082 for shams). Not all exceed 0.10 as the paper reports: the sham comparisons between the monkeys for OFC and dmPFC give corrected p = 0.082 (uncorrected 0.014 and 0.010), monkey T having more control trials per day.
figs1e-eyes-more-than-face and figs1e-eyes-vs-rdm-comparable¶
The stimulated monkeys made significantly more fixations per day on the partner's Eyes than on the non-eye Face, for days of all three stimulated areas.
The number of fixations on the partner's Eyes per day was comparable to the number of fixations on the RDM stimulus, for days of all three areas.
Day counts under the claims' readings (a) the day's total and (b) the mean per run. The events
counted: M1's fixations in the ROI, the M1-only events plus the M1 side of mutual events (the
Methods' "total number of gaze fixations ... of the stimulated monkey"), with M1-only events as
the variant (the code's event tables). The non-eye Face is the files' m1_neface table, which includes the Mouth, plus
mutual_face and mutual_mouth for mutual events (mutual_face is the face outside the Eyes and
the Mouth; data check). RDM fixations are the M1 events labelled Eyes in the Dot
runs (the RDM ROI is the Eyes ROI's location and size). FDR over the 15 tests of the claims file.
def roi_counts(evset):
e = events[(events.roi == "Eyes") & (events.who == "m1")]
if evset == "M1-only":
e = e[~e.mutual]
# non-eye Face: m1_neface (Mouth included), and mutual_face and mutual_mouth for mutual events
# (m1_wholeface = m1_eye + m1_neface; mutual_wholeface = mutual_eye + mutual_face + mutual_mouth)
sel = {"Eyes": (e.cond == "Social") & e.table.isin(["m1_eye", "mutual_eye"]),
"neFace": (e.cond == "Social") & e.table.isin(["m1_neface", "mutual_face", "mutual_mouth"]),
"RDM": (e.cond == "Dot") & (e.table == "m1_eye")}
runs_ = meta[meta.roi == "Eyes"]
out = {}
for name, s in sel.items():
cond = "Dot" if name == "RDM" else "Social"
per_run = e[s].groupby("run_id").size().reindex(runs_[runs_.cond == cond].run_id, fill_value=0)
rr = runs_.set_index("run_id").loc[per_run.index]
df = pd.DataFrame({"n": per_run.to_numpy(), "date": rr.date.to_numpy(), "region": rr.region.to_numpy(),
"subject": rr.subject.to_numpy()})
out[name] = df.groupby(["date", "region", "subject"]).n.agg(["sum", "mean"])
return pd.concat(out, axis=1).reset_index()
def roi_tests(c, reading):
col = "sum" if reading == "a" else "mean"
tests = []
for a in AREAS:
x = c[c.region == a]
tests.append(dict(test=f"Eyes vs neFace, {a}", kind="signrank", p=signrank(x[("Eyes", col)], x[("neFace", col)])[1],
median_diff=np.median(x[("Eyes", col)] - x[("neFace", col)])))
tests.append(dict(test=f"Eyes vs RDM, {a}", kind="signrank", p=signrank(x[("Eyes", col)], x[("RDM", col)])[1],
median_diff=np.median(x[("Eyes", col)] - x[("RDM", col)])))
for roi in ["Eyes", "neFace", "RDM"]:
for a, b in [("OFC", "dmPFC"), ("OFC", "ACCg"), ("dmPFC", "ACCg")]:
tests.append(dict(test=f"{roi}, {a} vs {b}", kind="ranksum",
p=ranksum(c[c.region == a][(roi, col)], c[c.region == b][(roi, col)]), median_diff=np.nan))
t = pd.DataFrame(tests)
t["p_fdr"] = multipletests(t.p, method="fdr_bh")[1]
return t
roi_res = {}
for evset in ["M1 and mutual", "M1-only"]:
c = roi_counts(evset)
for reading in ["a", "b"]:
roi_res[(evset, reading)] = roi_tests(c, reading)
c_main = roi_counts("M1 and mutual")
print("fixations per day (reading a, the day's total), M1 and mutual events:")
print(c_main.groupby("region")[[("Eyes", "sum"), ("neFace", "sum"), ("RDM", "sum")]].median())
fixations per day (reading a, the day's total), M1 and mutual events:
Eyes neFace RDM
sum sum sum
region
ACCg 382.0 171.0 180.0
OFC 384.0 179.0 211.0
dmPFC 408.0 134.0 177.0
for k, t in roi_res.items():
x = t[t.test.str.startswith("Eyes vs neFace")]
print(f"events {k[0]}, reading ({k[1]}):")
print(x[["test", "median_diff", "p", "p_fdr"]].to_string(index=False, float_format=lambda v: f"{v:.3g}"))
events M1 and mutual, reading (a):
test median_diff p p_fdr
Eyes vs neFace, OFC 230 6.28e-06 2.5e-05
Eyes vs neFace, dmPFC 258 5.6e-06 2.5e-05
Eyes vs neFace, ACCg 234 5.6e-06 2.5e-05
events M1 and mutual, reading (b):
test median_diff p p_fdr
Eyes vs neFace, OFC 21.8 1.49e-07 1.12e-06
Eyes vs neFace, dmPFC 21.4 1.49e-08 2.24e-07
Eyes vs neFace, ACCg 20.3 5.6e-06 2.8e-05
events M1-only, reading (a):
test median_diff p p_fdr
Eyes vs neFace, OFC 234 6.28e-06 3.14e-05
Eyes vs neFace, dmPFC 233 1.49e-08 1.12e-07
Eyes vs neFace, ACCg 225 1.49e-08 1.12e-07
events M1-only, reading (b):
test median_diff p p_fdr
Eyes vs neFace, OFC 20.4 6.28e-06 3.14e-05
Eyes vs neFace, dmPFC 21.1 5.6e-06 3.14e-05
Eyes vs neFace, ACCg 20.5 5.6e-06 3.14e-05
Result: reproduced. Under both day counts and both event sets, Eyes fixations exceed non-eye Face fixations for every area, with corrected p at most 3.1e-05, below 10^-4 (median day totals of 382 to 408 Eyes against 134 to 179 non-eye Face fixations).
for k, t in roi_res.items():
x = t[t.test.str.startswith("Eyes vs RDM")]
print(f"events {k[0]}, reading ({k[1]}):")
print(x[["test", "median_diff", "p", "p_fdr"]].to_string(index=False, float_format=lambda v: f"{v:.3g}"))
print("\nthe nine rank-sum tests between areas, events M1 and mutual, reading (a):")
print(roi_res[("M1 and mutual", "a")].query("kind == 'ranksum'")[["test", "p", "p_fdr"]].to_string(index=False, float_format=lambda v: f"{v:.3g}"))
events M1 and mutual, reading (a):
test median_diff p p_fdr
Eyes vs RDM, OFC 177 0.000121 0.000363
Eyes vs RDM, dmPFC 202 6.66e-06 2.5e-05
Eyes vs RDM, ACCg 241 0.000162 0.000406
events M1 and mutual, reading (b):
test median_diff p p_fdr
Eyes vs RDM, OFC 6.83 0.17 0.397
Eyes vs RDM, dmPFC 12.3 0.185 0.397
Eyes vs RDM, ACCg 14.1 0.279 0.524
events M1-only, reading (a):
test median_diff p p_fdr
Eyes vs RDM, OFC 157 0.00172 0.00429
Eyes vs RDM, dmPFC 175 5.49e-05 0.000206
Eyes vs RDM, ACCg 208 0.00131 0.00394
events M1-only, reading (b):
test median_diff p p_fdr
Eyes vs RDM, OFC 2.23 0.111 0.303
Eyes vs RDM, dmPFC 7.4 0.141 0.303
Eyes vs RDM, ACCg 11.3 0.185 0.348
the nine rank-sum tests between areas, events M1 and mutual, reading (a):
test p p_fdr
Eyes, OFC vs dmPFC 0.904 0.952
Eyes, OFC vs ACCg 0.755 0.952
Eyes, dmPFC vs ACCg 0.917 0.952
neFace, OFC vs dmPFC 0.0127 0.0273
neFace, OFC vs ACCg 0.333 0.554
neFace, dmPFC vs ACCg 0.246 0.462
RDM, OFC vs dmPFC 0.938 0.952
RDM, OFC vs ACCg 0.762 0.952
RDM, dmPFC vs ACCg 0.952 0.952
Result: not reproduced. Under the paper's reading (a), the day's total, Eyes fixations exceed RDM fixations for every area (corrected p = 0.00036, 0.000025 and 0.00041 for OFC, dmPFC and ACCg; median differences of 177 to 241 fixations per day), since a day holds about twice as many social runs as control runs. Under the code's reading (b), the mean per run, the counts are comparable (corrected p = 0.40, 0.40 and 0.52), as the paper reports. Whether the claims file's choice of reading (a) as the deciding one misreads the paper is asked of a reader (questions/000959.yaml); the outcome is judged on the claims file as it stands.
Social gaze distance: common tables¶
Kept trials of the Eyes days. Distances in pixels; the medians are also given in degrees
(DEG_PER_PX). Fixation readings: methods (M1-only and mutual events) and code (M1-only).
soc = trials[trials.kept & (trials.cond == "Social") & (trials.roi == "Eyes")].copy()
dot = trials[trials.kept & (trials.cond == "Dot") & (trials.roi == "Eyes")].copy()
FIX = {r: ms.m1_fixations(events, r) for r in ["methods", "code"]}
def distance_days(tr, reading, t1=ms.POST_GAZE_S, hemi=None, centre=("eyes_x", "eyes_y")):
tr = tr.copy()
tr["d"] = ms.trial_distance(tr, FIX[reading], centre=centre, t1=t1, hemi=hemi)
return tr, ms.day_means(tr, "d")
def show(tab, label="", degrees=True):
t = tab.copy()
if degrees:
t["median_diff_deg"] = t.median_diff * ms.DEG_PER_PX
print(label)
print(t.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
dist = {}
for reading in ["methods", "code"]:
for hemi in [None, "contra", "ipsi"]:
dist[(reading, hemi)] = distance_days(soc, reading, hemi=hemi)
tr0 = dist[("methods", None)][0]
print("social trials with no fixation in [0, 1.5) s:", int(tr0.d.isna().sum()), "of", len(tr0))
print("trials per day and type, quantiles:", tr0.groupby(["date", "type"]).size().quantile([0, 0.5, 1]).to_dict())
print("trial distance (pixels) range:", round(tr0.d.min(), 1), round(tr0.d.max(), 1))
social trials with no fixation in [0, 1.5) s: 1270 of 19474
trials per day and type, quantiles: {0.0: 50.0, 0.5: 123.5, 1.0: 170.0}
trial distance (pixels) range: 1.7 1494.1
fig2a-density-map-ofc¶
After OFC microstimulations, compared with shams, the stimulated monkey's fixations in the 1.5 s post-gaze epoch were more clustered around the partner monkey.
A 40 by 40 grid of 1-degree squares centred on the run's Eyes centre; monkey T's maps flipped horizontally; per day stimulation and sham maps; normalisation (1) z-score of each map, (2) difference divided by the day's total of stimulation and sham fixations. The measure is the day-averaged difference summed over the squares whose centres lie in (a) the whole Face ROI, (b) the Eyes ROI. Fixations: the Methods reading.
edges = np.arange(-20, 21, 1.0)
cent = (edges[:-1] + edges[1:]) / 2
fx = FIX["methods"].dropna(subset=["x"])
pairs = ms.window_pairs(soc, fx, 0.0, ms.POST_GAZE_S)
pairs["type"] = soc.loc[pairs.trial, "type"].to_numpy()
sgn = np.where(pairs.subject == "T", -1.0, 1.0) # negative dx is ipsilateral for both monkeys
pairs["dx"] = (pairs.x - pairs.eyes_x) * ms.DEG_PER_PX * sgn
pairs["dy"] = (pairs.y - pairs.eyes_y) * ms.DEG_PER_PX
print(f"post-gaze fixations: {len(pairs)}; outside the 40-degree grid: {int(((pairs.dx.abs() >= 20) | (pairs.dy.abs() >= 20)).sum())}")
gx, gy = np.meshgrid(cent, cent, indexing="xy")
dmaps, rows = {}, []
for (date, region, subject), g in pairs.groupby(["date", "region", "subject"]):
H = {t: np.histogram2d(g[g.type == t].dy, g[g.type == t].dx, bins=[edges, edges])[0] for t in ["stim", "sham"]}
z = {t: (H[t] - H[t].mean()) / H[t].std() for t in H}
norm = {"zscore": z["stim"] - z["sham"], "total": (H["stim"] - H["sham"]) / (H["stim"].sum() + H["sham"].sum())}
# squares inside the ROI, relative to the Eyes centre, in every social run of the day
# (three recalibrated days have two ROI sets; the intersection keeps the squares common to both)
masks = {"face": np.ones_like(gx, bool), "eyes": np.ones_like(gx, bool)}
for rid in soc[soc.date == date].run_id.unique():
ro = runs[rid]["rois"]
e = ms._centre(ro["roi_m1_eyes"])
for nm, key in [("face", "roi_m1_face"), ("eyes", "roi_m1_eyes")]:
r = ro[key]
xs = np.sort((r[[0, 2]] - e[0]) * ms.DEG_PER_PX * (-1 if subject == "T" else 1))
ys = np.sort((r[[1, 3]] - e[1]) * ms.DEG_PER_PX)
masks[nm] &= (gx >= xs[0]) & (gx <= xs[1]) & (gy >= ys[0]) & (gy <= ys[1])
for k, m in norm.items():
dmaps.setdefault((region, k), []).append(m)
rows.append(dict(date=date, region=region, subject=subject, norm=k,
face=m[masks["face"]].sum(), eyes=m[masks["eyes"]].sum(),
n_face_squares=masks["face"].sum(), n_eye_squares=masks["eyes"].sum()))
hm = pd.DataFrame(rows)
res = hm.groupby(["region", "norm"])[["face", "eyes"]].mean()
print("day-averaged stimulation minus sham density summed over the squares in each ROI:")
print(res.round(4))
print("by monkey:")
print(hm.groupby(["region", "subject", "norm"])[["face", "eyes"]].mean().round(4))
print("squares inside the ROIs per day, median:", hm.n_face_squares.median(), hm.n_eye_squares.median())
print("recalibrated days and their regions:", meta[meta.date.isin(nvar)].groupby("date").region.first().to_dict())
fig, ax = plt.subplots(1, 3, figsize=(12, 3.6))
for a, region in zip(ax, AREAS):
m = np.mean(dmaps[(region, "zscore")], 0)
im = a.imshow(m, origin="upper", extent=[-20, 20, 20, -20], cmap="RdBu_r", vmin=-0.6, vmax=0.6)
a.set_title(f"{region}, z-scored stim - sham")
print(f"{region}: plotted range {m.min():.2f} to {m.max():.2f}; squares beyond the colour scale {(np.abs(m) > 0.6).sum()}")
a.set_xlabel("deg from the Eyes centre (ipsi <- -> contra)")
plt.colorbar(im, ax=ax)
plt.show()
post-gaze fixations: 53542; outside the 40-degree grid: 10771
day-averaged stimulation minus sham density summed over the squares in each ROI:
face eyes
region norm
ACCg total -0.0048 -0.0013
zscore -2.5489 -0.2872
OFC total 0.0093 0.0046
zscore 7.4696 3.7107
dmPFC total -0.0027 -0.0021
zscore -2.5599 -1.9739
by monkey:
face eyes
region subject norm
ACCg L total -0.0055 -0.0011
zscore -2.9493 0.3962
T total -0.0040 -0.0016
zscore -2.0484 -1.1415
OFC L total 0.0100 0.0072
zscore 7.2513 5.0653
T total 0.0085 0.0013
zscore 7.7426 2.0174
dmPFC L total -0.0003 0.0001
zscore -1.8320 -1.4102
T total -0.0058 -0.0050
zscore -3.4698 -2.6786
squares inside the ROIs per day, median: 42.0 24.0
recalibrated days and their regions: {'08142019': 'dmPFC', '12102019': 'dmPFC', '12132019': 'ACCg'}
OFC: plotted range -1.80 to 1.75; squares beyond the colour scale 29
dmPFC: plotted range -1.59 to 1.33; squares beyond the colour scale 26
ACCg: plotted range -1.36 to 1.44; squares beyond the colour scale 24
Result: reproduced. The day-averaged OFC stimulation minus sham density summed over the whole Face is positive under both normalisations (7.47 z-scored, 0.0093 by total), and so is the sum over the Eyes (3.71 and 0.0046); it is positive in each monkey under each. The dmPFC and ACCg sums over the whole Face are negative under both. 10,771 of 53,542 post-gaze fixations lie outside the 40-degree grid. The three recalibrated days, whose ROI squares are now the intersection over their runs, are dmPFC and ACCg days, so the OFC sums do not depend on them.
fig2c-ofc-social-gaze-distance¶
OFC microstimulations decreased the social gaze distance, the mean distance of the stimulated monkey's fixations in the post-gaze epoch from the centre of the partner's Eyes, compared with shams; the effect held in each monkey.
Result: reproduced. With M1-only and mutual events, OFC stimulation lowered the day-mean social gaze distance (median stimulation minus sham -17.9 pixels, -0.70 degrees; p = 0.00021 over 27 days); monkey L gives p = 0.0067 (paper 0.008) and monkey T p = 0.016 (paper 0.002), both negative. The code's M1-only reading gives p = 0.00013, 0.0067 and 0.0034.
fig2c-dmpfc-accg-no-distance-effect¶
Microstimulations of dmPFC or ACCg did not change social gaze distance.
for reading in ["methods", "code"]:
show(day_table(dist[(reading, None)][1], areas=["dmPFC", "ACCg"]), f"fixations: {reading}")
fixations: methods region subject n median_diff mean_diff p median_diff_deg dmPFC both 27 -6.527 -6.849 0.2109 -0.2568 dmPFC L 15 -7.915 -11.64 0.1354 -0.3115 dmPFC T 12 -3.687 -0.857 0.9097 -0.1451 ACCg both 27 -4.869 -3.433 0.3865 -0.1916 ACCg L 15 -5.39 -12.73 0.08325 -0.2121 ACCg T 12 3.517 8.19 0.5693 0.1384 fixations: code region subject n median_diff mean_diff p median_diff_deg dmPFC both 27 -6.042 -5.587 0.2901 -0.2377 dmPFC L 15 -6.728 -10.19 0.2078 -0.2648 dmPFC T 12 -1.273 0.1655 1 -0.05011 ACCg both 27 -3.908 -4.089 0.2901 -0.1538 ACCg L 15 -5.52 -12.23 0.073 -0.2172 ACCg T 12 3.172 6.081 0.7334 0.1248
Result: reproduced. dmPFC gives p = 0.21 and ACCg p = 0.39 (paper 0.361 and 0.374); the code's reading gives 0.29 and 0.29.
fig2e-ofc-contralateral and fig2e-dmpfc-accg-no-hemifield-effect¶
The OFC effect on social gaze distance was stronger for fixations in the hemifield contralateral to the stimulated hemisphere; it was significant for contralateral fixations, not for ipsilateral ones, and the contralateral effect exceeded the ipsilateral one.
For dmPFC and ACCg there was no stimulation effect on social gaze distance in either hemifield, and no difference between hemifields.
def hemi_table(dsets, areas, subjects=(None, "L", "T")):
rows = []
for a in areas:
for s in subjects:
c = day_test(dsets["contra"], a, s)
i = day_test(dsets["ipsi"], a, s)
dc = dsets["contra"].assign(eff=lambda d: d.stim - d.sham)
di = dsets["ipsi"].assign(eff=lambda d: d.stim - d.sham)
m = dc.merge(di, on=["date", "subject", "region"], suffixes=("_c", "_i"))
m = m[(m.region == a) & ((m.subject == s) if s else True)].dropna(subset=["eff_c", "eff_i"])
n, p = signrank(m.eff_c, m.eff_i)
rows.append(dict(region=a, subject=s or "both", n_contra=c["n"], contra_median=c["median_diff"], contra_p=c["p"],
n_ipsi=i["n"], ipsi_median=i["median_diff"], ipsi_p=i["p"],
n_cvi=n, contra_minus_ipsi_median=np.median(m.eff_c - m.eff_i) if n else np.nan, cvi_p=p))
return pd.DataFrame(rows)
hemi_res = {r: hemi_table({h: dist[(r, h)][1] for h in ["contra", "ipsi"]}, AREAS) for r in ["methods", "code"]}
for r in ["methods", "code"]:
print(f"fixations: {r} (effects in pixels, stim - sham)")
print(hemi_res[r][hemi_res[r].region == "OFC"].to_string(index=False, float_format=lambda v: f"{v:.4g}"))
fixations: methods (effects in pixels, stim - sham) region subject n_contra contra_median contra_p n_ipsi ipsi_median ipsi_p n_cvi contra_minus_ipsi_median cvi_p OFC both 27 -21.21 0.0004787 27 -12.25 0.05847 27 -17.65 0.04626 OFC L 15 -33.52 0.01807 15 -12.16 0.1688 15 -17.65 0.073 OFC T 12 -21.21 0.01221 12 -12.79 0.2036 12 -18.86 0.3013 fixations: code (effects in pixels, stim - sham) region subject n_contra contra_median contra_p n_ipsi ipsi_median ipsi_p n_cvi contra_minus_ipsi_median cvi_p OFC both 27 -30.16 0.0003807 27 -12.25 0.0552 27 -22.44 0.02623 OFC L 15 -33.52 0.02155 15 -12.03 0.2078 15 -19.28 0.08325 OFC T 12 -28.03 0.01221 12 -12.79 0.1763 12 -35.62 0.1763
Result: reproduced. Contralateral: median effect -21.2 pixels, p = 0.00048; ipsilateral: -12.3 pixels, p = 0.058; contralateral against ipsilateral: median -17.7 pixels, p = 0.046 (paper < 0.001, 0.068 and 0.026). Per monkey the contralateral effect gives p = 0.018 (L) and 0.012 (T) (paper 0.015 and 0.012), and the hemifield difference p = 0.073 (L) and 0.30 (T) (paper 0.073 and 0.204). The code's M1-only reading gives the same pattern (0.00038, 0.055 and 0.026).
for r in ["methods", "code"]:
print(f"fixations: {r}")
x = hemi_res[r][(hemi_res[r].region != "OFC") & (hemi_res[r].subject == "both")]
print(x.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
fixations: methods region subject n_contra contra_median contra_p n_ipsi ipsi_median ipsi_p n_cvi contra_minus_ipsi_median cvi_p dmPFC both 27 -4.113 0.4846 27 -10.56 0.1004 27 25.63 0.2198 ACCg both 27 -7.447 0.5303 27 -4.717 0.5303 27 0.9568 0.8223 fixations: code region subject n_contra contra_median contra_p n_ipsi ipsi_median ipsi_p n_cvi contra_minus_ipsi_median cvi_p dmPFC both 27 -1.988 0.5944 27 -8.808 0.1111 27 23.34 0.2386 ACCg both 27 -7.447 0.427 27 -4.717 0.5149 27 6.338 0.9529
Result: reproduced. None of the six tests gives p < 0.05 (dmPFC 0.48, 0.10 and 0.22; ACCg 0.53, 0.53 and 0.82), under either fixation reading.
fig2f-nonsocial-no-distance-effect and fig2gh-nonsocial-no-hemifield-effect¶
In the non-social control condition, microstimulations of none of the three areas changed the distance of fixations from the centre of the RDM stimulus.
In the non-social control condition there was no stimulation effect in either hemifield, and no hemifield difference, for any area.
Dot runs of the Eyes days; the RDM centre is the centre of the run's roi_m1_eyes (the RDM ROI
has the Eyes ROI's location and size); M1 events only (there is no M2 gaze, so no mutual events).
Result: reproduced. No area changes the distance to the RDM centre (OFC p = 0.22, dmPFC 0.64, ACCg 0.50; paper 0.118, 0.719 and 0.302).
Result: reproduced. None of the nine tests over both monkeys gives p < 0.05 (smallest 0.16, OFC contralateral). Per monkey, outside the claim, dmPFC in monkey L gives p = 0.035 for contralateral fixations.
methods-ofc-distance-longer-windows¶
The OFC decrease of social gaze distance also held within 2 s and within 3 s after the trial for both monkeys combined and for monkey L, not for monkey T.
for reading in ["methods", "code"]:
for t1 in [2.0, 3.0]:
show(day_table(distance_days(soc, reading, t1=t1)[1], areas=["OFC"]), f"fixations {reading}, window [0, {t1:g}) s")
fixations methods, window [0, 2) s region subject n median_diff mean_diff p median_diff_deg OFC both 27 -19.68 -21.82 0.002967 -0.7743 OFC L 15 -24.45 -27.86 0.01245 -0.9619 OFC T 12 -5.187 -14.26 0.2334 -0.2041
fixations methods, window [0, 3) s region subject n median_diff mean_diff p median_diff_deg OFC both 27 -15.53 -18 0.003247 -0.611 OFC L 15 -18 -23.38 0.006714 -0.7081 OFC T 12 -10.82 -11.28 0.2036 -0.4258 fixations code, window [0, 2) s region subject n median_diff mean_diff p median_diff_deg OFC both 27 -22.89 -22.99 0.002248 -0.9008 OFC L 15 -23.65 -28.67 0.01025 -0.9305 OFC T 12 -10.2 -15.88 0.1763 -0.4014
fixations code, window [0, 3) s region subject n median_diff mean_diff p median_diff_deg OFC both 27 -16.37 -18.55 0.002469 -0.644 OFC L 15 -18.77 -23.78 0.006714 -0.7385 OFC T 12 -14.67 -12.01 0.1514 -0.5772
Result: reproduced. Within 2 s: both monkeys p = 0.0030, monkey L 0.012, monkey T 0.23; within 3 s: 0.0032, 0.0067 and 0.20 (paper 0.003, 0.008, 0.204 and 0.004, 0.007, 0.204), with stimulation below sham. The code's reading gives the same pattern.
methods-ofc-no-change-in-fixation-number¶
OFC microstimulations did not change the number of fixations within the partner's Eyes or whole Face in the post-gaze window, in either monkey, for windows of 1.5, 2 or 3 s.
def count_or_duration(tr, tables_, t1, what):
ev = events[events.table.isin(tables_) & (events.who == "m1")]
p = ms.window_pairs(tr, ev, 0.0, t1)
p["dur"] = (p.stop - p.start) * 1000.0
agg = p.groupby("trial").size() if what == "count" else p.groupby("trial").dur.sum()
tr = tr.copy()
tr["v"] = tr.index.map(agg).fillna(0).astype(float)
return ms.day_means(tr, "v")
SETS = {"a": {"Eyes": ["m1_eye"], "wholeFace": ["m1_wholeface"]},
"b": {"Eyes": ["m1_eye", "mutual_eye"], "wholeFace": ["m1_wholeface", "mutual_wholeface"]}}
rows = []
for reading in ["a", "b"]:
for roi in ["Eyes", "wholeFace"]:
for t1 in [1.5, 2.0, 3.0]:
dm = count_or_duration(soc, SETS[reading][roi], t1, "count")
for s in ["L", "T", None]:
r = day_test(dm, "OFC", s)
rows.append(dict(reading=reading, roi=roi, window=t1, **r))
fn = pd.DataFrame(rows)
print(fn.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
for reading in ["a", "b"]:
x = fn[(fn.reading == reading) & (fn.subject != "both")]
print(f"reading ({reading}): per-monkey tests with p < 0.05: {int((x.p < 0.05).sum())} of {len(x)}; smallest p {x.p.min():.3f}")
reading roi window region subject n median_diff mean_diff p
a Eyes 1.5 OFC L 15 0.02708 0.02263 0.1876
a Eyes 1.5 OFC T 12 0.02483 0.0195 0.3804
a Eyes 1.5 OFC both 27 0.02609 0.02124 0.08593
a Eyes 2 OFC L 15 0.03056 0.03039 0.2524
a Eyes 2 OFC T 12 0.02817 0.0162 0.4697
a Eyes 2 OFC both 27 0.03024 0.02408 0.1624
a Eyes 3 OFC L 15 0.01502 0.028 0.3303
a Eyes 3 OFC T 12 -0.0007299 -0.03528 0.5693
a Eyes 3 OFC both 27 0.01502 -0.0001291 0.8223
a wholeFace 1.5 OFC L 15 0.04378 0.04633 0.1205
a wholeFace 1.5 OFC T 12 0.0468 0.04041 0.1475
a wholeFace 1.5 OFC both 27 0.04378 0.0437 0.02459
a wholeFace 2 OFC L 15 0.0324 0.03647 0.3591
a wholeFace 2 OFC T 12 0.06302 0.04274 0.2036
a wholeFace 2 OFC both 27 0.04282 0.03926 0.1111
a wholeFace 3 OFC L 15 0.02689 0.04129 0.2524
a wholeFace 3 OFC T 12 -0.008894 -0.01473 0.8501
a wholeFace 3 OFC both 27 -0.003634 0.01639 0.5303
b Eyes 1.5 OFC L 15 0.03356 0.02756 0.2293
b Eyes 1.5 OFC T 12 0.01972 0.01448 0.3804
b Eyes 1.5 OFC both 27 0.02533 0.02175 0.1167
b Eyes 2 OFC L 15 0.02547 0.03604 0.1354
b Eyes 2 OFC T 12 0.02361 0.01014 0.3394
b Eyes 2 OFC both 27 0.02547 0.02453 0.08593
b Eyes 3 OFC L 15 0.04004 0.0365 0.3894
b Eyes 3 OFC T 12 -0.003626 -0.03772 0.4697
b Eyes 3 OFC both 27 0.007137 0.003513 0.9341
b wholeFace 1.5 OFC L 15 0.05585 0.05543 0.02557
b wholeFace 1.5 OFC T 12 0.02053 0.02252 0.6377
b wholeFace 1.5 OFC both 27 0.04106 0.0408 0.06194
b wholeFace 2 OFC L 15 0.02693 0.04824 0.1688
b wholeFace 2 OFC T 12 0.01927 0.02338 0.5693
b wholeFace 2 OFC both 27 0.02693 0.03719 0.1482
b wholeFace 3 OFC L 15 0.01076 0.06131 0.05536
b wholeFace 3 OFC T 12 -0.0665 -0.02581 0.6772
b wholeFace 3 OFC both 27 0.003335 0.02259 0.427
reading (a): per-monkey tests with p < 0.05: 0 of 12; smallest p 0.121
reading (b): per-monkey tests with p < 0.05: 1 of 12; smallest p 0.026
Result: reproduced. Under the code's reading (a), none of the 12 per-monkey tests gives p < 0.05 (smallest 0.12). Outside the claim, both monkeys combined give p = 0.025 for whole Face fixations within 1.5 s (more after stimulation), and under reading (b) monkey L gives p = 0.026 for the whole Face within 1.5 s.
methods-ofc-longer-face-looking¶
With OFC microstimulations the stimulated monkey spent more time looking within the partner's whole Face in the post-gaze epoch (48.9 ms more than with shams, both monkeys combined), and more weakly within the Eyes; dmPFC and ACCg showed no such effect.
for evs in ["a", "b"]:
rows = []
for roi in ["wholeFace", "Eyes"]:
dm = count_or_duration(soc, SETS[evs][roi], ms.POST_GAZE_S, "duration")
for a in AREAS:
for s in [None, "L", "T"]:
rows.append(dict(roi=roi, **day_test(dm, a, s)))
print(f"events ({evs}: {'M1-only, as the code' if evs == 'a' else 'M1-only and mutual'}); total duration in ms, "
"median_diff and mean_diff are readings (a) and (b) of the claim's 'diff'")
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.4g}"))
events (a: M1-only, as the code); total duration in ms, median_diff and mean_diff are readings (a) and (b) of the claim's 'diff'
roi region subject n median_diff mean_diff p
wholeFace OFC both 27 52.8 62.62 0.002469
wholeFace OFC L 15 65.07 85.27 0.006714
wholeFace OFC T 12 46.23 34.31 0.1514
wholeFace dmPFC both 27 -14.68 -8.156 0.804
wholeFace dmPFC L 15 -14.68 -13.72 0.5995
wholeFace dmPFC T 12 -7.692 -1.202 0.9097
wholeFace ACCg both 27 -2.152 -6.794 0.9341
wholeFace ACCg L 15 -6.068 -20.56 0.3028
wholeFace ACCg T 12 26.73 10.41 0.6221
Eyes OFC both 27 12.23 30.31 0.03853
Eyes OFC L 15 18.57 38.63 0.1205
Eyes OFC T 12 7.826 19.9 0.1514
Eyes dmPFC both 27 -11.68 -16.74 0.05847
Eyes dmPFC L 15 -27.47 -20.99 0.04791
Eyes dmPFC T 12 -4.048 -11.43 0.5693
Eyes ACCg both 27 -0.0292 -5.639 0.9906
Eyes ACCg L 15 -26.79 -26.48 0.2293
Eyes ACCg T 12 41.84 20.41 0.1099
events (b: M1-only and mutual); total duration in ms, median_diff and mean_diff are readings (a) and (b) of the claim's 'diff'
roi region subject n median_diff mean_diff p
wholeFace OFC both 27 63.76 54.18 0.0299
wholeFace OFC L 15 64.18 96.72 0.005371
wholeFace OFC T 12 10.54 1.01 0.8501
wholeFace dmPFC both 27 19.72 3.702 0.5303
wholeFace dmPFC L 15 32.66 10.26 0.4887
wholeFace dmPFC T 12 -9.983 -4.494 0.9697
wholeFace ACCg both 27 -10.94 -28.29 0.1624
wholeFace ACCg L 15 -10.94 -35.3 0.2078
wholeFace ACCg T 12 -13.98 -19.53 0.3804
Eyes OFC both 27 11.01 29.27 0.1937
Eyes OFC L 15 13.19 48.42 0.1205
Eyes OFC T 12 -2.477 5.337 0.9097
Eyes dmPFC both 27 -6.176 -14.67 0.2483
Eyes dmPFC L 15 -10.11 -14.24 0.2524
Eyes dmPFC T 12 3.907 -15.22 0.6772
Eyes ACCg both 27 -11.51 -18.49 0.2483
Eyes ACCg L 15 -11.51 -33.83 0.1354
Eyes ACCg T 12 -14.53 0.6773 0.9697
Result: partially reproduced. With the code's M1-only events, OFC whole Face looking is longer after stimulation (p = 0.0025; L 0.0067; T 0.15), the OFC Eyes test gives p = 0.039, and dmPFC and ACCg give p = 0.80 and 0.93, as the paper reports. Under reading (a) of the difference, the median over days, OFC whole Face is 52.8 ms (expected 48.9 ± 4.9) and monkey L 65.1 ms (66.6 ± 6.7), within tolerance, but monkey T is 46.2 ms (39.5 ± 4.0); reading (b), the mean, gives 62.6, 85.3 and 34.3 ms, none within tolerance. dmPFC's median, -14.7 ms, matches -14.5; ACCg's, -2.2 ms, does not match 3.4, nor does OFC Eyes, 12.2 ms against 16.5.
methods-mouth-stimulation-no-effect¶
When microstimulations were triggered by looking at the partner's Mouth (or at an RDM stimulus at the Mouth location), none of the three areas changed social or non-social gaze distance, in either hemifield.
Mouth days (monkey T, 5 per area). Social distance to the centre of the run's roi_m1_mouth;
non-social distance to the centre of the Dot run's roi_m1_mouth (the RDM at the Mouth
location). With 5 days the smallest two-sided exact signed-rank p is 0.0625.
msoc = trials[trials.kept & (trials.cond == "Social") & (trials.roi == "Mouth")]
mdot = trials[trials.kept & (trials.cond == "Dot") & (trials.roi == "Mouth")]
print("Mouth days: social trials", len(msoc), "control trials", len(mdot))
rows = []
mouth_dm = {}
for cond, tr, rd in [("social", msoc, "methods"), ("social, code reading", msoc, "code"), ("non-social", mdot, "code")]:
for h in [None, "contra", "ipsi"]:
dm = distance_days(tr, rd, hemi=h, centre=("mouth_x", "mouth_y"))[1]
mouth_dm[(cond, h)] = dm
for a in AREAS:
rows.append(dict(condition=cond, hemifield=h or "both", **day_test(dm, a, "T")))
mt = pd.DataFrame(rows)
print(mt.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
print("tests with p < 0.05:", int((mt.p < 0.05).sum()), "of", len(mt))
Mouth days: social trials 2444 control trials 1953
condition hemifield region subject n median_diff mean_diff p
social both OFC T 5 -0.2414 -1.878 1
social both dmPFC T 5 -5.766 2.848 1
social both ACCg T 5 -28.62 -24.82 0.1875
social contra OFC T 5 -6.264 -4.404 0.4375
social contra dmPFC T 5 -3.651 -6.338 0.8125
social contra ACCg T 5 -49.84 -33.84 0.1875
social ipsi OFC T 5 -8.384 -4.695 1
social ipsi dmPFC T 5 29.43 32.05 0.1875
social ipsi ACCg T 5 12.19 -5.026 0.8125
social, code reading both OFC T 5 -3 0.4998 1
social, code reading both dmPFC T 5 0.2902 2.452 1
social, code reading both ACCg T 5 -25.83 -22.01 0.1875
social, code reading contra OFC T 5 -1.583 -2.413 0.8125
social, code reading contra dmPFC T 5 2.66 -5.188 1
social, code reading contra ACCg T 5 -31.3 -29.73 0.3125
social, code reading ipsi OFC T 5 -1.856 -3.389 1
social, code reading ipsi dmPFC T 5 29.43 30.05 0.1875
social, code reading ipsi ACCg T 5 12.19 -4.323 0.625
non-social both OFC T 5 4.178 -2.567 1
non-social both dmPFC T 5 -42.23 -32.66 0.125
non-social both ACCg T 5 -16.18 -14.16 0.1875
non-social contra OFC T 5 -5.251 -4.743 1
non-social contra dmPFC T 5 -49.7 -25.65 0.3125
non-social contra ACCg T 5 -19.93 -9.602 0.625
non-social ipsi OFC T 5 -4.918 -25.38 1
non-social ipsi dmPFC T 5 -60.75 -70.78 0.0625
non-social ipsi ACCg T 5 -3.17 -5.299 0.8125
tests with p < 0.05: 0 of 27
Result: reproduced. None of the 27 tests on the Mouth days (5 days of monkey T per area) gives p < 0.05: social OFC, dmPFC and ACCg give p = 1, 1 and 0.19 (paper 1, 1 and 0.125), non-social 1, 0.125 and 0.19 (paper 1, 0.188 and 0.188). With 5 days the smallest possible two-sided p is 0.0625, so these tests could not find an effect at 0.05.
methods-eyes-vs-mouth-bootstrap¶
For contralateral fixations, the OFC stimulation effect on social gaze distance was significantly more negative in Eyes stimulation than in Mouth stimulation, but not for ipsilateral fixations or both hemifields combined; the difference was stronger using monkey T's Eyes days only.
Day effects are stimulation minus sham day means, under the Methods fixation reading (M1-only and mutual events) and, as a variant, the code's (M1-only events). 1,000 draws of 5 days with replacement; p readings (a) fraction of bootstrap medians at or above the Mouth median, (b) twice the smaller tail. Seed 959 (a standard choice).
rows = []
for reading, mkey in [("methods", "social"), ("code", "social, code reading")]:
rng = np.random.default_rng(959)
for h in [None, "contra", "ipsi"]:
e = dist[(reading, h)][1]
e = e[e.region == "OFC"].assign(eff=lambda d: d.stim - d.sham).dropna(subset=["eff"])
m = mouth_dm[(mkey, h)]
m = m[m.region == "OFC"].assign(eff=lambda d: d.stim - d.sham).dropna(subset=["eff"])
mouth_med = np.median(m.eff)
for pool in ["all", "T"]:
x = (e if pool == "all" else e[e.subject == "T"]).eff.to_numpy()
boots = np.median(rng.choice(x, size=(1000, 5), replace=True), axis=1)
pa = np.mean(boots >= mouth_med)
pb = min(1.0, 2 * min(np.mean(boots >= mouth_med), np.mean(boots <= mouth_med)))
rows.append(dict(fixations=reading, hemifield=h or "both", eyes_days=pool, n_eyes=len(x), n_mouth=len(m),
eyes_median=np.median(x), boot_median=np.median(boots), mouth_median=mouth_med,
p_a_one_sided=pa, p_b_two_sided=pb))
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.4g}"))
print("\ncontralateral, all Eyes days, p under readings (a) and (b) for five other seeds:")
for reading, mkey in [("methods", "social"), ("code", "social, code reading")]:
e = dist[(reading, "contra")][1]
e = e[e.region == "OFC"].assign(eff=lambda d: d.stim - d.sham).dropna(subset=["eff"]).eff.to_numpy()
m = mouth_dm[(mkey, "contra")]
mm = np.median((m[m.region == "OFC"].stim - m[m.region == "OFC"].sham).dropna())
for seed in range(1, 6):
b = np.median(np.random.default_rng(seed).choice(e, size=(1000, 5), replace=True), axis=1)
print(f" {reading} fixations, seed {seed}: (a) {np.mean(b >= mm):.3f} (b) {min(1, 2 * min(np.mean(b >= mm), np.mean(b <= mm))):.3f}")
fixations hemifield eyes_days n_eyes n_mouth eyes_median boot_median mouth_median p_a_one_sided p_b_two_sided
methods both all 27 5 -17.87 -17.87 -0.2414 0.043 0.086
methods both T 12 5 -14.02 -13.47 -0.2414 0.037 0.074
methods contra all 27 5 -21.21 -21.21 -6.264 0.078 0.156
methods contra T 12 5 -21.21 -21.21 -6.264 0.037 0.074
methods ipsi all 27 5 -12.25 -12.25 -8.384 0.295 0.59
methods ipsi T 12 5 -12.79 -12.25 -8.384 0.348 0.696
code both all 27 5 -20.38 -20.38 -3 0.05 0.1
code both T 12 5 -8.799 -8.799 -3 0.036 0.072
code contra all 27 5 -30.16 -30.16 -1.583 0.049 0.098
code contra T 12 5 -28.03 -25.89 -1.583 0.006 0.012
code ipsi all 27 5 -12.25 -12.25 -1.856 0.295 0.59
code ipsi T 12 5 -12.79 -12.25 -1.856 0.348 0.696
contralateral, all Eyes days, p under readings (a) and (b) for five other seeds:
methods fixations, seed 1: (a) 0.064 (b) 0.128
methods fixations, seed 2: (a) 0.069 (b) 0.138
methods fixations, seed 3: (a) 0.059 (b) 0.118
methods fixations, seed 4: (a) 0.070 (b) 0.140
methods fixations, seed 5: (a) 0.070 (b) 0.140
code fixations, seed 1: (a) 0.046 (b) 0.092
code fixations, seed 2: (a) 0.040 (b) 0.080
code fixations, seed 3: (a) 0.040 (b) 0.080
code fixations, seed 4: (a) 0.049 (b) 0.098
code fixations, seed 5: (a) 0.043 (b) 0.086
Result: not reproduced. With the Methods fixations, the bootstrap medians of the OFC Eyes effects for contralateral fixations are not significantly below the Mouth median using all Eyes days (p = 0.078 one-sided, 0.156 two-sided; 0.059 to 0.070 one-sided over five other seeds), while for both hemifields combined the one-sided p is 0.043, where the paper finds no difference. With monkey T's Eyes days the contralateral p is 0.037 one-sided (0.074 two-sided) and the combined p is also 0.037. The ipsilateral p values are 0.30 and 0.35. With the code's M1-only fixations, the contralateral p is 0.049 one-sided (0.040 to 0.049 over five other seeds) and 0.006 with monkey T's days, closer to the paper's 0.009 and 0.005, but both hemifields combined gives 0.050 (all days) and 0.036 (monkey T), and the paper's pattern needs the combined p not below 0.05 for both. Neither reading, nor either fixation set, gives the paper's pattern.
Inter-looking interval¶
Per trial, the time from the trial to the start of the first M1 whole Face event (M1-only or mutual, Methods, Gaze latency analyses) that starts within 5 s after it; NaN when there is none.
wf_m1 = events[(events.who == "m1") & events.table.isin(["m1_wholeface", "mutual_wholeface"])]
wf_m2 = events[(events.who == "m2") & events.table.isin(["m2_wholeface", "mutual_wholeface"])]
soc["ili"] = ms.first_latency(soc, wf_m1)
dot["ili"] = ms.first_latency(dot, events[(events.who == "m1") & (events.table == "m1_wholeface")])
ili_dm = ms.day_means(soc, "ili")
print("social trials with an inter-looking interval:", int(soc.ili.notna().sum()), "of", len(soc),
"; range (s):", round(soc.ili.min(), 3), round(soc.ili.max(), 3),
"; intervals of 0 (an event starting at the trial sample):", int((soc.ili == 0).sum()))
social trials with an inter-looking interval: 11589 of 19474 ; range (s): 0.001 4.999 ; intervals of 0 (an event starting at the trial sample): 0
fig3b-ofc-inter-looking-interval¶
OFC microstimulations decreased the inter-looking interval, the latency of the stimulated monkey to look back at the partner's whole Face within 5 s after the trial, compared with shams; dmPFC and ACCg microstimulations did not.
show(day_table(ili_dm), "inter-looking interval (s), day level", degrees=False)
v = soc.copy()
v["ili_m1only"] = ms.first_latency(v, events[(events.who == "m1") & (events.table == "m1_wholeface")])
show(day_table(ms.day_means(v, "ili_m1only")), "variant: M1-only whole Face events", degrees=False)
p0 = ms.window_pairs(v, wf_m1, 0.0, ms.LATENCY_WINDOW_S).groupby("trial").lag.min()
v["ili_with0"] = v.index.map(p0).astype(float)
print("trials whose first whole Face event starts at the trial sample itself:", int((v.ili_with0 == 0).sum()))
show(day_table(ms.day_means(v, "ili_with0"), subjects=(None,)), "variant: events starting at the trial sample counted (latency 0)", degrees=False)
inter-looking interval (s), day level region subject n median_diff mean_diff p OFC both 27 -0.06928 -0.1072 0.04356 OFC L 15 -0.06928 -0.09927 0.2524 OFC T 12 -0.07442 -0.1172 0.1099 dmPFC both 27 0.02292 -0.01849 0.7676 dmPFC L 15 -0.002927 -0.03938 0.8904 dmPFC T 12 0.1168 0.007634 0.6772 ACCg both 27 0.02053 0.04014 0.3361 ACCg L 15 -0.003057 0.003093 0.7615 ACCg T 12 0.02241 0.08646 0.2334
variant: M1-only whole Face events region subject n median_diff mean_diff p OFC both 27 -0.09935 -0.1297 0.002248 OFC L 15 -0.07347 -0.08706 0.1876 OFC T 12 -0.2362 -0.1831 0.002441 dmPFC both 27 0.05669 0.005792 0.4996 dmPFC L 15 0.05669 -0.01896 0.7197 dmPFC T 12 0.07132 0.03673 0.5186 ACCg both 27 0.05175 0.05235 0.2687 ACCg L 15 0.1006 0.04511 0.4212 ACCg T 12 0.0005409 0.0614 0.5186 trials whose first whole Face event starts at the trial sample itself: 91 variant: events starting at the trial sample counted (latency 0) region subject n median_diff mean_diff p OFC both 27 -0.03758 -0.1089 0.05847 dmPFC both 27 -0.0008149 -0.03068 0.9341 ACCg both 27 0.006374 0.04043 0.3865
Result: reproduced. OFC stimulation shortened the inter-looking interval (median day difference -0.069 s, p = 0.044; paper 0.035), while dmPFC (p = 0.77) and ACCg (p = 0.34) did not (paper 0.792 and 0.291). Neither monkey alone is significant (L 0.25, T 0.11; paper 0.188 and 0.092). The result rests on near-threshold details: counting an event that starts at the trial sample as a latency of 0 gives OFC p = 0.058, and M1-only events give p = 0.0022.
ili-nonsocial-no-effect¶
The inter-looking interval did not change with stimulation in the non-social control condition, for any area.
Result: reproduced. No area changes the interval in the control runs (OFC p = 0.70, dmPFC 0.13, ACCg 0.84; paper 0.773, 0.080 and 0.943).
ili-missing-trials-41¶
On 41% of microstimulation and sham trials per day on average, the stimulated monkey did not look back at the partner's whole Face within the 5 s window.
miss = soc.assign(m=soc.ili.isna()).groupby(["date", "region", "subject"]).m.mean().reset_index()
print(f"(a) all {miss.date.nunique()} Eyes days: {100 * miss.m.mean():.1f}%")
print(f"(b) the {miss[miss.region == 'OFC'].date.nunique()} OFC days: {100 * miss[miss.region == 'OFC'].m.mean():.1f}%")
print("by area and monkey (%):")
print((100 * miss.groupby(["region", "subject"]).m.mean()).round(1))
print(f"pooled over trials: {100 * soc.ili.isna().mean():.1f}%")
(a) all 81 Eyes days: 41.3%
(b) the 27 OFC days: 39.9%
by area and monkey (%):
region subject
ACCg L 38.8
T 45.2
OFC L 38.0
T 42.2
dmPFC L 41.9
T 43.2
Name: m, dtype: float64
pooled over trials: 40.5%
Result: reproduced. Averaged over the 81 Eyes days, 41.3% of trials have no whole Face event within 5 s (reading a; expected 41 ± 4.1); over the 27 OFC days, 39.9% (reading b).
figs3b-ili-trial-level¶
With all trials of each area pooled across days, OFC microstimulations still decreased the inter-looking interval, while dmPFC and ACCg did not.
def trial_level(tr, col):
rows = []
for a in AREAS:
for s in [None, "L", "T"]:
x = tr[(tr.region == a) & ((tr.subject == s) if s else True)]
st, sh = x[x.type == "stim"][col].dropna(), x[x.type == "sham"][col].dropna()
rows.append(dict(region=a, subject=s or "both", n_stim=len(st), n_sham=len(sh),
median_stim=st.median(), median_sham=sh.median(),
mean_stim=st.mean(), mean_sham=sh.mean(), p=ranksum(st, sh)))
return pd.DataFrame(rows)
print(trial_level(soc, "ili").to_string(index=False, float_format=lambda v: f"{v:.4g}"))
region subject n_stim n_sham median_stim median_sham mean_stim mean_sham p OFC both 1992 1971 1.193 1.359 1.589 1.694 0.01665 OFC L 1164 1131 1.342 1.525 1.718 1.83 0.04909 OFC T 828 840 0.9435 1.115 1.408 1.51 0.1359 dmPFC both 1846 1918 1.34 1.361 1.688 1.729 0.3933 dmPFC L 1028 1090 1.559 1.664 1.844 1.904 0.2464 dmPFC T 818 828 1.099 0.9683 1.493 1.498 0.775 ACCg both 1910 1952 1.346 1.344 1.715 1.691 0.6215 ACCg L 1112 1163 1.563 1.546 1.846 1.834 0.9494 ACCg T 798 789 1.091 0.9489 1.533 1.48 0.304
Result: reproduced. Pooling trials, OFC stimulation shortened the interval (median 1.19 against 1.36 s, p = 0.017; paper 0.010), while dmPFC (p = 0.39) and ACCg (p = 0.62) did not (paper 0.301 and 0.602).
Reciprocation latency¶
Whole Face events of M1 and of M2 that start within 5 s after the trial, mutual events in both lists; for each run of consecutive M2 events followed by an M1 event, the time from the run's first M2 event to that M1 event; the trial's value is the mean over sequences.
soc["recip"] = ms.reciprocation(soc, wf_m1, wf_m2, first="m2")
print("social trials with a reciprocation latency:", int(soc.recip.notna().sum()), "of", len(soc),
"; range (s):", round(soc.recip.min(), 3), round(soc.recip.max(), 3))
social trials with a reciprocation latency: 2689 of 19474 ; range (s): 0.001 4.894
figs3c-reciprocation-day-level-none¶
At the day level, microstimulations did not significantly reduce the reciprocation latency in any area; on 86% of trials per day on average there was no sequence of the partner looking at the stimulated monkey followed by the stimulated monkey looking back.
show(day_table(ms.day_means(soc, "recip")), "reciprocation latency (s), day level", degrees=False)
rm = soc.assign(m=soc.recip.isna()).groupby(["date", "region"]).m.mean().reset_index()
print(f"\ntrials with no sequence: (a) all {rm.date.nunique()} days {100 * rm.m.mean():.1f}%; "
f"(b) OFC days {100 * rm[rm.region == 'OFC'].m.mean():.1f}%")
print("days with no sequence in a trial type, by area:",
ms.day_means(soc, "recip").groupby("region")[["stim", "sham"]].apply(lambda d: d.isna().any(axis=1).sum()).to_dict())
reciprocation latency (s), day level
region subject n median_diff mean_diff p
OFC both 27 -0.08724 -0.1452 0.1111
OFC L 15 -0.1051 -0.1528 0.1876
OFC T 12 -7.689e-05 -0.1357 0.4697
dmPFC both 27 -0.1042 -0.04121 0.5619
dmPFC L 15 0.06931 0.04912 0.5995
dmPFC T 12 -0.1925 -0.1541 0.07715
ACCg both 27 -0.07884 -0.06129 0.3242
ACCg L 15 0.002315 0.05722 0.804
ACCg T 12 -0.1043 -0.2094 0.06396
trials with no sequence: (a) all 81 days 86.4%; (b) OFC days 85.1%
days with no sequence in a trial type, by area: {'ACCg': 0, 'OFC': 0, 'dmPFC': 0}
Result: reproduced. No area changes the day-level reciprocation latency (OFC p = 0.11, dmPFC 0.56, ACCg 0.32; paper 0.130, 0.701 and 0.400), and 86.4% of trials per day have no M2-then-M1 sequence (reading a; expected 86 ± 8.6; 85.1% on the OFC days).
fig3d-ofc-reciprocation-latency¶
With trials pooled across days, OFC microstimulations decreased the latency of the stimulated monkey to reciprocate the partner's gaze to his face; dmPFC and ACCg did not.
print(trial_level(soc, "recip").to_string(index=False, float_format=lambda v: f"{v:.4g}"))
region subject n_stim n_sham median_stim median_sham mean_stim mean_sham p OFC both 482 503 0.8385 1.009 1.191 1.338 0.02277 OFC L 254 250 0.7994 0.894 1.183 1.382 0.1457 OFC T 228 253 0.9099 1.124 1.2 1.296 0.07758 dmPFC both 410 434 0.9561 0.98 1.287 1.31 0.7205 dmPFC L 228 254 1.107 0.9109 1.401 1.312 0.3109 dmPFC T 182 180 0.8586 1.013 1.145 1.308 0.08648 ACCg both 421 439 0.9856 1.027 1.3 1.311 0.3622 ACCg L 255 257 1.031 1.005 1.352 1.254 0.7924 ACCg T 166 182 0.7911 1.06 1.221 1.391 0.1056
Result: reproduced. Pooling trials, OFC stimulation shortened the reciprocation latency (median 0.84 against 1.01 s, p = 0.023; paper 0.011), while dmPFC (p = 0.72) and ACCg (p = 0.36) did not (paper 0.777 and 0.368); per monkey, L gives p = 0.15 and T p = 0.078 (paper 0.074 and 0.079).
partner-monkey-unchanged¶
Microstimulations in none of the areas changed the partner monkey's social gaze distance, inter-looking interval or reciprocation latency.
M2 distance: M2's fixations starting in [0, 1.5) s, distance to the centre of roi_m2_eyes
(M1's Eyes in M2's coordinates). M2 inter-looking interval, readings (a) start of the second
minus start of the first M2 whole Face event in the 5-s window, (b) start of the second minus
end of the first. M2 reciprocation latency, readings (a) from the start of the first M1 event of
a run of M1 whole Face events to the next M2 event (the mirror of M1's measure), (b) from the
trial time to the first M2 event that follows an M1 event.
rows = []
for rd in ["methods", "code"]:
t = soc.copy()
t["d2"] = ms.trial_distance(t, ms.m2_fixations(events, rd).dropna(subset=["x"]), centre=("m2_eyes_x", "m2_eyes_y"))
for a in AREAS:
rows.append(dict(measure=f"M2 distance ({rd})", **day_test(ms.day_means(t, "d2"), a)))
p2 = ms.window_pairs(soc, wf_m2, 0.0, ms.LATENCY_WINDOW_S).sort_values(["trial", "start"])
p2["k"] = p2.groupby("trial").cumcount()
first2 = p2[p2.k == 0].set_index("trial")
second2 = p2[p2.k == 1].set_index("trial")
t = soc.copy()
t["ili2a"] = t.index.map(second2.start - first2.start.reindex(second2.index)).astype(float)
t["ili2b"] = t.index.map(second2.start - first2.stop.reindex(second2.index)).astype(float)
t["rec2a"] = ms.reciprocation(soc, wf_m1, wf_m2, first="m1")
p1 = ms.window_pairs(soc, wf_m1, 0.0, ms.LATENCY_WINDOW_S).groupby("trial").start.min()
q = p2.assign(f1=p2.trial.map(p1)).query("start > f1").groupby("trial").lag.min()
t["rec2b"] = t.index.map(q).astype(float)
for col, lab in [("ili2a", "M2 inter-looking interval (a)"), ("ili2b", "M2 inter-looking interval (b)"),
("rec2a", "M2 reciprocation latency (a)"), ("rec2b", "M2 reciprocation latency (b)")]:
for a in AREAS:
rows.append(dict(measure=lab, **day_test(ms.day_means(t, col), a)))
pm = pd.DataFrame(rows)
print(pm.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
print("trials with a value:", {c: int(t[c].notna().sum()) for c in ["ili2a", "ili2b", "rec2a", "rec2b"]}, "of", len(t))
print("tests with p < 0.05:", pm[pm.p < 0.05][["measure", "region", "p"]].to_dict("records"))
measure region subject n median_diff mean_diff p
M2 distance (methods) OFC both 27 8.801 11.66 0.1225
M2 distance (methods) dmPFC both 27 -14.07 -5.961 0.4698
M2 distance (methods) ACCg both 27 15.15 14.06 0.0619
M2 distance (code) OFC both 27 10.26 13.22 0.1111
M2 distance (code) dmPFC both 27 -13.73 -5.918 0.3735
M2 distance (code) ACCg both 27 11.75 14.63 0.06923
M2 inter-looking interval (a) OFC both 27 0.02486 -0.004503 0.9906
M2 inter-looking interval (a) dmPFC both 27 0.02619 -0.06115 0.9153
M2 inter-looking interval (a) ACCg both 27 -0.08403 -0.08026 0.1937
M2 inter-looking interval (b) OFC both 27 0.05699 0.0109 0.9906
M2 inter-looking interval (b) dmPFC both 27 0.01357 -0.0547 0.9153
M2 inter-looking interval (b) ACCg both 27 -0.108 -0.1523 0.3361
M2 reciprocation latency (a) OFC both 27 0.06506 0.01443 0.5619
M2 reciprocation latency (a) dmPFC both 27 0.1144 0.1287 0.04356
M2 reciprocation latency (a) ACCg both 27 -0.02447 0.01002 0.9906
M2 reciprocation latency (b) OFC both 27 0.1052 0.09665 0.4553
M2 reciprocation latency (b) dmPFC both 27 -0.141 0.0258 0.9529
M2 reciprocation latency (b) ACCg both 27 0.006837 0.03378 0.9153
trials with a value: {'ili2a': 2104, 'ili2b': 2104, 'rec2a': 2798, 'rec2b': 2797} of 19474
tests with p < 0.05: [{'measure': 'M2 reciprocation latency (a)', 'region': 'dmPFC', 'p': 0.04355788230895996}]
Result: partially reproduced. M2's social gaze distance (OFC p = 0.12, dmPFC 0.47, ACCg 0.062; paper 0.118, 0.302 and 0.064) and M2's inter-looking interval under both readings (smallest p 0.19) do not change. M2's reciprocation latency does not change under reading (b) (smallest p 0.46), but under reading (a), from the start of the M1 look, dmPFC gives p = 0.044 (median 0.11 s longer after stimulation). The criterion asks that no test give p < 0.05 under every reading; 1 of the 18 tests does (6 distance tests and 3 for each of the four other measures).
epoch-days-76-of-81¶
76 of the 81 days had at least 90 microstimulations and entered the epoch analysis; the five excluded were 1 OFC day of monkey L and, of monkey T, 1 OFC, 2 dmPFC and 1 ACCg day.
allsoc = trials[(trials.cond == "Social") & (trials.roi == "Eyes")]
ns = pd.DataFrame({
"a_kept": allsoc[allsoc.kept & (allsoc.type == "stim")].groupby("date").size(),
"b_all": allsoc[allsoc.type == "stim"].groupby("date").size()}).fillna(0).astype(int)
ns = ns.join(meta.drop_duplicates("date").set_index("date")[["subject", "region"]])
for col in ["a_kept", "b_all"]:
ok = ns[col] >= 90
print(f"reading ({col}): days with at least 90 stimulations {int(ok.sum())} of {len(ns)}; excluded:")
print(ns[~ok].sort_values(["subject", "region"]).to_string())
named = ["08122019", "12222021", "11052021", "12172021", "11042021"]
print("\nthe five days the code names:")
print(ns.loc[named].to_string())
print("stimulations per day, quantiles (a):", ns.a_kept.quantile([0, 0.1, 0.5, 1]).to_dict())
reading (a_kept): days with at least 90 stimulations 76 of 81; excluded:
a_kept b_all subject region
date
08122019 87 96 L OFC
11042021 50 60 T ACCg
12222021 80 86 T OFC
11052021 64 75 T dmPFC
12172021 67 76 T dmPFC
reading (b_all): days with at least 90 stimulations 77 of 81; excluded:
a_kept b_all subject region
date
11042021 50 60 T ACCg
12222021 80 86 T OFC
11052021 64 75 T dmPFC
12172021 67 76 T dmPFC
the five days the code names:
a_kept b_all subject region
date
08122019 87 96 L OFC
12222021 80 86 T OFC
11052021 64 75 T dmPFC
12172021 67 76 T dmPFC
11042021 50 60 T ACCg
stimulations per day, quantiles (a): {0.0: 50.0, 0.1: 97.0, 0.5: 121.0, 1.0: 170.0}
Result: reproduced. Counting kept stimulations (reading a), 76 of the 81 days have at least 90, and the five others are the days the paper's code names (the paper gives only their composition): 08122019 (L, OFC, 87), 12222021 (T, OFC, 80), 11052021 and 12172021 (T, dmPFC, 64 and 67) and 11042021 (T, ACCg, 50). Counting all stimulations delivered (reading b) gives 77 days, the L OFC day having 96.
Causal decomposition claims (blocked)¶
figs4a-causal-strength-no-stim-effect, figs4b-causal-strength-no-epoch-effect,
ofc-slope-unchanged-by-epoch, fig4c-dmpfc-slope-late-greater, fig4d-dmpfc-slope-permutation and
dmpfc-slope-positive all rest on the relative causal strength of each trial, from the causal
decomposition of Yang, Peng and Huang (2018), which is built on ensemble empirical mode
decomposition (EEMD). The Yang et al. code (causal_decomposition.m, release
causal_decomposition_20180901 at github.com/accyang/causal-decomposition-analysis) normalises both
series, runs one EEMD of each, then for each IMF removes it from each series and runs one EEMD of
each again: 2 + 2 x (number of IMFs) EEMDs per trial. Its eemd.m sifts a fixed 10 times per IMF
(maxSift = 10, the default). Neither the noise level nor the ensemble size has a default in
that code; its example call uses noise 0.35 and 1,000 ensembles. The claims file's reading (a)
takes noise 0.2 and 1,000 ensembles, and reading (b), the code's settings, is taken here as the
example's 0.35 and 1,000, so both readings need 1,000 ensembles and cost the same (the noise
level does not change the work).
The cell below times one EEMD of one 150-bin series with PyEMD set to the Yang code's sifting (10 fixed sifts), at 1,000 ensembles on one core and on four processes, and counts the trials each claim needs: every social Eyes trial for figs4a and figs4b, and the first 90 stimulation trials of each of the 76 epoch days for the four slope claims (by area). It also gives a floor for any implementation: the number of cubic-spline envelopes the decomposition fits, times the measured cost of one scipy spline fit on 150 points.
from PyEMD import EEMD, EMD
rng_c = np.random.default_rng(4)
sig = np.cumsum(rng_c.normal(size=150))
sig = (sig - sig.mean()) / sig.std()
cost = {}
for label, kw in [("one core", dict(parallel=False)), ("four processes", dict(parallel=True, processes=4))]:
eemd = EEMD(trials=1000, noise_width=0.2, ext_EMD=EMD(FIXE=10), **kw)
eemd.noise_seed(4)
t0 = time.time()
imfs = eemd.eemd(sig)
cost[label] = time.time() - t0
n_imf = imfs.shape[0]
per_eemd = cost["one core"]
n_eemd = 2 + 2 * n_imf # Yang et al. (2018), causal_decomposition.m
per_trial = per_eemd * n_eemd
print(f"EEMD of one 150-bin series, 1,000 ensembles, 10 fixed sifts: {cost['one core']:.1f} s on one core, "
f"{cost['four processes']:.1f} s on four processes; {n_imf} IMFs; {n_eemd} EEMDs per trial; "
f"{per_trial:.0f} CPU s per trial")
from scipy.interpolate import CubicSpline
xk, tt = np.linspace(0, 149, 25), np.arange(150.0)
t0 = time.time()
for _ in range(1000):
CubicSpline(xk, rng_c.normal(size=25))(tt)
spline_s = (time.time() - t0) / 1000
n_splines = 1000 * (n_imf - 1) * 10 * 2 * n_eemd # ensembles x IMFs x sifts x two envelopes x EEMDs
print(f"one spline fit and evaluation: {spline_s * 1e6:.0f} us; splines per trial {n_splines:,}; "
f"floor per trial {n_splines * spline_s:.0f} s at this cost, {n_splines * 10e-6:.0f} s at 10 us per spline")
epoch_days = ns[ns.a_kept >= 90]
n90 = epoch_days.groupby("region").size() * 90
need = {"figs4a, figs4b: every social Eyes trial": len(soc),
"slope claims: first 90 stimulations of the 76 epoch days": int(n90.sum())}
need.update({f" of which {r} days ({int(n90[r] / 90)})": int(n90[r]) for r in AREAS})
speed4 = cost["one core"] / cost["four processes"]
rows = [dict(subset=k, trials=v, cpu_hours=v * per_trial / 3600, wall_hours_4_processes=v * per_trial / 3600 / speed4,
floor_cpu_hours_10us=v * n_splines * 10e-6 / 3600) for k, v in need.items()]
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.1f}"))
print("notebook limit: 1 hour of wall time")
EEMD of one 150-bin series, 1,000 ensembles, 10 fixed sifts: 33.0 s on one core, 14.1 s on four processes; 7 IMFs; 16 EEMDs per trial; 528 CPU s per trial
one spline fit and evaluation: 242 us; splines per trial 1,920,000; floor per trial 464 s at this cost, 19 s at 10 us per spline
subset trials cpu_hours wall_hours_4_processes floor_cpu_hours_10us
figs4a, figs4b: every social Eyes trial 19474 2855.2 1217.2 103.9
slope claims: first 90 stimulations of the 76 epoch days 6840 1002.9 427.5 36.5
of which OFC days (25) 2250 329.9 140.6 12.0
of which dmPFC days (25) 2250 329.9 140.6 12.0
of which ACCg days (26) 2340 343.1 146.3 12.5
notebook limit: 1 hour of wall time
Result: blocked (too_large). At the Yang code's settings, which both readings share (1,000 ensembles, 10 fixed sifts), the smallest subset any of the six claims needs, the OFC epoch days of ofc-slope-unchanged-by-epoch, takes the CPU and wall hours printed above, and figs4a needs every social trial. Each is far beyond the one-hour limit of a notebook, and the spline floor shows that a compiled implementation at 10 us per spline would still need at least 12 CPU hours for the smallest subset; that fits in an hour only on a dozen or more cores, and on the four processes timed here it does not. The relative causal strength is therefore not computed and the six claims that use it are blocked.
fig5a-trigger-duration-unchanged¶
Microstimulations of OFC or dmPFC did not change the duration of the Eyes fixation that triggered the trial.
The M1 Eyes event in progress at the trial time (start at or before it, end at or after it). Readings: (a) M1-only and mutual Eyes events, (b) M1-only. The trials with no Eyes event in progress are broken down by area and trial type, and a variant counts an Eyes event that starts within 20 ms after the trial as the trigger.
trig = {}
for rd, tabs in [("a", ["m1_eye", "mutual_eye"]), ("b", ["m1_eye"])]:
ev = events[(events.who == "m1") & events.table.isin(tabs)]
p = ms.window_pairs(soc, ev, -10.0, 0.0005)
p = p[p.stop >= p.start - p.lag]
p["dur"] = (p.stop - p.start) * 1000
t = soc.copy()
t["trig"] = t.index.map(p.groupby("trial").dur.max()).astype(float)
trig[rd] = t
print(f"reading ({rd}): trials with an Eyes event in progress {int(t.trig.notna().sum())} of {len(t)}; "
f"median duration {t.trig.median():.0f} ms")
show(day_table(ms.day_means(t, "trig"), subjects=(None,)), "duration (ms), day level", degrees=False)
# the trials with no Eyes event in progress (reading a)
t = trig["a"].assign(none=lambda d: d.trig.isna())
print("\nfraction of trials with no Eyes event in progress, by area and trial type:")
print(t.groupby(["region", "type"]).none.agg(["mean", "size"]).unstack("type").round(3))
print("by monkey and trial type:")
print(t.groupby(["subject", "type"]).none.mean().unstack("type").round(3))
dn = t.groupby(["date", "region", "subject", "type"]).none.mean().unstack("type").reset_index()
print("stimulation minus sham fraction per day, signed-rank over days:")
show(day_table(dn, subjects=(None,)), "", degrees=False)
out = t[t.none]
ev_a = events[(events.who == "m1") & events.table.isin(["m1_eye", "mutual_eye"])]
nxt = ms.window_pairs(out, ev_a, 0.0005, 2.0).groupby("trial").lag.min().reindex(out.index)
print(f"of the {len(out)} trials with none: an Eyes event starts within 20 ms after the trial in "
f"{int((nxt <= 0.02).sum())}, within 2 s in {int(nxt.notna().sum())}")
anyev = events[(events.who == "m1") & events.table.isin(["m1_everywhere", "m1_neface", "mutual_face", "mutual_mouth"])]
pa = ms.window_pairs(out, anyev, -10.0, 0.0005)
pa = pa[pa.stop >= pa.start - pa.lag].drop_duplicates("trial")
print(" inside another M1 event at the trial time:", len(pa), "; by table:", pa.table.value_counts().to_dict())
# variant: an Eyes event that starts within 20 ms after the trial counts as the trigger
ev = ev_a
p2 = ms.window_pairs(soc, ev, -10.0, 0.0205)
p2 = p2[(p2.stop >= p2.start - p2.lag) | (p2.lag > 0)]
p2["dur"] = (p2.stop - p2.start) * 1000
tv = soc.copy()
tv["trig"] = tv.index.map(p2.sort_values("lag").groupby("trial").dur.first()).astype(float)
print(f"\nvariant, Eyes event in progress or starting within 20 ms (reading a): {int(tv.trig.notna().sum())} trials")
show(day_table(ms.day_means(tv, "trig"), subjects=(None,)), "duration (ms), day level", degrees=False)
reading (a): trials with an Eyes event in progress 13779 of 19474; median duration 367 ms duration (ms), day level region subject n median_diff mean_diff p OFC both 27 8.148 14.94 0.6109 dmPFC both 27 31.04 46.01 0.01866 ACCg both 27 3.798 25.05 0.3012
reading (b): trials with an Eyes event in progress 12647 of 19474; median duration 358 ms
duration (ms), day level
region subject n median_diff mean_diff p
OFC both 27 9.957 24.14 0.3608
dmPFC both 27 28.34 46.25 0.02802
ACCg both 27 39.24 36.58 0.1349
fraction of trials with no Eyes event in progress, by area and trial type:
mean size
type sham stim sham stim
region
ACCg 0.300 0.291 3282 3217
OFC 0.295 0.301 3294 3242
dmPFC 0.294 0.273 3240 3199
by monkey and trial type:
type sham stim
subject
L 0.362 0.354
T 0.212 0.205
stimulation minus sham fraction per day, signed-rank over days:
region subject n median_diff mean_diff p
OFC both 27 0.0045 0.00642 0.5781
dmPFC both 27 -0.01885 -0.02408 0.06923
ACCg both 27 -0.01325 -0.005868 0.5303
of the 5695 trials with none: an Eyes event starts within 20 ms after the trial in 609, within 2 s in 1778
inside another M1 event at the trial time: 4227 ; by table: {'m1_everywhere': 2036, 'm1_neface': 1600, 'mutual_face': 558, 'mutual_mouth': 33}
variant, Eyes event in progress or starting within 20 ms (reading a): 14388 trials
duration (ms), day level
region subject n median_diff mean_diff p
OFC both 27 15.04 10.49 0.5781
dmPFC both 27 29.99 45.19 0.02
ACCg both 27 3.858 25.41 0.3012
Result: partially reproduced. The duration of the Eyes event in progress at the trial does not differ for OFC (p = 0.61 with mutual events, 0.36 without), as the paper reports, but for dmPFC it is longer on stimulation trials (median day difference 31 ms, p = 0.019; 28 ms, p = 0.028 without mutual events), where the paper reports p = 0.269. The 5,695 trials with no Eyes event in progress (29%) are spread evenly over areas and trial types (27.3% to 30.1%; stimulation minus sham per day, dmPFC p = 0.069, OFC 0.58); 4,227 of them fall inside another M1 fixation, mostly outside the face or on the non-eye face, and 609 have an Eyes event starting within 20 ms after the trial. Counting those 609 as the trigger leaves the result unchanged (dmPFC p = 0.020, OFC 0.58).
fig5b-pupil-unchanged¶
Pupil size did not differ between microstimulation and sham trials in any 10 ms bin of the 700 ms after trial onset, for any area.
pup = np.vstack([r["pupil"] for r in runs])
offset = np.concatenate([[0], np.cumsum([len(r["trials"]) for r in runs])[:-1]])
sel = soc.assign(prow=offset[soc.run_id.to_numpy()] + soc.row.to_numpy())
sel = sel[sel.in_face_700]
print("social trials with gaze in the whole Face for 700 ms:", len(sel), "of", len(soc))
P = pd.DataFrame(pup[sel.prow.to_numpy()], index=sel.index)
P[["date", "region", "subject", "type"]] = sel[["date", "region", "subject", "type"]]
dmp = P.groupby(["date", "region", "subject", "type"]).mean()
rows = []
for a in AREAS:
x = dmp.xs(a, level="region")
st = x.xs("stim", level="type")
sh = x.xs("sham", level="type")
common = st.index.intersection(sh.index)
ps = np.array([signrank(st.loc[common, b], sh.loc[common, b])[1] for b in range(pup.shape[1])])
rows.append(dict(region=a, days=len(common), bins=len(ps), bins_p_below_05=int((ps < 0.05).sum()),
min_p=ps.min(), bins_stim_larger=int((st.loc[common].median() > sh.loc[common].median()).sum())))
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.4g}"))
print("pupil (a.u.) of the selected trials, range:", np.nanmin(P.iloc[:, :70].to_numpy()), np.nanmax(P.iloc[:, :70].to_numpy()))
social trials with gaze in the whole Face for 700 ms: 2198 of 19474
region days bins bins_p_below_05 min_p bins_stim_larger OFC 27 70 0 0.3126 29 dmPFC 27 70 0 0.5149 33 ACCg 27 70 0 0.1056 10 pupil (a.u.) of the selected trials, range: 1112.7 5312.1
Result: reproduced. Of the 70 bins of 10 ms, none gives p < 0.05 in any area (smallest p 0.31, 0.51 and 0.11 for OFC, dmPFC and ACCg), on the 2,198 trials whose gaze stayed in the whole Face for 700 ms.
fig5cd-saccade-counts-unchanged¶
The numbers of microsaccades and macrosaccades in the post-gaze epoch did not differ between microstimulation and sham trials for OFC or dmPFC.
The files hold no saccade events (data check). Saccades are detected here with a Python port of
the Otero-Millan et al. (2014) unsupervised clustering detector, from the authors' MATLAB code at
a pinned commit, called as the paper's own code calls it (script/run_saccade_detection.m and
script/detect_microsaccades.m of the paper's repository at cc0bf994): M1's gaze in degrees
from the first valid time sample, 1 kHz, NaN samples as blinks (microstim.detect_saccades).
This makes the claim with_inputs for that software, as the claims file says. A saccade is a
microsaccade when it lies strictly within an M1 fixation of the files (all M1 fixations, M1-only
and mutual), reading (a), the claims file's; reading (b) is the paper's code
(script/label_saccades.m): its start at least 50 samples after a fixation's start and at least
50 before its end. Counts per trial are of saccades starting in [0, 1.5) s after the trial; day
means per trial type; signed-rank over days.
sac = []
for k, r in enumerate(runs):
if r["saccades"] is not None and len(r["saccades"]):
sac.append(r["saccades"].assign(run_id=k))
sac = pd.concat(sac, ignore_index=True).merge(meta[["run_id", "date", "subject", "region", "roi", "cond"]], on="run_id")
sac["duration_ms"] = (sac.end - sac.start + 1)
sm = meta[meta.cond == "Social"]
print(f"social runs: {len(sm)}; saccades detected: {len(sac)} ({len(sac) / len(sm):.0f} per run); "
f"velocity peaks per run, median {sm.sacc_n_peaks.median():.0f}; candidate events after limits and validity, "
f"median {sm.sacc_n_candidates.median():.0f}; clusters chosen: {sm.sacc_n_clusters.value_counts().sort_index().to_dict()}")
print("runs where the MATLAB merge of overlapping saccades would not pair starts and ends (merged as a union here):",
int(sm.sacc_merged_as_union.sum()), "of", len(sm))
print("saccade duration (ms) quantiles:", sac.duration_ms.quantile([0.05, 0.25, 0.5, 0.75, 0.95]).to_dict())
print("amplitude (deg) quantiles:", sac.amplitude.quantile([0.05, 0.5, 0.95]).round(2).to_dict(),
"; peak velocity (deg/s) quantiles:", sac.peak_velocity.quantile([0.05, 0.5, 0.95]).round(0).to_dict())
# micro or macro: (a) strictly within an M1 fixation, (b) start at least 50 ms after a fixation's
# start and 50 ms before its end (label_saccades.m)
fixall = events[(events.who == "m1") & events.table.isin(["m1_all", "mutual_wholeface"]) & (events.cond == "Social")]
sp = ms.window_pairs(sac.assign(time=sac.start_time), fixall, -30.0, 0.0005)
sp["sac_start"] = sp.start - sp.lag
sp["sac_end"] = sac.loc[sp.trial, "end_time"].to_numpy()
in_a = sp[(sp.start < sp.sac_start) & (sp.sac_end < sp.stop)].trial.unique()
in_b = sp[(sp.sac_start - sp.start >= 0.050) & (sp.stop - sp.sac_start >= 0.050)].trial.unique()
sac["micro_a"] = sac.index.isin(in_a)
sac["micro_b"] = sac.index.isin(in_b)
print(f"microsaccades: reading (a) {int(sac.micro_a.sum())} ({sac.micro_a.mean():.3f}), "
f"reading (b) {int(sac.micro_b.sum())} ({sac.micro_b.mean():.3f}) of {len(sac)}")
ssoc = sac[sac.roi == "Eyes"]
wp = ms.window_pairs(soc, ssoc.assign(start=ssoc.start_time), 0.0, ms.POST_GAZE_S)
sc = soc.copy()
for rd in ["a", "b"]:
sc[f"micro_{rd}"] = sc.index.map(wp[wp[f"micro_{rd}"]].groupby("trial").size()).fillna(0)
sc[f"macro_{rd}"] = sc.index.map(wp[~wp[f"micro_{rd}"]].groupby("trial").size()).fillna(0)
print("saccades per post-gaze epoch, means (a):", sc[["micro_a", "macro_a"]].mean().round(3).to_dict(),
"(b):", sc[["micro_b", "macro_b"]].mean().round(3).to_dict())
rows = []
for rd in ["a", "b"]:
for what in ["micro", "macro"]:
for a in AREAS:
rows.append(dict(reading=rd, saccades=what, **day_test(ms.day_means(sc, f"{what}_{rd}"), a)))
sct = pd.DataFrame(rows)
print(sct.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
for rd in ["a", "b"]:
x = sct[(sct.reading == rd) & sct.region.isin(["OFC", "dmPFC"])]
print(f"reading ({rd}): OFC and dmPFC tests with p < 0.05: {int((x.p < 0.05).sum())} of {len(x)}")
print("by monkey, reading (a):")
rows = [dict(saccades=w, **day_test(ms.day_means(sc, f"{w}_a"), a, s)) for w in ["micro", "macro"] for a in ["OFC", "dmPFC"] for s in ["L", "T"]]
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.4g}"))
social runs: 1179; saccades detected: 454642 (386 per run); velocity peaks per run, median 1384; candidate events after limits and validity, median 557; clusters chosen: {2.0: 817, 3.0: 175, 4.0: 187}
runs where the MATLAB merge of overlapping saccades would not pair starts and ends (merged as a union here): 202 of 1179
saccade duration (ms) quantiles: {0.05: 41.0, 0.25: 77.0, 0.5: 147.0, 0.75: 383.0, 0.95: 1062.0}
amplitude (deg) quantiles: {0.05: 3.02, 0.5: 13.71, 0.95: 34.58} ; peak velocity (deg/s) quantiles: {0.05: 274.0, 0.5: 620.0, 0.95: 1150.0}
microsaccades: reading (a) 79753 (0.175), reading (b) 257222 (0.566) of 454642
saccades per post-gaze epoch, means (a): {'micro_a': 0.516, 'macro_a': 1.968} (b): {'micro_b': 1.275, 'macro_b': 1.208}
reading saccades region subject n median_diff mean_diff p
a micro OFC both 27 0.01739 0.02366 0.06923
a micro dmPFC both 27 -0.0001139 -0.00283 0.8889
a micro ACCg both 27 0.03018 0.08268 0.008889
a macro OFC both 27 -0.0694 -0.06774 0.04356
a macro dmPFC both 27 -0.03392 -0.04047 0.2386
a macro ACCg both 27 -0.05 -0.09441 0.003878
b micro OFC both 27 0.007194 -0.01335 0.7317
b micro dmPFC both 27 -0.06384 -0.03284 0.1482
b micro ACCg both 27 0.03465 0.03891 0.2792
b macro OFC both 27 0.007895 -0.03073 0.6938
b macro dmPFC both 27 -0.007204 -0.01046 0.9341
b macro ACCg both 27 -0.04695 -0.05064 0.01506
reading (a): OFC and dmPFC tests with p < 0.05: 1 of 4
reading (b): OFC and dmPFC tests with p < 0.05: 0 of 4
by monkey, reading (a):
saccades region subject n median_diff mean_diff p
micro OFC L 15 0.04143 0.03589 0.0946
micro OFC T 12 0.002692 0.008371 0.3804
micro dmPFC L 15 -0.02967 -0.009476 0.7197
micro dmPFC T 12 0.003876 0.005479 0.4131
macro OFC L 15 -0.1222 -0.09905 0.06372
macro OFC T 12 -0.03574 -0.0286 0.1514
macro dmPFC L 15 0.07597 -0.02574 0.7615
macro dmPFC T 12 -0.08487 -0.05889 0.1514
# the detector's velocity on one social run, during the files' fixations, against its 5 deg/s
# limit for the start and end of a saccade
import io
import h5py
rid0 = int(sm.run_id.iloc[0])
with h5py.File(io.BytesIO(requests.get(f"{ms.API}/assets/{meta.asset_id[rid0]}/download/", timeout=300).content), "r") as f:
g0 = f["acquisition/EyeTracking_m1/spatialseries/data"][:] * ms.DEG_PER_PX
t0_ = f["acquisition/t/data"][:]
v0 = ms._differentiate(g0, 1000.0)
vel0 = np.hypot(v0[:, 0], v0[:, 1])
ev0 = events[(events.run_id == rid0) & (events.table == "m1_all")]
ok_t = ~np.isnan(t0_)
idx = np.flatnonzero(ok_t)
a0 = idx[np.searchsorted(t0_[ok_t], ev0.start.to_numpy())]
b0 = idx[np.minimum(np.searchsorted(t0_[ok_t], ev0.stop.to_numpy()), ok_t.sum() - 1)]
infix = np.concatenate([vel0[a:b] for a, b in zip(a0, b0)])
lay0 = ms.file_layout(meta.asset_id[rid0])
print("groups of one social file,", meta.path[rid0])
for g, names in lay0.items():
print(f" {g}: {names}")
print("tables named like saccade in any group of that file:",
[n for names in lay0.values() if names for n in names if "sacc" in n.lower()])
print("run", meta.path[rid0], "; fixations", len(ev0))
print("velocity (deg/s) inside fixations, quantiles 5/25/50/75/95%:", np.nanpercentile(infix, [5, 25, 50, 75, 95]).round(1),
"; fraction under 5 deg/s:", round(float(np.nanmean(infix < 5)), 4))
print("velocity (deg/s) of all samples, median:", round(float(np.nanmedian(vel0)), 1),
"; fraction under 5 deg/s:", round(float(np.nanmean(vel0 < 5)), 4))
groups of one social file, sub-L/sub-L_ses-20190718_obj-jvhdtl_behavior.nwb /: ['acquisition', 'analysis', 'file_create_date', 'general', 'identifier', 'intervals', 'processing', 'scratch', 'session_description', 'session_start_time', 'specifications', 'stimulus', 'timestamps_reference_time'] intervals: ['m1_all', 'm1_everywhere', 'm1_eye', 'm1_mouth', 'm1_neface', 'm1_nenmface', 'm1_wholeface', 'm2_all', 'm2_everywhere', 'm2_eye', 'm2_mouth', 'm2_neface', 'm2_nenmface', 'm2_wholeface', 'mutual_face', 'mutual_wholeface'] processing: [] acquisition: ['EyeTracking_m1', 'EyeTracking_m2', 'pupil_m1', 'pupil_m2', 'roi_definitions', 'stimulation_events', 't'] scratch: ['stimulated_region', 'stimulation_condition', 'stimulation_roi'] analysis: [] tables named like saccade in any group of that file: [] run sub-L/sub-L_ses-20190718_obj-jvhdtl_behavior.nwb ; fixations 463 velocity (deg/s) inside fixations, quantiles 5/25/50/75/95%: [ 4.6 11.6 19.7 32.5 127.8] ; fraction under 5 deg/s: 0.0538 velocity (deg/s) of all samples, median: 20.3 ; fraction under 5 deg/s: 0.0526
# check of the port on synthetic gaze: fixations joined by saccades of 0.3 to 15 deg with a
# raised-cosine profile and main-sequence durations (2.2 ms per deg + 21 ms), 60 s at 1 kHz,
# with white noise of three levels; the last gives a velocity inside fixations like the data's
rng_s = np.random.default_rng(5959)
def synth_gaze(noise, T=60.0, rate=1000):
n = int(T * rate)
x, y, true = np.zeros(n), np.zeros(n), []
t, px, py = 1500, 0.0, 0.0
while t < n - 1500:
amp, ang = rng_s.uniform(0.3, 15), rng_s.uniform(0, 2 * np.pi)
dur = int(round(2.2 * amp + 21))
prof = (1 - np.cos(np.pi * np.linspace(0, 1, dur))) / 2
nx, ny = px + amp * np.cos(ang), py + amp * np.sin(ang)
if abs(nx) > 15 or abs(ny) > 15:
nx, ny = px - amp * np.cos(ang), py - amp * np.sin(ang)
x[t:t + dur], y[t:t + dur] = px + (nx - px) * prof, py + (ny - py) * prof
true.append((t, t + dur - 1))
px, py = nx, ny
t2 = t + dur + int(rng_s.uniform(250, 900))
x[t + dur:t2], y[t + dur:t2] = px, py
t = t2
x[t:], y[t:] = px, py
return x + rng_s.normal(0, noise, n), y + rng_s.normal(0, noise, n), np.array(true)
rows = []
for noise in [0.005, 0.02, 0.05]:
x, y, tr = synth_gaze(noise)
d = ms.detect_saccades(x, y)
v = ms._differentiate(np.column_stack([x, y]), 1000.0)
sp = np.hypot(v[:, 0], v[:, 1])
still = np.ones(len(x), bool)
for a, b in tr:
still[a - 5:b + 6] = False
still[:1000] = False
found = np.array([((d.start <= b) & (d.end >= a)).any() for a, b in tr])
n_true = np.array([((tr[:, 0] <= e) & (tr[:, 1] >= s)).sum() for s, e in zip(d.start, d.end)])
rows.append(dict(noise_deg=noise, fix_velocity_median=np.median(sp[still]), fix_under_5=np.mean(sp[still] < 5),
true=len(tr), detected=len(d), true_found=found.mean(), detected_matching_none=(n_true == 0).mean(),
detected_spanning_2_or_more=(n_true >= 2).mean(), true_duration_median_ms=np.median(tr[:, 1] - tr[:, 0] + 1),
detected_duration_median_ms=np.median(d.end - d.start + 1),
detected_duration_95_ms=np.percentile(d.end - d.start + 1, 95)))
print(pd.DataFrame(rows).to_string(index=False, float_format=lambda v: f"{v:.3g}"))
noise_deg fix_velocity_median fix_under_5 true detected true_found detected_matching_none detected_spanning_2_or_more true_duration_median_ms detected_duration_median_ms detected_duration_95_ms
0.005 1.96 0.989 93 93 1 0 0 38 38 49.4
0.02 7.84 0.248 95 87 0.916 0 0 39 47 61.7
0.05 19.8 0.043 90 78 0.867 0 0 37.5 72.5 172
Result: partially reproduced. With the claims file's reading (a), a microsaccade lying strictly within an M1 fixation, microsaccade counts do not change for OFC (p = 0.069) or dmPFC (p = 0.89), and macrosaccade counts do not change for dmPFC (p = 0.24), but OFC stimulation trials have fewer macrosaccades (median day difference -0.069 per epoch, p = 0.044; paper 0.055, also the closest to significance of its four tests). With the code's reading (b), none of the four tests gives p < 0.05 (smallest 0.15). The detector, run as the paper's code calls it, finds 386 saccades per run with a median duration of 147 ms (95th percentile 1,062 ms): on these data the gaze velocity inside fixations has a median of about 20 deg/s and is under the detector's 5 deg/s saccade limit for only about 5% of samples, so many detected saccades extend over several fixations. In 202 of the 1,179 social runs the MATLAB code's merge of overlapping saccades would not pair starts and ends; they are merged as a union here. The port was not run against the MATLAB code, which is not available here. On synthetic gaze with little noise (velocity inside fixations 2.0 deg/s median) it finds every one of 93 saccades with their durations (median 38 ms detected and true); at white noise that gives the data's velocity inside fixations (19.8 deg/s median, 4.3% under 5 deg/s) it finds 86.7% and their median duration doubles (72.5 against 37.5 ms; 95th percentile 172 ms), with no detection spanning two saccades. The data's durations are longer still (median 147 ms, 95th percentile 1,062 ms), so their noise is not only white. The outcome therefore rests on an assumption: that this detection, on the files' gaze unfiltered as the paper's code passes it, stands in for the paper's, whose detected saccades are not reported. Reading (b) gives a different count of failing tests, so the outcome also depends on the micro/macro labelling.
fig5e-h-saccade-kinematics-unchanged¶
Microstimulations of OFC or dmPFC did not change macrosaccade peak velocity, the slope of peak velocity on amplitude, or that slope by saccade direction.
The first macrosaccade starting in the post-gaze epoch of each trial (both readings of the
macrosaccade). Day mean peak velocity per trial type, signed-rank. Per day, a least-squares
line of peak velocity on amplitude for stimulation and for sham trials (at least 3 saccades
each), the slope difference tested against zero by signed-rank; the permutation test shuffles
the trial labels 1,000 times within each day, takes the median slope difference over days for
each shuffle, and gives the fraction of shuffled medians beyond the real median, on its side.
Direction groups II, IC, CI and CC by the hemifields of the start and end positions (against the
centre of roi_m2_face), slope difference against zero per group. Seed 5 (a standard choice).
def slope(x, y):
x, y = np.asarray(x, float), np.asarray(y, float)
if len(x) < 3 or np.ptp(x) == 0:
return np.nan
return np.polyfit(x, y, 1)[0]
def hemi_of(x, m2x, subject):
right = x > m2x
return np.where((subject == "T") == right, "I", "C")
rng_k = np.random.default_rng(5)
kin_rows, perm_rows, dir_rows = [], [], []
for rd in ["a", "b"]:
mac = wp[~wp[f"micro_{rd}"]].sort_values(["trial", "lag"]).drop_duplicates("trial")
fm = soc.loc[mac.trial, ["date", "region", "subject", "type", "m2face_x"]].assign(
pv=mac.peak_velocity.to_numpy(), amp=mac.amplitude.to_numpy(),
d=[a + b for a, b in zip(hemi_of(mac.x0.to_numpy(), soc.loc[mac.trial, "m2face_x"].to_numpy(), soc.loc[mac.trial, "subject"].to_numpy()),
hemi_of(mac.x1.to_numpy(), soc.loc[mac.trial, "m2face_x"].to_numpy(), soc.loc[mac.trial, "subject"].to_numpy()))])
print(f"reading ({rd}): trials with a first macrosaccade {len(fm)} of {len(soc)}; direction groups {fm.d.value_counts().to_dict()}")
for a in AREAS:
kin_rows.append(dict(reading=rd, test="peak velocity", **day_test(ms.day_means(fm, "pv"), a)))
sd = []
for (date, region, subject), g in fm.groupby(["date", "region", "subject"]):
s1, s0 = slope(g[g.type == "stim"].amp, g[g.type == "stim"].pv), slope(g[g.type == "sham"].amp, g[g.type == "sham"].pv)
null = []
lab = g.type.to_numpy()
for _ in range(1000):
pl = rng_k.permutation(lab)
null.append(slope(g.amp[pl == "stim"], g.pv[pl == "stim"]) - slope(g.amp[pl == "sham"], g.pv[pl == "sham"]))
sd.append(dict(date=date, region=region, subject=subject, diff=s1 - s0, null=np.array(null)))
for dg, gg in g.groupby("d"):
dir_rows.append(dict(reading=rd, date=date, region=region, group=dg,
diff=slope(gg[gg.type == "stim"].amp, gg[gg.type == "stim"].pv) - slope(gg[gg.type == "sham"].amp, gg[gg.type == "sham"].pv)))
sd = pd.DataFrame(sd)
for a in AREAS:
x = sd[(sd.region == a) & sd["diff"].notna()]
n, p = signrank(x["diff"])
kin_rows.append(dict(reading=rd, test="slope difference against zero", region=a, subject="both", n=n,
median_diff=np.median(x["diff"]), mean_diff=np.mean(x["diff"]), p=p))
real = np.median(x["diff"])
nullm = np.nanmedian(np.vstack(x.null.to_numpy()), axis=0)
pp = np.mean(nullm >= real) if real >= 0 else np.mean(nullm <= real)
kin_rows.append(dict(reading=rd, test="slope difference, permutation", region=a, subject="both", n=len(x),
median_diff=real, mean_diff=np.nan, p=pp))
kin = pd.DataFrame(kin_rows)
print(kin.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
dr = pd.DataFrame(dir_rows).dropna(subset=["diff"])
dtab = pd.DataFrame([dict(reading=rd, region=a, group=gname, n=signrank(g["diff"])[0], median_diff=np.median(g["diff"]), p=signrank(g["diff"])[1])
for (rd, a, gname), g in dr.groupby(["reading", "region", "group"])])
print("\nslope difference by direction group, against zero:")
print(dtab.to_string(index=False, float_format=lambda v: f"{v:.4g}"))
for rd in ["a", "b"]:
x = pd.concat([kin[(kin.reading == rd) & kin.region.isin(["OFC", "dmPFC"])].p,
dtab[(dtab.reading == rd) & dtab.region.isin(["OFC", "dmPFC"])].p])
print(f"reading ({rd}): OFC and dmPFC tests with p < 0.05: {int((x < 0.05).sum())} of {len(x)}")
reading (a): trials with a first macrosaccade 15961 of 19474; direction groups {'II': 6555, 'CC': 3926, 'IC': 2994, 'CI': 2486}
reading (b): trials with a first macrosaccade 11354 of 19474; direction groups {'II': 5651, 'IC': 2481, 'CC': 2082, 'CI': 1140}
reading test region subject n median_diff mean_diff p
a peak velocity OFC both 27 -19.23 -11.73 0.05208
a peak velocity dmPFC both 27 -5.567 0.1022 0.8779
a peak velocity ACCg both 27 15.67 16.77 0.07723
a slope difference against zero OFC both 27 -0.5536 -0.04129 0.6964
a slope difference, permutation OFC both 27 -0.5536 NaN 0.284
a slope difference against zero dmPFC both 27 1.142 0.5306 0.5944
a slope difference, permutation dmPFC both 27 1.142 NaN 0.186
a slope difference against zero ACCg both 27 1.3 2.973 0.09539
a slope difference, permutation ACCg both 27 1.3 NaN 0.163
b peak velocity OFC both 27 -7.27 8.86 0.9717
b peak velocity dmPFC both 27 2.727 19.75 0.9906
b peak velocity ACCg both 27 -0.2571 15.15 0.6617
b slope difference against zero OFC both 27 1.242 2.212 0.1775
b slope difference, permutation OFC both 27 1.242 NaN 0.131
b slope difference against zero dmPFC both 27 1.971 4.571 0.0299
b slope difference, permutation dmPFC both 27 1.971 NaN 0.064
b slope difference against zero ACCg both 27 -0.2473 -3.598 0.8408
b slope difference, permutation ACCg both 27 -0.2473 NaN 0.442
slope difference by direction group, against zero:
reading region group n median_diff p
a ACCg CC 27 2.637 0.03621
a ACCg CI 24 4.751 0.04568
a ACCg IC 25 -0.9082 0.5965
a ACCg II 27 5.761 0.0552
a OFC CC 26 -0.9642 0.7454
a OFC CI 25 -5.251 0.07098
a OFC IC 26 0.5785 0.3403
a OFC II 27 0.3992 0.9341
a dmPFC CC 24 -1.043 0.7898
a dmPFC CI 26 -1.775 0.5153
a dmPFC IC 27 0.6524 0.5781
a dmPFC II 27 2.02 0.2792
b ACCg CC 26 0.3166 0.7835
b ACCg CI 23 0.8953 0.5399
b ACCg IC 23 1.812 0.8229
b ACCg II 26 -0.6742 0.7078
b OFC CC 26 3.907 0.6348
b OFC CI 23 -1.576 0.3604
b OFC IC 20 2.073 0.04844
b OFC II 24 2.117 0.2643
b dmPFC CC 25 -5.467 0.03668
b dmPFC CI 23 1.574 0.9881
b dmPFC IC 25 3.497 0.07098
b dmPFC II 22 2.108 0.6102
reading (a): OFC and dmPFC tests with p < 0.05: 0 of 14
reading (b): OFC and dmPFC tests with p < 0.05: 3 of 14
Result: reproduced. With reading (a) of the macrosaccade, none of the 14 OFC and dmPFC tests gives p < 0.05: peak velocity OFC p = 0.052 and dmPFC 0.88 (paper 0.848 and 1), slope difference against zero 0.70 and 0.59 (paper 0.665, 0.904), permutation 0.28 and 0.19 (paper 0.515, 0.164), and by direction group the smallest is 0.071 (OFC CI). With the code's reading (b), 3 of the 14 do (dmPFC slope difference 0.030, OFC IC 0.048, dmPFC CC 0.037), so the result depends on how microsaccades are told from macrosaccades. It also rests on the assumption stated for fig5cd-saccade-counts-unchanged, that this detection on the files' unfiltered gaze stands in for the paper's; ACCg, outside the claim, gives p below 0.05 for two direction groups under (a).
Summary¶
| Claim | Outcome |
|---|---|
| design-days-per-area | reproduced |
| design-sessions-per-day | partially reproduced |
| fig1e-stim-counts-comparable | reproduced |
| figs1d-stim-counts-nonsocial | reproduced |
| figs1e-eyes-more-than-face | reproduced |
| figs1e-eyes-vs-rdm-comparable | not reproduced |
| fig2a-density-map-ofc | reproduced |
| fig2c-ofc-social-gaze-distance | reproduced |
| fig2c-dmpfc-accg-no-distance-effect | reproduced |
| fig2e-ofc-contralateral | reproduced |
| fig2e-dmpfc-accg-no-hemifield-effect | reproduced |
| fig2f-nonsocial-no-distance-effect | reproduced |
| fig2gh-nonsocial-no-hemifield-effect | reproduced |
| methods-ofc-distance-longer-windows | reproduced |
| methods-ofc-no-change-in-fixation-number | reproduced |
| methods-ofc-longer-face-looking | partially reproduced |
| methods-mouth-stimulation-no-effect | reproduced |
| methods-eyes-vs-mouth-bootstrap | not reproduced |
| fig3b-ofc-inter-looking-interval | reproduced |
| ili-nonsocial-no-effect | reproduced |
| ili-missing-trials-41 | reproduced |
| figs3b-ili-trial-level | reproduced |
| figs3c-reciprocation-day-level-none | reproduced |
| fig3d-ofc-reciprocation-latency | reproduced |
| partner-monkey-unchanged | partially reproduced |
| epoch-days-76-of-81 | reproduced |
| figs4a, figs4b, ofc-slope, fig4c, fig4d, dmpfc-slope-positive | blocked, too_large |
| fig5a-trigger-duration-unchanged | partially reproduced |
| fig5b-pupil-unchanged | reproduced |
| fig5cd-saccade-counts-unchanged (with outside software) | partially reproduced |
| fig5e-h-saccade-kinematics-unchanged (with outside software) | reproduced |
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. It starts from the fixation events in the
files, not from the paper's own pipeline. Trial selection, fixation sets, hemifield split
and the measures follow the claims file, which takes them from the paper and its code; the
details the claims file leaves open are taken as stated in microstim.py. Statistics use
scipy's signed-rank and rank-sum tests, which can differ from MATLAB's in the last digits for
ties and for n above 15. The causal decomposition was not run (its cost is timed in the notebook). Saccades were detected
with a Python port of the Otero-Millan et al. (2014) detector, called as the paper's code calls
it; the port uses scipy's k-means, and in 202 social runs merges overlapping saccades as a union
where the MATLAB code would stop. Both saccade outcomes rest on an assumption: that this detection,
on the files' unfiltered gaze, stands in for the paper's. The port was not run against the MATLAB
code (MATLAB is not available here); on synthetic gaze it recovers known saccades and their
durations, and at the data's noise level it lengthens them, as on the data.
What was found about the dataset. The dataset holds every day of the paper (81 Eyes and
15 Mouth days), with one file per run and no run deposited twice. The files carry the
region, triggering ROI and condition; the event tables are consistent with each other; the
trial times fall inside an M1 event in the triggering ROI on 71% of social Eyes trials against
5 to 20% for times shifted by 0.5 to 2 s, consistent with a shared clock. The files' own
laterality labels match the claims file's hemifield rule for 99.98% of M1 events. Monkey L
has 6.0 control runs per day where the paper reports 5. 6.3% of the M1 Eyes events of social runs and 10.7% of
those of Dot runs have a mean position outside the Eyes rectangle of their file, and the files do not say how events were
assigned to ROIs. Neither the files nor the paper give the padding of the Eyes ROI. The files hold no saccade
events. mutual_face is the face outside the Eyes and the Mouth while m1_neface includes the
Mouth; the table descriptions do not say so.
print(f"notebook run time: {(time.time() - t_start) / 60:.1f} min")
notebook run time: 20.9 min