Reproducing Chowdhury, Glaser & Miller (2020) from dandiset 000127¶
Paper: Chowdhury RH, Glaser JI, Miller LE (2020). Area 2 of primary somatosensory cortex encodes kinematics of the whole arm. eLife 9:e48198. doi:10.7554/eLife.48198 (open access; the full text and methods were used).
Data: DANDI:000127, version 0.220113.0359, streamed from the archive.
What can and cannot be tested¶
The paper has two experiments. Eight sessions from three monkeys compare reaching in two workspaces (Figures 4 to 7). Four sessions from two monkeys compare active reaches with passive bumps of the hand (Figures 8 to 11).
The dandiset holds one session: monkey H ("Han") in the active/passive task on 2017-12-07, as packaged for the Neural Latents Benchmark. So:
| Claim | Paper | Here |
|---|---|---|
fig4-whole-arm-fits-better to fig7-pd-shift |
Two-workspace experiment | Blocked: those sessions are not in the dandiset |
fig8-hand-speed-similar |
Hand speed is similar in active and passive trials in the first 120 ms | Tested |
fig8-arm-kinematics-separable |
Whole-arm kinematics tell active from passive trials (89%) | Tested with a stand-in |
fig8-neurons-separable |
Many neurons tell active from passive trials above chance | Tested |
fig9-whole-arm-fits-better |
The whole-arm model tends to beat the hand-only model, with many ties | Tested with a stand-in |
fig10-inconsistent-across-conditions |
Many neurons are fit poorly by one model across both conditions | Tested with a stand-in |
fig11-separability-vs-consistency |
The better a neuron separates the conditions, the worse one model fits it | Tested with a stand-in |
The stand-in. The paper's whole-arm model uses the positions and velocities of markers on the hand and elbow. The marker data is not in the dandiset. It does hold the seven joint angles and their velocities computed from those markers, so the whole-arm model here is hand kinematics plus joint kinematics. The paper's appendix reports that, on the two-workspace data, a joint-kinematics model gave predictions similar to its whole-arm model, though somewhat worse.
Everything below is from one session of one monkey, where the paper has four sessions from two.
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
import area2
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning, module="joblib")
plt.rcParams.update({"figure.dpi": 100, "axes.spines.top": False, "axes.spines.right": False})
ACTIVE, PASSIVE = "k", "#c0392b"
rng = np.random.default_rng(127)
Load the data and check it against its description¶
area2.load_session streams the file and reduces each successful trial to the 120 ms after
movement onset: the average hand and joint kinematics, and each neuron's spike count.
session = area2.load_session()
trials, passive = session.trials, session.trials.passive.to_numpy()
for name, value in session.sanity.items():
print(f"{name}: {value}")
sampling rate (Hz): 1000.0
recording length (s): 2223.4
trials in file: 826
successful trials with a movement onset: 364
trials analysed: 364
neurons: 65
hand position range (m): [-0.149, 0.14]
hand speed, 99th percentile (m/s): 0.463
joint angle range (deg): [-72, 106]
joint columns in file: columns=['shoulder_adduction', 'shoulder_rotation', 'shoulder_flexion', 'elbow_flexion', 'radial_pronation', 'wrist_flexion', 'wrist_abduction']
movement directions (deg): [0.0, 45.0, 90.0, 135.0, 180.0, 225.0, 270.0, 315.0]
onset minus bump time (s): min, median, max: [-0.049, -0.003, 0.027]
trials by result and bump: {'A/no bump': 187, 'A/bump': 203, 'I/no bump': 41, 'I/bump': 31, 'R/no bump': 193, 'R/bump': 171}
series share a timebase: True
by_direction = trials.groupby(["passive", "direction"]).size().unstack(fill_value=0)
by_direction.index = ["active", "passive"]
print(f"{(~passive).sum()} active and {passive.sum()} passive trials")
by_direction
193 active and 171 passive trials
| direction | 0.0 | 45.0 | 90.0 | 135.0 | 180.0 | 225.0 | 270.0 | 315.0 |
|---|---|---|---|---|---|---|---|---|
| active | 33 | 28 | 22 | 27 | 19 | 26 | 20 | 18 |
| passive | 27 | 14 | 18 | 21 | 19 | 22 | 25 | 25 |
Checks against the file's own descriptions and the paper:
- Units. The descriptions of hand position and velocity say centimetres, while the stated unit is metres with a conversion of 0.01. Applying the conversion gives a workspace about 0.3 m across and speeds under 1 m/s, which is plausible for a monkey reaching on a manipulandum, so the values are stored in centimetres and are converted to metres here.
- Directions. This session has eight target and bump directions. The paper's methods describe four.
- Onset. For passive trials the recorded movement onset is within a few milliseconds of the bump time for the typical trial (range -49 to +27 ms), so it marks the bump, as the analysis needs.
- Trial types. In the file's result column, "R" marks rewarded trials; only those are used, as in the paper.
Claim fig8-hand-speed-similar¶
"Average movement speed was similar between active and passive trials in this 120 ms window."
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4))
t = area2.SPEED_TIMES * 1000
for mask, color, label in ((~passive, ACTIVE, "active"), (passive, PASSIVE, "passive")):
mean, sd = session.speed[mask].mean(axis=0), session.speed[mask].std(axis=0)
axes[0].plot(t, mean, color=color, label=label)
axes[0].fill_between(t, mean - sd, mean + sd, color=color, alpha=0.15, lw=0)
axes[0].axvspan(0, area2.WINDOW * 1000, color="0.9", zorder=0)
axes[0].set_xlabel("Time from movement onset (ms)"); axes[0].set_ylabel("Hand speed (m/s)")
axes[0].legend(frameon=False)
bins = np.linspace(0, 0.45, 31)
axes[1].hist(trials.mean_speed[~passive], bins=bins, color=ACTIVE, alpha=0.5, label="active")
axes[1].hist(trials.mean_speed[passive], bins=bins, color=PASSIVE, alpha=0.5, label="passive")
axes[1].set_xlabel("Mean hand speed in the window (m/s)"); axes[1].set_ylabel("Trials")
plt.tight_layout()
speed = trials.groupby("passive").mean_speed.agg(["mean", "std", "size"]).round(3)
speed.index = ["active", "passive"]
print(speed.to_string())
a, p = trials.mean_speed[~passive], trials.mean_speed[passive]
print(f"\nPassive is {100 * (p.mean() / a.mean() - 1):.0f}% faster on average "
f"(Welch t-test p = {stats.ttest_ind(a, p, equal_var=False).pvalue:.1e}); "
f"difference in units of the pooled SD: {(p.mean() - a.mean()) / np.sqrt((a.var() + p.var()) / 2):.2f}")
hand_sep = area2.separability(trials[area2.HAND].to_numpy(), passive, rng)
chance = max(passive.mean(), 1 - passive.mean())
print(f"Trial type predicted from hand position and velocity alone: {hand_sep.mean():.0%} "
f"(chance {chance:.0%})")
speed_sep = area2.separability(trials.mean_speed.to_numpy(), passive, rng)
print(f"Trial type predicted from mean hand speed alone: {speed_sep.mean():.0%}")
mean std size active 0.183 0.039 193 passive 0.207 0.055 171 Passive is 13% faster on average (Welch t-test p = 3.5e-06); difference in units of the pooled SD: 0.50
Trial type predicted from hand position and velocity alone: 58% (chance 53%) Trial type predicted from mean hand speed alone: 58%
Result: reproduced. In the 120 ms window the two speed traces overlap, and they separate only afterwards, as in the paper's Figure 8B. The speeds are not identical: passive movements are 13% faster on average, half a standard deviation, and that difference is statistically clear. "Similar" in the paper is a qualitative statement with no threshold, and it holds in that sense.
Claim fig8-arm-kinematics-separable¶
"A linear discriminant analysis (LDA) classifier could predict the movement type 89% of the time, using only the whole-arm kinematics in the analysis window." (An average over the paper's four sessions.)
joint_sep = area2.separability(trials[area2.JOINT_KINEMATICS].to_numpy(), passive, rng)
arm_sep = area2.separability(trials[area2.HAND + area2.JOINT_KINEMATICS].to_numpy(), passive, rng)
pd.DataFrame({
"covariates": ["hand position and velocity", "joint angles and velocities", "hand and joints"],
"accuracy": [hand_sep.mean(), joint_sep.mean(), arm_sep.mean()],
"sd across folds": [hand_sep.std(), joint_sep.std(), arm_sep.std()],
}).round(3).assign(chance=round(chance, 3))
| covariates | accuracy | sd across folds | chance | |
|---|---|---|---|---|
| 0 | hand position and velocity | 0.584 | 0.054 | 0.53 |
| 1 | joint angles and velocities | 0.897 | 0.030 | 0.53 |
| 2 | hand and joints | 0.907 | 0.032 | 0.53 |
Result: reproduced. Joint kinematics tell active from passive trials 90% of the time in this session, against the paper's 89% averaged over four sessions, where hand kinematics alone manage 58%. The features are joint angles and velocities, not the paper's hand and elbow markers.
Claim fig8-neurons-separable¶
"Many neurons had an above chance separability index." The separability index of a neuron is the cross-validated accuracy of a linear discriminant predicting the trial type from that neuron's spike count in the window.
Neurons firing below 1 Hz in the window are left out of this and the model analyses. With a 120 ms window they have almost no spikes, and model fits for them are meaningless. The paper does not state a firing-rate criterion.
rate = session.counts.mean(axis=0) / area2.WINDOW
keep = rate >= area2.MIN_RATE
counts = session.counts[:, keep].astype(float)
print(f"{keep.sum()} of {len(keep)} neurons fire at {area2.MIN_RATE:.0f} Hz or more in the window")
sep_folds = np.array([area2.separability(counts[:, i], passive, rng) for i in range(counts.shape[1])])
sep = sep_folds.mean(axis=1)
sep_p = np.array([area2.corrected_t_test(f - chance) for f in sep_folds])
above = (sep_p < 0.05) & (sep > chance)
order = np.argsort(sep)
fig, ax = plt.subplots(figsize=(9, 3))
ax.bar(range(len(sep)), sep[order], color=np.where(above[order], "#3b6fb6", "0.75"))
ax.axhline(chance, color="k", ls="--", lw=1)
ax.set_ylim(0.4, 0.85); ax.set_xlabel("Neuron (sorted)"); ax.set_ylabel("Separability index")
plt.tight_layout()
print(f"{above.sum()} of {len(sep)} neurons are above chance ({chance:.0%}) at p < 0.05; "
f"median {np.median(sep):.2f}, highest {sep.max():.2f}")
50 of 65 neurons fire at 1 Hz or more in the window
20 of 50 neurons are above chance (53%) at p < 0.05; median 0.56, highest 0.80
Result: reproduced. 20 of 50 neurons separate the trial types above chance. The test is not corrected for 50 comparisons, but it is conservative, and 20 is far more than false positives would give.
Models of neural activity¶
As in the paper: a Poisson generalised linear model per neuron predicts its spike count in the window from kinematics, scored by pseudo-R2 on held-out trials over 20 repeats of 5-fold cross-validation, with models compared by the corrected resampled t-test.
- Hand-only: hand position and velocity in the plane (5 parameters; the paper's has 7, in three dimensions).
- Whole-arm stand-in: hand kinematics plus seven joint angles and velocities (19 parameters; the paper's whole-arm model has 13).
Two deviations from the paper's method, both to keep the fits stable at these low counts: a small ridge penalty on the weights, and pseudo-R2 for one trial type is computed on the pooled held-out predictions of a whole repeat instead of on single folds.
models = {"hand-only": trials[area2.HAND].to_numpy(),
"whole-arm": trials[area2.HAND + area2.JOINT_KINEMATICS].to_numpy()}
cv = area2.cross_validate(models, counts, passive, within="whole-arm", rng=rng)
pr2 = {name: scores.mean(axis=0) for name, scores in cv["fold"].items()}
pd.DataFrame(pr2).describe().loc[["min", "25%", "50%", "75%", "max"]].round(3)
| hand-only | whole-arm | |
|---|---|---|
| min | -0.060 | -0.069 |
| 25% | -0.019 | 0.023 |
| 50% | 0.028 | 0.089 |
| 75% | 0.104 | 0.260 |
| max | 0.617 | 0.620 |
Claim fig9-whole-arm-fits-better¶
"The whole-arm model again tended to out-perform the hand-only model. However, there were many more neurons for which the difference between models was insignificant compared to the two-workspace experiment."
difference = cv["fold"]["whole-arm"] - cv["fold"]["hand-only"]
p_models = np.array([area2.corrected_t_test(difference[:, i]) for i in range(counts.shape[1])])
better = (p_models < 0.05) & (difference.mean(axis=0) > 0)
worse = (p_models < 0.05) & (difference.mean(axis=0) < 0)
fig, ax = plt.subplots(figsize=(4.2, 4.2))
ax.scatter(pr2["hand-only"][~(better | worse)], pr2["whole-arm"][~(better | worse)],
facecolors="none", edgecolors="0.4", s=24, label="no significant difference")
ax.scatter(pr2["hand-only"][better | worse], pr2["whole-arm"][better | worse],
color="#3b6fb6", s=24, label="significant difference")
lim = [-0.15, 0.8]
ax.plot(lim, lim, "k--", lw=1); ax.set_xlim(lim); ax.set_ylim(lim)
ax.set_xlabel("Hand-only pseudo-R2"); ax.set_ylabel("Whole-arm pseudo-R2")
ax.legend(frameon=False, fontsize=8, loc="lower right")
plt.tight_layout()
print(f"Whole-arm significantly better: {better.sum()} of {len(better)} neurons; "
f"hand-only significantly better: {worse.sum()}; no significant difference: {(~(better | worse)).sum()}")
print(f"Median pseudo-R2: hand-only {np.median(pr2['hand-only']):.3f}, whole-arm {np.median(pr2['whole-arm']):.3f}")
Whole-arm significantly better: 15 of 50 neurons; hand-only significantly better: 0; no significant difference: 35 Median pseudo-R2: hand-only 0.028, whole-arm 0.089
Result: reproduced. The whole-arm stand-in is significantly better for 15 of 50 neurons and the hand-only model for none; for 35 there is no significant difference. That is 70% without a difference, against 17% (50 of 288) in the paper's two-workspace experiment.
The comparison with the two-workspace experiment rests on the paper's number, since those sessions are not in the dandiset. The paper applies a Bonferroni correction to its model comparisons, which is not applied here; it would lower the 15. The counts here and below move by one or two with a different random split of the trials.
Claim fig10-inconsistent-across-conditions¶
"Many neurons in the active/passive task did not have this consistent relationship, indicated by the many neurons with negative pR2 values for the model trained over both conditions." The comparison is between the whole-arm model trained on both trial types and the same model trained on one, each scored on held-out trials of that type.
both = {t: cv["both"][("whole-arm", t)].mean(axis=0) for t in ("active", "passive")}
within = {t: cv["within"][t].mean(axis=0) for t in ("active", "passive")}
fig, axes = plt.subplots(1, 2, figsize=(8.4, 4.2))
table = []
for ax, t, color in zip(axes, ("active", "passive"), (ACTIVE, PASSIVE)):
ax.scatter(within[t], both[t], s=22, color=color, alpha=0.7)
lim = [-0.5, 0.8]
ax.plot(lim, lim, "k--", lw=1); ax.axhline(0, color="0.8", lw=1); ax.axvline(0, color="0.8", lw=1)
ax.set_xlim(lim); ax.set_ylim(lim); ax.set_title(f"Scored on {t} trials", fontsize=10)
ax.set_xlabel(f"Pseudo-R2, trained on {t} trials only"); ax.set_ylabel("Pseudo-R2, trained on both")
below = int((both[t] < lim[0]).sum())
if below:
ax.text(0.02, 0.02, f"{below} below axis", transform=ax.transAxes, fontsize=8)
table.append(dict(scored_on=t, neurons=len(both[t]),
negative_trained_on_both=int((both[t] < 0).sum()),
negative_trained_within=int((within[t] < 0).sum()),
positive_within_but_negative_both=int(((within[t] > 0) & (both[t] < 0)).sum()),
median_both=round(float(np.median(both[t])), 3),
median_within=round(float(np.median(within[t])), 3),
wilcoxon_p=stats.wilcoxon(within[t], both[t]).pvalue))
plt.tight_layout()
for t in ("active", "passive"):
print(f"{t}: pseudo-R2 range, trained on both {both[t].min():.2f} to {both[t].max():.2f}; "
f"trained within {within[t].min():.2f} to {within[t].max():.2f}")
log = area2.FIT_LOG
print(f"Model fits: {log['fits']}, not converged: {log['not_converged']}, "
f"with no spikes in the training trials: {log['no_training_spikes']}")
pd.DataFrame(table).set_index("scored_on")
active: pseudo-R2 range, trained on both -0.89 to 0.71; trained within -0.19 to 0.78 passive: pseudo-R2 range, trained on both -0.20 to 0.71; trained within -0.18 to 0.71 Model fits: 20000, not converged: 0, with no spikes in the training trials: 0
| neurons | negative_trained_on_both | negative_trained_within | positive_within_but_negative_both | median_both | median_within | wilcoxon_p | |
|---|---|---|---|---|---|---|---|
| scored_on | |||||||
| active | 50 | 19 | 17 | 5 | 0.026 | 0.084 | 0.017575 |
| passive | 50 | 23 | 25 | 3 | 0.025 | -0.003 | 0.323028 |
Result: partially reproduced. Read literally the claim holds: 19 of 50 neurons scored on active trials and 23 of 50 scored on passive trials have negative pseudo-R2 for the model trained on both conditions. But about as many are negative for models trained within one condition (17 and 25), so most of these are neurons that no kinematic model fits, not neurons with an inconsistent relationship. Only 5 and 3 neurons are positive when trained within a condition and negative when trained on both, which is the pattern the paper's argument needs.
There is a modest cost of training on both conditions for active trials, and none for passive trials.
The pseudo-R2 for one trial type is measured against that type's own mean firing in the held-out trials. Measured against a mean taken over both types, fewer negative values would be expected; that was not computed here. The paper does not say which it used.
Claim fig11-separability-vs-consistency¶
"Model consistency actually correlated negatively with the separability index." Consistency is the pseudo-R2 of the whole-arm model trained on both trial types.
def bootstrap_ci(x, y, n=2000):
idx = rng.integers(0, len(x), size=(n, len(x)))
r = [stats.spearmanr(x[i], y[i])[0] for i in idx]
return np.nanpercentile(r, [2.5, 97.5])
fig, axes = plt.subplots(1, 2, figsize=(8.4, 3.6), sharey=True)
table = []
for ax, t, color in zip(axes, ("active", "passive"), (ACTIVE, PASSIVE)):
ax.scatter(sep, both[t], s=22, color=color, alpha=0.7)
ax.axhline(0, color="0.8", lw=1); ax.axvline(chance, color="k", ls="--", lw=1)
ax.set_ylim(-1.0, 0.8); ax.set_xlabel("Separability index"); ax.set_title(f"Scored on {t} trials", fontsize=10)
rho, p = stats.spearmanr(sep, both[t])
lo, hi = bootstrap_ci(sep, both[t])
table.append(dict(scored_on=t, spearman=round(rho, 2), ci_low=round(lo, 2), ci_high=round(hi, 2), p=p,
pearson=round(stats.pearsonr(sep, both[t])[0], 2)))
axes[0].set_ylabel("Pseudo-R2, trained on both")
plt.tight_layout()
overall = stats.spearmanr(sep, pr2["whole-arm"])
print(f"Scored on all trials together: Spearman {overall[0]:.2f} (p = {overall[1]:.2f})")
pd.DataFrame(table).set_index("scored_on")
Scored on all trials together: Spearman 0.14 (p = 0.32)
| spearman | ci_low | ci_high | p | pearson | |
|---|---|---|---|---|---|
| scored_on | |||||
| active | -0.38 | -0.60 | -0.11 | 0.006108 | -0.54 |
| passive | -0.13 | -0.42 | 0.16 | 0.365415 | -0.06 |
Result: partially reproduced. Scored on active trials, the correlation is negative, as the paper reports: neurons that separate the trial types better are fit worse by a single model. Scored on passive trials, or on all trials together, it cannot be told apart from zero. The paper describes a negative correlation without restricting it to one trial type.
The Spearman correlation is the one to read. The Pearson value for active trials leans on a single neuron.
Summary¶
| Claim | Outcome | Note |
|---|---|---|
fig4-whole-arm-fits-better, fig5-consistent-across-workspaces, fig6-tuning-curves, fig7-pd-shift |
Blocked: data missing | The two-workspace sessions are not in the dandiset. |
fig8-hand-speed-similar |
Reproduced | Traces overlap in the window; passive is 13% faster on average. |
fig8-arm-kinematics-separable |
Reproduced | 90% here, 89% in the paper. |
fig8-neurons-separable |
Reproduced | 20 of 50 neurons above chance. |
fig9-whole-arm-fits-better |
Reproduced | Better for 15 of 50 neurons, worse for none. |
fig10-inconsistent-across-conditions |
Partially reproduced | Many negative fits, but as many when trained within a condition. |
fig11-separability-vs-consistency |
Partially reproduced | Negative in active trials only. |
How this differs from the paper's analysis
- One session from one monkey, where the paper has four from two.
- The whole-arm model is hand plus joint kinematics, because the marker positions the paper's model uses are not in the dandiset.
- The hand-only model is in two dimensions, the session has eight movement directions, neurons below 1 Hz are left out, the fits carry a small ridge penalty, and no Bonferroni correction is applied.
- Pseudo-R2 for one trial type is computed on a repeat's pooled held-out predictions, against that type's own mean firing.
With one session and a substituted model, the weaker results for Figures 10 and 11 cannot be attributed to any one of these.
What this attempt found about the dandiset¶
- It holds one of the twelve sessions analysed in the paper, and none from the two-workspace experiment, so four of the paper's main findings cannot be tested from it.
- The positions of the arm markers, which the paper's whole-arm model is built on, are not included; the file does hold joint angles and muscle lengths derived from them.
- The text descriptions of hand position and velocity give centimetres while the stated unit is metres. The stated unit, with the file's conversion factor, is the correct one.
- The dandiset's metadata on the archive names the paper without saying which part of it the data covers.