diff --git a/analysis/matlab/make_boxa_variations.m b/analysis/matlab/make_boxa_variations.m
new file mode 100644
index 0000000..5d131a1
--- /dev/null
+++ b/analysis/matlab/make_boxa_variations.m
@@ -0,0 +1,73 @@
+function make_boxa_variations()
+%MAKE_BOXA_VARIATIONS Variation folders for the two Box-A pooling groupings.
+%   MAKE_BOXA_VARIATIONS() writes variations/<grouping>_<window>/ -- each with
+%   data.csv + analyze.m (a copy of variation_analyze.m) + result.txt, fitting
+%   the paper LME (behavior ~ stim + day + stim:day + (1|rat)) on the success
+%   COUNT -- for:
+%     boxa_a2   B2           vs A2 + Box-A     ("b2 vs a2+a")
+%     boxa_b2   B2 + Box-A   vs A2             ("b2+a vs a2")
+%   over windows d0_5 (0-5), d6_10 (6-10), d0_10 (0-10). Box-A is the
+%   single-animal Electrode-Box-A condition, pooled into the control (boxa_a2)
+%   or the treatment (boxa_b2). Also writes variations/boxa_summary.csv with the
+%   residual / Satterthwaite / honest random-slope interaction p-values.
+
+thisDir = fileparts(mfilename('fullpath'));
+template = fullfile(thisDir, 'variation_analyze.m');
+root = fullfile(thisDir, 'variations');
+if ~exist(root, 'dir'); mkdir(root); end
+
+Tc = tdcs_load_data();
+B2 = 'Electrode-Box-B2'; A2 = 'Electrode-Box-A2'; BOXA = 'Electrode-Box-A';
+
+% grouping name | treatment groups (stim=1) | control groups (stim=0)
+groupings = {
+    'boxa_a2', {B2},        {A2, BOXA}
+    'boxa_b2', {B2, BOXA},  {A2}
+    };
+windows = {'d0_5', [0 5]; 'd6_10', [6 10]; 'd0_10', [0 10]};
+
+allGroup = cellstr(Tc.group);
+rows = {};
+for gi = 1:size(groupings, 1)
+    gname = groupings{gi, 1};
+    inTreat = ismember(allGroup, groupings{gi, 2});
+    inCtrl = ismember(allGroup, groupings{gi, 3});
+    for wi = 1:size(windows, 1)
+        wname = windows{wi, 1}; w = windows{wi, 2};
+        sel = (inTreat | inCtrl) & Tc.day >= w(1) & Tc.day <= w(2);
+        D = Tc(sel, :);
+        stim = double(inTreat(sel));
+        Dout = table(string(D.subject), string(D.group), D.day, D.success, ...
+            D.total, stim, 'VariableNames', ...
+            {'subject', 'group', 'day', 'success', 'total', 'stim'});
+
+        vname = [gname '_' wname];
+        folder = fullfile(root, vname);
+        if ~exist(folder, 'dir'); mkdir(folder); end
+        writetable(Dout, fullfile(folder, 'data.csv'));
+        copyfile(template, fullfile(folder, 'analyze.m'));
+
+        r = localRun(fullfile(folder, 'analyze.m'));
+        rows(end + 1, :) = {vname, gname, wname, r.nRats, r.nObs, ...
+            r.interP, r.interPsatt, r.interPrs, r.interEst, r.stimP}; %#ok<AGROW>
+        fprintf('  %-14s N=%d obs=%3d  int p: res=%.3f satt=%.3f rs=%-5s  est=%+.2f\n', ...
+            vname, r.nRats, r.nObs, r.interP, r.interPsatt, localNum(r.interPrs), r.interEst);
+    end
+end
+
+S = cell2table(rows, 'VariableNames', {'variation', 'grouping', 'window', ...
+    'nRats', 'nObs', 'interaction_p', 'interaction_p_satt', 'interaction_p_rs', ...
+    'interaction_est', 'stim_p'});
+writetable(S, fullfile(root, 'boxa_summary.csv'));
+fprintf('\nWrote %d folders + variations/boxa_summary.csv\n', size(rows, 1));
+
+end
+
+function r = localRun(scriptPath)
+run(scriptPath);
+r = VARRESULT; %#ok<NODEF>  (defined by the analyze.m script just run)
+end
+
+function s = localNum(x)
+if isnan(x); s = 'n/a'; else; s = sprintf('%.3f', x); end
+end
diff --git a/analysis/matlab/variations/boxa_a2_d0_10/analyze.m b/analysis/matlab/variations/boxa_a2_d0_10/analyze.m
new file mode 100644
index 0000000..82da5ed
--- /dev/null
+++ b/analysis/matlab/variations/boxa_a2_d0_10/analyze.m
@@ -0,0 +1,98 @@
+% Variation analysis -- the paper's linear mixed model on the successful-reach
+% COUNT, fit on this folder's curated data subset.
+%
+%   model: behavior ~ stim + day + stim:day + (1|rat)
+%     behavior = successful reaches (count per session)
+%     stim     = 1 for the treatment group(s), 0 for the control group(s)
+%     day      = training day within this window (0 = first analyzed day)
+%     rat      = subject (random intercept)
+%
+% Self-contained: reads data.csv beside this script and writes result.txt.
+% Run headless from this folder with:   matlab -batch "analyze"
+% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
+
+here = fileparts(mfilename('fullpath'));
+if isempty(here); here = pwd; end
+vname = regexprep(here, '.*[/\\]', '');            % folder name = variation id
+
+D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
+
+tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
+    'VariableNames', {'behavior', 'day', 'stim', 'rat'});
+m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
+C = m.Coefficients; A = anova(m); ci = coefCI(m);
+As = anova(m, 'DFMethod', 'satterthwaite');           % Satterthwaite denominator DF
+
+% Honest test: refit with a per-animal random SLOPE so the interaction DF
+% collapses toward the animal count (guarded -- may not converge in short windows).
+rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
+wst = warning('off', 'all');
+try
+    mr = fitlme(tbl, 'behavior ~ stim + day + stim:day + (day|rat)');
+    Ar = anova(mr, 'DFMethod', 'satterthwaite');
+    ri = strcmp(Ar.Term, 'day:stim');
+    rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
+catch
+end
+warning(wst);
+
+gi = @(t) find(strcmp(C.Name, t), 1);
+ga = @(t) find(strcmp(A.Term, t), 1);
+gs = @(t) find(strcmp(As.Term, t), 1);
+row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f  p=%.4g    p=%.4g (df=%.0f)\n', nm, ...
+    C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
+    As.pValue(gs(t)), As.DF2(gs(t)));
+
+maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
+maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
+if abs(maxT - maxC) > 2
+    cov = '** WARNING: unequal day coverage -- interaction may be confounded. **';
+else
+    cov = '(equal day coverage over this window)';
+end
+
+ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
+if pI >= 0.05
+    verdict = 'n.s. -- slopes parallel (no differential learning rate)';
+elseif eI > 0
+    verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
+else
+    verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
+end
+
+bar = repmat('=', 1, 78);
+raw = regexprep(evalc('disp(m)'), '</?strong>', '');
+s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
+s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)\n')];
+s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
+s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
+s = [s sprintf('control   (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
+s = [s sprintf('N = %d rats, %d sessions   raw day coverage: treat %d..%d, control %d..%d\n', ...
+    numel(unique(D.subject)), height(D), minT, maxT, minC, maxC)];
+s = [s sprintf('%s\n\n%s\nFULL MODEL SUMMARY -- fitlme\n%s\n%s\n', cov, bar, bar, raw)];
+s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (resid)', 'Satterthwaite: p (df)', repmat('-', 1, 76))];
+s = [s row('stim x day (interaction)', 'day:stim')];
+s = [s row('day (learning)', 'day')];
+s = [s row('stim (main, window start)', 'stim')];
+s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
+if rsOk
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
+else
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): did not converge for this window.\n')];
+end
+s = [s sprintf(['  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
+    '   model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
+s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
+s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
+
+fprintf('%s', s);
+fid = fopen(fullfile(here, 'result.txt'), 'w');
+fprintf(fid, '%s', s);
+fclose(fid);
+
+% Machine-readable handoff for the summary table (see make_variations.m).
+VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
+    'nObs', height(D), 'interP', pI, 'interEst', eI, ...
+    'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
+    'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
+    'covEqual', abs(maxT - maxC) <= 2);
diff --git a/analysis/matlab/variations/boxa_a2_d0_10/data.csv b/analysis/matlab/variations/boxa_a2_d0_10/data.csv
new file mode 100644
index 0000000..4974845
--- /dev/null
+++ b/analysis/matlab/variations/boxa_a2_d0_10/data.csv
@@ -0,0 +1,73 @@
+subject,group,day,success,total,stim
+Banh-mi-1,Electrode-Box-B2,0,13,84,1
+Banh-mi-1,Electrode-Box-B2,1,35,86,1
+Banh-mi-1,Electrode-Box-B2,2,47,110,1
+Banh-mi-1,Electrode-Box-B2,3,65,140,1
+Banh-mi-1,Electrode-Box-B2,4,70,127,1
+Banh-mi-1,Electrode-Box-B2,5,102,142,1
+Banh-mi-1,Electrode-Box-B2,6,90,131,1
+Banh-mi-1,Electrode-Box-B2,7,109,148,1
+Banh-mi-1,Electrode-Box-B2,8,104,137,1
+Banh-mi-1,Electrode-Box-B2,9,119,150,1
+Banh-mi-1,Electrode-Box-B2,10,121,158,1
+Banh-mi-2,Electrode-Box-A2,0,25,97,0
+Banh-mi-2,Electrode-Box-A2,1,22,101,0
+Banh-mi-2,Electrode-Box-A2,2,22,119,0
+Banh-mi-2,Electrode-Box-A2,3,29,118,0
+Banh-mi-2,Electrode-Box-A2,4,27,136,0
+Banh-mi-2,Electrode-Box-A2,5,43,146,0
+Banh-mi-2,Electrode-Box-A2,6,74,146,0
+Banh-mi-2,Electrode-Box-A2,7,70,148,0
+Banh-mi-2,Electrode-Box-A2,8,65,130,0
+Banh-mi-2,Electrode-Box-A2,9,79,151,0
+Banh-mi-2,Electrode-Box-A2,10,93,152,0
+Egg-tart-1,Electrode-Box-B2,0,7,56,1
+Egg-tart-1,Electrode-Box-B2,1,16,78,1
+Egg-tart-1,Electrode-Box-B2,2,23,103,1
+Egg-tart-1,Electrode-Box-B2,3,63,120,1
+Egg-tart-1,Electrode-Box-B2,4,69,132,1
+Egg-tart-1,Electrode-Box-B2,5,83,136,1
+Egg-tart-1,Electrode-Box-B2,6,71,142,1
+Egg-tart-1,Electrode-Box-B2,7,79,138,1
+Egg-tart-1,Electrode-Box-B2,8,98,142,1
+Egg-tart-1,Electrode-Box-B2,9,89,139,1
+Egg-tart-1,Electrode-Box-B2,10,96,143,1
+Egg-tart-2,Electrode-Box-A2,0,9,32,0
+Egg-tart-2,Electrode-Box-A2,1,2,38,0
+Egg-tart-2,Electrode-Box-A2,2,31,93,0
+Egg-tart-2,Electrode-Box-A2,3,44,101,0
+Egg-tart-2,Electrode-Box-A2,4,54,131,0
+Egg-tart-2,Electrode-Box-A2,5,84,139,0
+Egg-tart-2,Electrode-Box-A2,6,85,145,0
+Egg-tart-2,Electrode-Box-A2,7,79,143,0
+Egg-tart-2,Electrode-Box-A2,8,76,131,0
+Egg-tart-2,Electrode-Box-A2,9,88,149,0
+Egg-tart-2,Electrode-Box-A2,10,81,151,0
+Khoai-tay-1,Electrode-Box-A,0,14,62,0
+Khoai-tay-1,Electrode-Box-A,1,22,83,0
+Khoai-tay-1,Electrode-Box-A,2,6,79,0
+Khoai-tay-1,Electrode-Box-A,3,28,97,0
+Khoai-tay-1,Electrode-Box-A,4,60,134,0
+Khoai-tay-1,Electrode-Box-A,5,75,138,0
+Khoai-tay-1,Electrode-Box-A,6,81,137,0
+Khoai-tay-1,Electrode-Box-A,7,78,147,0
+Khoai-tay-1,Electrode-Box-A,8,91,132,0
+Khoai-tay-1,Electrode-Box-A,9,99,146,0
+Khoai-tay-1,Electrode-Box-A,10,94,143,0
+Root-beer-1,Electrode-Box-B2,0,11,85,1
+Root-beer-1,Electrode-Box-B2,1,18,76,1
+Root-beer-1,Electrode-Box-B2,2,40,105,1
+Root-beer-1,Electrode-Box-B2,3,55,134,1
+Root-beer-1,Electrode-Box-B2,4,75,136,1
+Root-beer-1,Electrode-Box-B2,5,64,133,1
+Root-beer-1,Electrode-Box-B2,6,104,139,1
+Root-beer-1,Electrode-Box-B2,7,98,148,1
+Root-beer-1,Electrode-Box-B2,8,81,145,1
+Root-beer-1,Electrode-Box-B2,9,89,156,1
+Root-beer-1,Electrode-Box-B2,10,105,158,1
+Root-beer-2,Electrode-Box-A2,0,22,74,0
+Root-beer-2,Electrode-Box-A2,1,31,87,0
+Root-beer-2,Electrode-Box-A2,2,49,134,0
+Root-beer-2,Electrode-Box-A2,3,31,89,0
+Root-beer-2,Electrode-Box-A2,4,60,140,0
+Root-beer-2,Electrode-Box-A2,5,84,147,0
diff --git a/analysis/matlab/variations/boxa_a2_d0_10/result.txt b/analysis/matlab/variations/boxa_a2_d0_10/result.txt
new file mode 100644
index 0000000..9d79605
--- /dev/null
+++ b/analysis/matlab/variations/boxa_a2_d0_10/result.txt
@@ -0,0 +1,68 @@
+==============================================================================
+VARIATION: boxa_a2_d0_10
+==============================================================================
+model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)
+day = training day within window (0 = first analyzed day)
+treatment (stim=1): Electrode-Box-B2
+control   (stim=0): Electrode-Box-A, Electrode-Box-A2
+N = 7 rats, 72 sessions   raw day coverage: treat 0..10, control 0..10
+(equal day coverage over this window)
+
+==============================================================================
+FULL MODEL SUMMARY -- fitlme
+==============================================================================
+
+Linear mixed-effects model fit by ML
+
+Model information:
+    Number of observations              72
+    Fixed effects coefficients           4
+    Random effects coefficients          7
+    Covariance parameters                2
+
+Formula:
+    behavior ~ 1 + day*stim + (1 | rat)
+
+Model fit statistics:
+    AIC      BIC       LogLikelihood    Deviance
+    590.8    604.46    -289.4           578.8   
+
+Fixed effects coefficients (95% CIs):
+    Name                   Estimate    SE         tStat     DF    pValue    
+    {'(Intercept)'}        15.599       3.8531    4.0485    68    0.00013444
+    {'day'        }        8.3258      0.69175    12.036    68    1.6397e-18
+    {'stim'       }        6.7494       5.8389    1.1559    68       0.25175
+    {'day:stim'   }        1.1985       1.0141    1.1818    68        0.2414
+
+
+    Lower       Upper 
+      7.9104    23.288
+      6.9454    9.7061
+     -4.9019    18.401
+    -0.82514    3.2221
+
+Random effects covariance parameters (95% CIs):
+Group: rat (7 Levels)
+    Name1                  Name2                  Type           Estimate
+    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0       
+
+
+    Lower    Upper
+    NaN      NaN  
+
+Group: Error
+    Name               Estimate    Lower     Upper 
+    {'Res Std'}        13.471      11.441    15.861
+
+
+effect                     t(df) / F(df1)     p (resid)    Satterthwaite: p (df)
+----------------------------------------------------------------------------
+stim x day (interaction)   t(68)=  1.18 F(1)=  1.397  p=0.2414    p=0.2412 (df=72)
+day (learning)             t(68)= 12.04 F(1)=144.862  p=1.64e-18    p=6.575e-19 (df=72)
+stim (main, window start)  t(68)=  1.16 F(1)=  1.336  p=0.2517    p=0.2515 (df=72)
+interaction 95% CI: [-0.83, +3.22]
+HONEST LME (per-animal random slope, day|rat): interaction F(1,14.0)=0.68, p=0.4245
+  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope
+   model above is the honest learning-rate test -- DF collapses toward the animal count.)
+INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.2414, slope diff=+1.20)
+Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
diff --git a/analysis/matlab/variations/boxa_a2_d0_5/analyze.m b/analysis/matlab/variations/boxa_a2_d0_5/analyze.m
new file mode 100644
index 0000000..82da5ed
--- /dev/null
+++ b/analysis/matlab/variations/boxa_a2_d0_5/analyze.m
@@ -0,0 +1,98 @@
+% Variation analysis -- the paper's linear mixed model on the successful-reach
+% COUNT, fit on this folder's curated data subset.
+%
+%   model: behavior ~ stim + day + stim:day + (1|rat)
+%     behavior = successful reaches (count per session)
+%     stim     = 1 for the treatment group(s), 0 for the control group(s)
+%     day      = training day within this window (0 = first analyzed day)
+%     rat      = subject (random intercept)
+%
+% Self-contained: reads data.csv beside this script and writes result.txt.
+% Run headless from this folder with:   matlab -batch "analyze"
+% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
+
+here = fileparts(mfilename('fullpath'));
+if isempty(here); here = pwd; end
+vname = regexprep(here, '.*[/\\]', '');            % folder name = variation id
+
+D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
+
+tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
+    'VariableNames', {'behavior', 'day', 'stim', 'rat'});
+m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
+C = m.Coefficients; A = anova(m); ci = coefCI(m);
+As = anova(m, 'DFMethod', 'satterthwaite');           % Satterthwaite denominator DF
+
+% Honest test: refit with a per-animal random SLOPE so the interaction DF
+% collapses toward the animal count (guarded -- may not converge in short windows).
+rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
+wst = warning('off', 'all');
+try
+    mr = fitlme(tbl, 'behavior ~ stim + day + stim:day + (day|rat)');
+    Ar = anova(mr, 'DFMethod', 'satterthwaite');
+    ri = strcmp(Ar.Term, 'day:stim');
+    rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
+catch
+end
+warning(wst);
+
+gi = @(t) find(strcmp(C.Name, t), 1);
+ga = @(t) find(strcmp(A.Term, t), 1);
+gs = @(t) find(strcmp(As.Term, t), 1);
+row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f  p=%.4g    p=%.4g (df=%.0f)\n', nm, ...
+    C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
+    As.pValue(gs(t)), As.DF2(gs(t)));
+
+maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
+maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
+if abs(maxT - maxC) > 2
+    cov = '** WARNING: unequal day coverage -- interaction may be confounded. **';
+else
+    cov = '(equal day coverage over this window)';
+end
+
+ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
+if pI >= 0.05
+    verdict = 'n.s. -- slopes parallel (no differential learning rate)';
+elseif eI > 0
+    verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
+else
+    verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
+end
+
+bar = repmat('=', 1, 78);
+raw = regexprep(evalc('disp(m)'), '</?strong>', '');
+s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
+s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)\n')];
+s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
+s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
+s = [s sprintf('control   (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
+s = [s sprintf('N = %d rats, %d sessions   raw day coverage: treat %d..%d, control %d..%d\n', ...
+    numel(unique(D.subject)), height(D), minT, maxT, minC, maxC)];
+s = [s sprintf('%s\n\n%s\nFULL MODEL SUMMARY -- fitlme\n%s\n%s\n', cov, bar, bar, raw)];
+s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (resid)', 'Satterthwaite: p (df)', repmat('-', 1, 76))];
+s = [s row('stim x day (interaction)', 'day:stim')];
+s = [s row('day (learning)', 'day')];
+s = [s row('stim (main, window start)', 'stim')];
+s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
+if rsOk
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
+else
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): did not converge for this window.\n')];
+end
+s = [s sprintf(['  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
+    '   model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
+s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
+s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
+
+fprintf('%s', s);
+fid = fopen(fullfile(here, 'result.txt'), 'w');
+fprintf(fid, '%s', s);
+fclose(fid);
+
+% Machine-readable handoff for the summary table (see make_variations.m).
+VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
+    'nObs', height(D), 'interP', pI, 'interEst', eI, ...
+    'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
+    'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
+    'covEqual', abs(maxT - maxC) <= 2);
diff --git a/analysis/matlab/variations/boxa_a2_d0_5/data.csv b/analysis/matlab/variations/boxa_a2_d0_5/data.csv
new file mode 100644
index 0000000..4090301
--- /dev/null
+++ b/analysis/matlab/variations/boxa_a2_d0_5/data.csv
@@ -0,0 +1,43 @@
+subject,group,day,success,total,stim
+Banh-mi-1,Electrode-Box-B2,0,13,84,1
+Banh-mi-1,Electrode-Box-B2,1,35,86,1
+Banh-mi-1,Electrode-Box-B2,2,47,110,1
+Banh-mi-1,Electrode-Box-B2,3,65,140,1
+Banh-mi-1,Electrode-Box-B2,4,70,127,1
+Banh-mi-1,Electrode-Box-B2,5,102,142,1
+Banh-mi-2,Electrode-Box-A2,0,25,97,0
+Banh-mi-2,Electrode-Box-A2,1,22,101,0
+Banh-mi-2,Electrode-Box-A2,2,22,119,0
+Banh-mi-2,Electrode-Box-A2,3,29,118,0
+Banh-mi-2,Electrode-Box-A2,4,27,136,0
+Banh-mi-2,Electrode-Box-A2,5,43,146,0
+Egg-tart-1,Electrode-Box-B2,0,7,56,1
+Egg-tart-1,Electrode-Box-B2,1,16,78,1
+Egg-tart-1,Electrode-Box-B2,2,23,103,1
+Egg-tart-1,Electrode-Box-B2,3,63,120,1
+Egg-tart-1,Electrode-Box-B2,4,69,132,1
+Egg-tart-1,Electrode-Box-B2,5,83,136,1
+Egg-tart-2,Electrode-Box-A2,0,9,32,0
+Egg-tart-2,Electrode-Box-A2,1,2,38,0
+Egg-tart-2,Electrode-Box-A2,2,31,93,0
+Egg-tart-2,Electrode-Box-A2,3,44,101,0
+Egg-tart-2,Electrode-Box-A2,4,54,131,0
+Egg-tart-2,Electrode-Box-A2,5,84,139,0
+Khoai-tay-1,Electrode-Box-A,0,14,62,0
+Khoai-tay-1,Electrode-Box-A,1,22,83,0
+Khoai-tay-1,Electrode-Box-A,2,6,79,0
+Khoai-tay-1,Electrode-Box-A,3,28,97,0
+Khoai-tay-1,Electrode-Box-A,4,60,134,0
+Khoai-tay-1,Electrode-Box-A,5,75,138,0
+Root-beer-1,Electrode-Box-B2,0,11,85,1
+Root-beer-1,Electrode-Box-B2,1,18,76,1
+Root-beer-1,Electrode-Box-B2,2,40,105,1
+Root-beer-1,Electrode-Box-B2,3,55,134,1
+Root-beer-1,Electrode-Box-B2,4,75,136,1
+Root-beer-1,Electrode-Box-B2,5,64,133,1
+Root-beer-2,Electrode-Box-A2,0,22,74,0
+Root-beer-2,Electrode-Box-A2,1,31,87,0
+Root-beer-2,Electrode-Box-A2,2,49,134,0
+Root-beer-2,Electrode-Box-A2,3,31,89,0
+Root-beer-2,Electrode-Box-A2,4,60,140,0
+Root-beer-2,Electrode-Box-A2,5,84,147,0
diff --git a/analysis/matlab/variations/boxa_a2_d0_5/result.txt b/analysis/matlab/variations/boxa_a2_d0_5/result.txt
new file mode 100644
index 0000000..e88076d
--- /dev/null
+++ b/analysis/matlab/variations/boxa_a2_d0_5/result.txt
@@ -0,0 +1,68 @@
+==============================================================================
+VARIATION: boxa_a2_d0_5
+==============================================================================
+model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)
+day = training day within window (0 = first analyzed day)
+treatment (stim=1): Electrode-Box-B2
+control   (stim=0): Electrode-Box-A, Electrode-Box-A2
+N = 7 rats, 42 sessions   raw day coverage: treat 0..5, control 0..5
+(equal day coverage over this window)
+
+==============================================================================
+FULL MODEL SUMMARY -- fitlme
+==============================================================================
+
+Linear mixed-effects model fit by ML
+
+Model information:
+    Number of observations              42
+    Fixed effects coefficients           4
+    Random effects coefficients          7
+    Covariance parameters                2
+
+Formula:
+    behavior ~ 1 + day*stim + (1 | rat)
+
+Model fit statistics:
+    AIC       BIC       LogLikelihood    Deviance
+    339.96    350.39    -163.98          327.96  
+
+Fixed effects coefficients (95% CIs):
+    Name                   Estimate    SE        tStat        DF    pValue    
+    {'(Intercept)'}           10.06    4.3445       2.3155    38      0.026083
+    {'day'        }          10.543    1.4349       7.3472    38    8.3775e-09
+    {'stim'       }        -0.55159    6.6363    -0.083116    38        0.9342
+    {'day:stim'   }          4.6762    2.1919       2.1334    38       0.03941
+
+
+    Lower      Upper 
+     1.2645    18.855
+      7.638    13.448
+    -13.986    12.883
+     0.2389    9.1135
+
+Random effects covariance parameters (95% CIs):
+Group: rat (7 Levels)
+    Name1                  Name2                  Type           Estimate
+    {'(Intercept)'}        {'(Intercept)'}        {'std'}        0       
+
+
+    Lower    Upper
+    NaN      NaN  
+
+Group: Error
+    Name               Estimate    Lower     Upper 
+    {'Res Std'}        12.006      9.6941    14.868
+
+
+effect                     t(df) / F(df1)     p (resid)    Satterthwaite: p (df)
+----------------------------------------------------------------------------
+stim x day (interaction)   t(38)=  2.13 F(1)=  4.551  p=0.03941    p=0.03878 (df=42)
+day (learning)             t(38)=  7.35 F(1)= 53.982  p=8.377e-09    p=4.651e-09 (df=42)
+stim (main, window start)  t(38)= -0.08 F(1)=  0.007  p=0.9342    p=0.9342 (df=42)
+interaction 95% CI: [+0.24, +9.11]
+HONEST LME (per-animal random slope, day|rat): interaction F(1,7.0)=2.87, p=0.1341
+  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope
+   model above is the honest learning-rate test -- DF collapses toward the animal count.)
+INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.03941, slope diff=+4.68)
+Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
diff --git a/analysis/matlab/variations/boxa_a2_d6_10/analyze.m b/analysis/matlab/variations/boxa_a2_d6_10/analyze.m
new file mode 100644
index 0000000..82da5ed
--- /dev/null
+++ b/analysis/matlab/variations/boxa_a2_d6_10/analyze.m
@@ -0,0 +1,98 @@
+% Variation analysis -- the paper's linear mixed model on the successful-reach
+% COUNT, fit on this folder's curated data subset.
+%
+%   model: behavior ~ stim + day + stim:day + (1|rat)
+%     behavior = successful reaches (count per session)
+%     stim     = 1 for the treatment group(s), 0 for the control group(s)
+%     day      = training day within this window (0 = first analyzed day)
+%     rat      = subject (random intercept)
+%
+% Self-contained: reads data.csv beside this script and writes result.txt.
+% Run headless from this folder with:   matlab -batch "analyze"
+% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
+
+here = fileparts(mfilename('fullpath'));
+if isempty(here); here = pwd; end
+vname = regexprep(here, '.*[/\\]', '');            % folder name = variation id
+
+D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
+
+tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
+    'VariableNames', {'behavior', 'day', 'stim', 'rat'});
+m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
+C = m.Coefficients; A = anova(m); ci = coefCI(m);
+As = anova(m, 'DFMethod', 'satterthwaite');           % Satterthwaite denominator DF
+
+% Honest test: refit with a per-animal random SLOPE so the interaction DF
+% collapses toward the animal count (guarded -- may not converge in short windows).
+rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
+wst = warning('off', 'all');
+try
+    mr = fitlme(tbl, 'behavior ~ stim + day + stim:day + (day|rat)');
+    Ar = anova(mr, 'DFMethod', 'satterthwaite');
+    ri = strcmp(Ar.Term, 'day:stim');
+    rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
+catch
+end
+warning(wst);
+
+gi = @(t) find(strcmp(C.Name, t), 1);
+ga = @(t) find(strcmp(A.Term, t), 1);
+gs = @(t) find(strcmp(As.Term, t), 1);
+row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f  p=%.4g    p=%.4g (df=%.0f)\n', nm, ...
+    C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
+    As.pValue(gs(t)), As.DF2(gs(t)));
+
+maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
+maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
+if abs(maxT - maxC) > 2
+    cov = '** WARNING: unequal day coverage -- interaction may be confounded. **';
+else
+    cov = '(equal day coverage over this window)';
+end
+
+ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
+if pI >= 0.05
+    verdict = 'n.s. -- slopes parallel (no differential learning rate)';
+elseif eI > 0
+    verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
+else
+    verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
+end
+
+bar = repmat('=', 1, 78);
+raw = regexprep(evalc('disp(m)'), '</?strong>', '');
+s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
+s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)\n')];
+s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
+s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
+s = [s sprintf('control   (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
+s = [s sprintf('N = %d rats, %d sessions   raw day coverage: treat %d..%d, control %d..%d\n', ...
+    numel(unique(D.subject)), height(D), minT, maxT, minC, maxC)];
+s = [s sprintf('%s\n\n%s\nFULL MODEL SUMMARY -- fitlme\n%s\n%s\n', cov, bar, bar, raw)];
+s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (resid)', 'Satterthwaite: p (df)', repmat('-', 1, 76))];
+s = [s row('stim x day (interaction)', 'day:stim')];
+s = [s row('day (learning)', 'day')];
+s = [s row('stim (main, window start)', 'stim')];
+s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
+if rsOk
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
+else
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): did not converge for this window.\n')];
+end
+s = [s sprintf(['  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
+    '   model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
+s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
+s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
+
+fprintf('%s', s);
+fid = fopen(fullfile(here, 'result.txt'), 'w');
+fprintf(fid, '%s', s);
+fclose(fid);
+
+% Machine-readable handoff for the summary table (see make_variations.m).
+VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
+    'nObs', height(D), 'interP', pI, 'interEst', eI, ...
+    'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
+    'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
+    'covEqual', abs(maxT - maxC) <= 2);
diff --git a/analysis/matlab/variations/boxa_a2_d6_10/data.csv b/analysis/matlab/variations/boxa_a2_d6_10/data.csv
new file mode 100644
index 0000000..b12c77a
--- /dev/null
+++ b/analysis/matlab/variations/boxa_a2_d6_10/data.csv
@@ -0,0 +1,31 @@
+subject,group,day,success,total,stim
+Banh-mi-1,Electrode-Box-B2,6,90,131,1
+Banh-mi-1,Electrode-Box-B2,7,109,148,1
+Banh-mi-1,Electrode-Box-B2,8,104,137,1
+Banh-mi-1,Electrode-Box-B2,9,119,150,1
+Banh-mi-1,Electrode-Box-B2,10,121,158,1
+Banh-mi-2,Electrode-Box-A2,6,74,146,0
+Banh-mi-2,Electrode-Box-A2,7,70,148,0
+Banh-mi-2,Electrode-Box-A2,8,65,130,0
+Banh-mi-2,Electrode-Box-A2,9,79,151,0
+Banh-mi-2,Electrode-Box-A2,10,93,152,0
+Egg-tart-1,Electrode-Box-B2,6,71,142,1
+Egg-tart-1,Electrode-Box-B2,7,79,138,1
+Egg-tart-1,Electrode-Box-B2,8,98,142,1
+Egg-tart-1,Electrode-Box-B2,9,89,139,1
+Egg-tart-1,Electrode-Box-B2,10,96,143,1
+Egg-tart-2,Electrode-Box-A2,6,85,145,0
+Egg-tart-2,Electrode-Box-A2,7,79,143,0
+Egg-tart-2,Electrode-Box-A2,8,76,131,0
+Egg-tart-2,Electrode-Box-A2,9,88,149,0
+Egg-tart-2,Electrode-Box-A2,10,81,151,0
+Khoai-tay-1,Electrode-Box-A,6,81,137,0
+Khoai-tay-1,Electrode-Box-A,7,78,147,0
+Khoai-tay-1,Electrode-Box-A,8,91,132,0
+Khoai-tay-1,Electrode-Box-A,9,99,146,0
+Khoai-tay-1,Electrode-Box-A,10,94,143,0
+Root-beer-1,Electrode-Box-B2,6,104,139,1
+Root-beer-1,Electrode-Box-B2,7,98,148,1
+Root-beer-1,Electrode-Box-B2,8,81,145,1
+Root-beer-1,Electrode-Box-B2,9,89,156,1
+Root-beer-1,Electrode-Box-B2,10,105,158,1
diff --git a/analysis/matlab/variations/boxa_a2_d6_10/result.txt b/analysis/matlab/variations/boxa_a2_d6_10/result.txt
new file mode 100644
index 0000000..ff825b2
--- /dev/null
+++ b/analysis/matlab/variations/boxa_a2_d6_10/result.txt
@@ -0,0 +1,68 @@
+==============================================================================
+VARIATION: boxa_a2_d6_10
+==============================================================================
+model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)
+day = training day within window (0 = first analyzed day)
+treatment (stim=1): Electrode-Box-B2
+control   (stim=0): Electrode-Box-A, Electrode-Box-A2
+N = 6 rats, 30 sessions   raw day coverage: treat 6..10, control 6..10
+(equal day coverage over this window)
+
+==============================================================================
+FULL MODEL SUMMARY -- fitlme
+==============================================================================
+
+Linear mixed-effects model fit by ML
+
+Model information:
+    Number of observations              30
+    Fixed effects coefficients           4
+    Random effects coefficients          6
+    Covariance parameters                2
+
+Formula:
+    behavior ~ 1 + day*stim + (1 | rat)
+
+Model fit statistics:
+    AIC       BIC       LogLikelihood    Deviance
+    231.74    240.15    -109.87          219.74  
+
+Fixed effects coefficients (95% CIs):
+    Name                   Estimate    SE        tStat      DF    pValue    
+    {'(Intercept)'}        75.867      5.1856      14.63    26    4.6141e-14
+    {'day'        }        3.1667      1.4985     2.1132    26      0.044334
+    {'stim'       }        12.667      7.3335     1.7272    26      0.095989
+    {'day:stim'   }             1      2.1192    0.47188    26       0.64095
+
+
+    Lower       Upper 
+      65.208    86.526
+    0.086483    6.2469
+     -2.4075    27.741
+      -3.356     5.356
+
+Random effects covariance parameters (95% CIs):
+Group: rat (6 Levels)
+    Name1                  Name2                  Type           Estimate
+    {'(Intercept)'}        {'(Intercept)'}        {'std'}        6.3444  
+
+
+    Lower     Upper 
+    2.9639    13.581
+
+Group: Error
+    Name               Estimate    Lower     Upper 
+    {'Res Std'}        8.2076      6.1852    10.891
+
+
+effect                     t(df) / F(df1)     p (resid)    Satterthwaite: p (df)
+----------------------------------------------------------------------------
+stim x day (interaction)   t(26)=  0.47 F(1)=  0.223  p=0.6409    p=0.6413 (df=24)
+day (learning)             t(26)=  2.11 F(1)=  4.466  p=0.04433    p=0.04517 (df=24)
+stim (main, window start)  t(26)=  1.73 F(1)=  2.983  p=0.09599    p=0.1083 (df=13)
+interaction 95% CI: [-3.36, +5.36]
+HONEST LME (per-animal random slope, day|rat): interaction F(1,6.0)=0.18, p=0.6872
+  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope
+   model above is the honest learning-rate test -- DF collapses toward the animal count.)
+INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.6409, slope diff=+1.00)
+Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
diff --git a/analysis/matlab/variations/boxa_b2_d0_10/analyze.m b/analysis/matlab/variations/boxa_b2_d0_10/analyze.m
new file mode 100644
index 0000000..82da5ed
--- /dev/null
+++ b/analysis/matlab/variations/boxa_b2_d0_10/analyze.m
@@ -0,0 +1,98 @@
+% Variation analysis -- the paper's linear mixed model on the successful-reach
+% COUNT, fit on this folder's curated data subset.
+%
+%   model: behavior ~ stim + day + stim:day + (1|rat)
+%     behavior = successful reaches (count per session)
+%     stim     = 1 for the treatment group(s), 0 for the control group(s)
+%     day      = training day within this window (0 = first analyzed day)
+%     rat      = subject (random intercept)
+%
+% Self-contained: reads data.csv beside this script and writes result.txt.
+% Run headless from this folder with:   matlab -batch "analyze"
+% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
+
+here = fileparts(mfilename('fullpath'));
+if isempty(here); here = pwd; end
+vname = regexprep(here, '.*[/\\]', '');            % folder name = variation id
+
+D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
+
+tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
+    'VariableNames', {'behavior', 'day', 'stim', 'rat'});
+m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
+C = m.Coefficients; A = anova(m); ci = coefCI(m);
+As = anova(m, 'DFMethod', 'satterthwaite');           % Satterthwaite denominator DF
+
+% Honest test: refit with a per-animal random SLOPE so the interaction DF
+% collapses toward the animal count (guarded -- may not converge in short windows).
+rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
+wst = warning('off', 'all');
+try
+    mr = fitlme(tbl, 'behavior ~ stim + day + stim:day + (day|rat)');
+    Ar = anova(mr, 'DFMethod', 'satterthwaite');
+    ri = strcmp(Ar.Term, 'day:stim');
+    rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
+catch
+end
+warning(wst);
+
+gi = @(t) find(strcmp(C.Name, t), 1);
+ga = @(t) find(strcmp(A.Term, t), 1);
+gs = @(t) find(strcmp(As.Term, t), 1);
+row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f  p=%.4g    p=%.4g (df=%.0f)\n', nm, ...
+    C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
+    As.pValue(gs(t)), As.DF2(gs(t)));
+
+maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
+maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
+if abs(maxT - maxC) > 2
+    cov = '** WARNING: unequal day coverage -- interaction may be confounded. **';
+else
+    cov = '(equal day coverage over this window)';
+end
+
+ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
+if pI >= 0.05
+    verdict = 'n.s. -- slopes parallel (no differential learning rate)';
+elseif eI > 0
+    verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
+else
+    verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
+end
+
+bar = repmat('=', 1, 78);
+raw = regexprep(evalc('disp(m)'), '</?strong>', '');
+s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
+s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)\n')];
+s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
+s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
+s = [s sprintf('control   (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
+s = [s sprintf('N = %d rats, %d sessions   raw day coverage: treat %d..%d, control %d..%d\n', ...
+    numel(unique(D.subject)), height(D), minT, maxT, minC, maxC)];
+s = [s sprintf('%s\n\n%s\nFULL MODEL SUMMARY -- fitlme\n%s\n%s\n', cov, bar, bar, raw)];
+s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (resid)', 'Satterthwaite: p (df)', repmat('-', 1, 76))];
+s = [s row('stim x day (interaction)', 'day:stim')];
+s = [s row('day (learning)', 'day')];
+s = [s row('stim (main, window start)', 'stim')];
+s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
+if rsOk
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
+else
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): did not converge for this window.\n')];
+end
+s = [s sprintf(['  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
+    '   model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
+s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
+s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
+
+fprintf('%s', s);
+fid = fopen(fullfile(here, 'result.txt'), 'w');
+fprintf(fid, '%s', s);
+fclose(fid);
+
+% Machine-readable handoff for the summary table (see make_variations.m).
+VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
+    'nObs', height(D), 'interP', pI, 'interEst', eI, ...
+    'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
+    'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
+    'covEqual', abs(maxT - maxC) <= 2);
diff --git a/analysis/matlab/variations/boxa_b2_d0_10/data.csv b/analysis/matlab/variations/boxa_b2_d0_10/data.csv
new file mode 100644
index 0000000..dd7e209
--- /dev/null
+++ b/analysis/matlab/variations/boxa_b2_d0_10/data.csv
@@ -0,0 +1,73 @@
+subject,group,day,success,total,stim
+Banh-mi-1,Electrode-Box-B2,0,13,84,1
+Banh-mi-1,Electrode-Box-B2,1,35,86,1
+Banh-mi-1,Electrode-Box-B2,2,47,110,1
+Banh-mi-1,Electrode-Box-B2,3,65,140,1
+Banh-mi-1,Electrode-Box-B2,4,70,127,1
+Banh-mi-1,Electrode-Box-B2,5,102,142,1
+Banh-mi-1,Electrode-Box-B2,6,90,131,1
+Banh-mi-1,Electrode-Box-B2,7,109,148,1
+Banh-mi-1,Electrode-Box-B2,8,104,137,1
+Banh-mi-1,Electrode-Box-B2,9,119,150,1
+Banh-mi-1,Electrode-Box-B2,10,121,158,1
+Banh-mi-2,Electrode-Box-A2,0,25,97,0
+Banh-mi-2,Electrode-Box-A2,1,22,101,0
+Banh-mi-2,Electrode-Box-A2,2,22,119,0
+Banh-mi-2,Electrode-Box-A2,3,29,118,0
+Banh-mi-2,Electrode-Box-A2,4,27,136,0
+Banh-mi-2,Electrode-Box-A2,5,43,146,0
+Banh-mi-2,Electrode-Box-A2,6,74,146,0
+Banh-mi-2,Electrode-Box-A2,7,70,148,0
+Banh-mi-2,Electrode-Box-A2,8,65,130,0
+Banh-mi-2,Electrode-Box-A2,9,79,151,0
+Banh-mi-2,Electrode-Box-A2,10,93,152,0
+Egg-tart-1,Electrode-Box-B2,0,7,56,1
+Egg-tart-1,Electrode-Box-B2,1,16,78,1
+Egg-tart-1,Electrode-Box-B2,2,23,103,1
+Egg-tart-1,Electrode-Box-B2,3,63,120,1
+Egg-tart-1,Electrode-Box-B2,4,69,132,1
+Egg-tart-1,Electrode-Box-B2,5,83,136,1
+Egg-tart-1,Electrode-Box-B2,6,71,142,1
+Egg-tart-1,Electrode-Box-B2,7,79,138,1
+Egg-tart-1,Electrode-Box-B2,8,98,142,1
+Egg-tart-1,Electrode-Box-B2,9,89,139,1
+Egg-tart-1,Electrode-Box-B2,10,96,143,1
+Egg-tart-2,Electrode-Box-A2,0,9,32,0
+Egg-tart-2,Electrode-Box-A2,1,2,38,0
+Egg-tart-2,Electrode-Box-A2,2,31,93,0
+Egg-tart-2,Electrode-Box-A2,3,44,101,0
+Egg-tart-2,Electrode-Box-A2,4,54,131,0
+Egg-tart-2,Electrode-Box-A2,5,84,139,0
+Egg-tart-2,Electrode-Box-A2,6,85,145,0
+Egg-tart-2,Electrode-Box-A2,7,79,143,0
+Egg-tart-2,Electrode-Box-A2,8,76,131,0
+Egg-tart-2,Electrode-Box-A2,9,88,149,0
+Egg-tart-2,Electrode-Box-A2,10,81,151,0
+Khoai-tay-1,Electrode-Box-A,0,14,62,1
+Khoai-tay-1,Electrode-Box-A,1,22,83,1
+Khoai-tay-1,Electrode-Box-A,2,6,79,1
+Khoai-tay-1,Electrode-Box-A,3,28,97,1
+Khoai-tay-1,Electrode-Box-A,4,60,134,1
+Khoai-tay-1,Electrode-Box-A,5,75,138,1
+Khoai-tay-1,Electrode-Box-A,6,81,137,1
+Khoai-tay-1,Electrode-Box-A,7,78,147,1
+Khoai-tay-1,Electrode-Box-A,8,91,132,1
+Khoai-tay-1,Electrode-Box-A,9,99,146,1
+Khoai-tay-1,Electrode-Box-A,10,94,143,1
+Root-beer-1,Electrode-Box-B2,0,11,85,1
+Root-beer-1,Electrode-Box-B2,1,18,76,1
+Root-beer-1,Electrode-Box-B2,2,40,105,1
+Root-beer-1,Electrode-Box-B2,3,55,134,1
+Root-beer-1,Electrode-Box-B2,4,75,136,1
+Root-beer-1,Electrode-Box-B2,5,64,133,1
+Root-beer-1,Electrode-Box-B2,6,104,139,1
+Root-beer-1,Electrode-Box-B2,7,98,148,1
+Root-beer-1,Electrode-Box-B2,8,81,145,1
+Root-beer-1,Electrode-Box-B2,9,89,156,1
+Root-beer-1,Electrode-Box-B2,10,105,158,1
+Root-beer-2,Electrode-Box-A2,0,22,74,0
+Root-beer-2,Electrode-Box-A2,1,31,87,0
+Root-beer-2,Electrode-Box-A2,2,49,134,0
+Root-beer-2,Electrode-Box-A2,3,31,89,0
+Root-beer-2,Electrode-Box-A2,4,60,140,0
+Root-beer-2,Electrode-Box-A2,5,84,147,0
diff --git a/analysis/matlab/variations/boxa_b2_d0_10/result.txt b/analysis/matlab/variations/boxa_b2_d0_10/result.txt
new file mode 100644
index 0000000..0e12cda
--- /dev/null
+++ b/analysis/matlab/variations/boxa_b2_d0_10/result.txt
@@ -0,0 +1,68 @@
+==============================================================================
+VARIATION: boxa_b2_d0_10
+==============================================================================
+model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)
+day = training day within window (0 = first analyzed day)
+treatment (stim=1): Electrode-Box-A, Electrode-Box-B2
+control   (stim=0): Electrode-Box-A2
+N = 7 rats, 72 sessions   raw day coverage: treat 0..10, control 0..10
+(equal day coverage over this window)
+
+==============================================================================
+FULL MODEL SUMMARY -- fitlme
+==============================================================================
+
+Linear mixed-effects model fit by ML
+
+Model information:
+    Number of observations              72
+    Fixed effects coefficients           4
+    Random effects coefficients          7
+    Covariance parameters                2
+
+Formula:
+    behavior ~ 1 + day*stim + (1 | rat)
+
+Model fit statistics:
+    AIC       BIC      LogLikelihood    Deviance
+    587.74    601.4    -287.87          575.74  
+
+Fixed effects coefficients (95% CIs):
+    Name                   Estimate    SE         tStat      DF    pValue    
+    {'(Intercept)'}        17.308       5.4187     3.1941    68     0.0021263
+    {'day'        }        8.0028      0.78816     10.154    68    2.9164e-15
+    {'stim'       }        1.8623       7.1264    0.26132    68       0.79463
+    {'day:stim'   }         1.604      0.98592     1.6269    68       0.10838
+
+
+    Lower       Upper 
+      6.4953    28.121
+      6.4301    9.5756
+     -12.358    16.083
+    -0.36337    3.5714
+
+Random effects covariance parameters (95% CIs):
+Group: rat (7 Levels)
+    Name1                  Name2                  Type           Estimate
+    {'(Intercept)'}        {'(Intercept)'}        {'std'}        6.0474  
+
+
+    Lower     Upper 
+    2.8281    12.931
+
+Group: Error
+    Name               Estimate    Lower    Upper 
+    {'Res Std'}        12.425      10.46    14.759
+
+
+effect                     t(df) / F(df1)     p (resid)    Satterthwaite: p (df)
+----------------------------------------------------------------------------
+stim x day (interaction)   t(68)=  1.63 F(1)=  2.647  p=0.1084    p=0.1084 (df=68)
+day (learning)             t(68)= 10.15 F(1)=103.100  p=2.916e-15    p=2.368e-15 (df=69)
+stim (main, window start)  t(68)=  0.26 F(1)=  0.068  p=0.7946    p=0.7968 (df=18)
+interaction 95% CI: [-0.36, +3.57]
+HONEST LME (per-animal random slope, day|rat): interaction F(1,22.2)=2.05, p=0.1666
+  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope
+   model above is the honest learning-rate test -- DF collapses toward the animal count.)
+INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.1084, slope diff=+1.60)
+Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
diff --git a/analysis/matlab/variations/boxa_b2_d0_5/analyze.m b/analysis/matlab/variations/boxa_b2_d0_5/analyze.m
new file mode 100644
index 0000000..82da5ed
--- /dev/null
+++ b/analysis/matlab/variations/boxa_b2_d0_5/analyze.m
@@ -0,0 +1,98 @@
+% Variation analysis -- the paper's linear mixed model on the successful-reach
+% COUNT, fit on this folder's curated data subset.
+%
+%   model: behavior ~ stim + day + stim:day + (1|rat)
+%     behavior = successful reaches (count per session)
+%     stim     = 1 for the treatment group(s), 0 for the control group(s)
+%     day      = training day within this window (0 = first analyzed day)
+%     rat      = subject (random intercept)
+%
+% Self-contained: reads data.csv beside this script and writes result.txt.
+% Run headless from this folder with:   matlab -batch "analyze"
+% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
+
+here = fileparts(mfilename('fullpath'));
+if isempty(here); here = pwd; end
+vname = regexprep(here, '.*[/\\]', '');            % folder name = variation id
+
+D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
+
+tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
+    'VariableNames', {'behavior', 'day', 'stim', 'rat'});
+m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
+C = m.Coefficients; A = anova(m); ci = coefCI(m);
+As = anova(m, 'DFMethod', 'satterthwaite');           % Satterthwaite denominator DF
+
+% Honest test: refit with a per-animal random SLOPE so the interaction DF
+% collapses toward the animal count (guarded -- may not converge in short windows).
+rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
+wst = warning('off', 'all');
+try
+    mr = fitlme(tbl, 'behavior ~ stim + day + stim:day + (day|rat)');
+    Ar = anova(mr, 'DFMethod', 'satterthwaite');
+    ri = strcmp(Ar.Term, 'day:stim');
+    rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
+catch
+end
+warning(wst);
+
+gi = @(t) find(strcmp(C.Name, t), 1);
+ga = @(t) find(strcmp(A.Term, t), 1);
+gs = @(t) find(strcmp(As.Term, t), 1);
+row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f  p=%.4g    p=%.4g (df=%.0f)\n', nm, ...
+    C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
+    As.pValue(gs(t)), As.DF2(gs(t)));
+
+maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
+maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
+if abs(maxT - maxC) > 2
+    cov = '** WARNING: unequal day coverage -- interaction may be confounded. **';
+else
+    cov = '(equal day coverage over this window)';
+end
+
+ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
+if pI >= 0.05
+    verdict = 'n.s. -- slopes parallel (no differential learning rate)';
+elseif eI > 0
+    verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
+else
+    verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
+end
+
+bar = repmat('=', 1, 78);
+raw = regexprep(evalc('disp(m)'), '</?strong>', '');
+s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
+s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)\n')];
+s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
+s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
+s = [s sprintf('control   (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
+s = [s sprintf('N = %d rats, %d sessions   raw day coverage: treat %d..%d, control %d..%d\n', ...
+    numel(unique(D.subject)), height(D), minT, maxT, minC, maxC)];
+s = [s sprintf('%s\n\n%s\nFULL MODEL SUMMARY -- fitlme\n%s\n%s\n', cov, bar, bar, raw)];
+s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (resid)', 'Satterthwaite: p (df)', repmat('-', 1, 76))];
+s = [s row('stim x day (interaction)', 'day:stim')];
+s = [s row('day (learning)', 'day')];
+s = [s row('stim (main, window start)', 'stim')];
+s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
+if rsOk
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
+else
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): did not converge for this window.\n')];
+end
+s = [s sprintf(['  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
+    '   model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
+s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
+s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
+
+fprintf('%s', s);
+fid = fopen(fullfile(here, 'result.txt'), 'w');
+fprintf(fid, '%s', s);
+fclose(fid);
+
+% Machine-readable handoff for the summary table (see make_variations.m).
+VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
+    'nObs', height(D), 'interP', pI, 'interEst', eI, ...
+    'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
+    'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
+    'covEqual', abs(maxT - maxC) <= 2);
diff --git a/analysis/matlab/variations/boxa_b2_d0_5/data.csv b/analysis/matlab/variations/boxa_b2_d0_5/data.csv
new file mode 100644
index 0000000..bb2761a
--- /dev/null
+++ b/analysis/matlab/variations/boxa_b2_d0_5/data.csv
@@ -0,0 +1,43 @@
+subject,group,day,success,total,stim
+Banh-mi-1,Electrode-Box-B2,0,13,84,1
+Banh-mi-1,Electrode-Box-B2,1,35,86,1
+Banh-mi-1,Electrode-Box-B2,2,47,110,1
+Banh-mi-1,Electrode-Box-B2,3,65,140,1
+Banh-mi-1,Electrode-Box-B2,4,70,127,1
+Banh-mi-1,Electrode-Box-B2,5,102,142,1
+Banh-mi-2,Electrode-Box-A2,0,25,97,0
+Banh-mi-2,Electrode-Box-A2,1,22,101,0
+Banh-mi-2,Electrode-Box-A2,2,22,119,0
+Banh-mi-2,Electrode-Box-A2,3,29,118,0
+Banh-mi-2,Electrode-Box-A2,4,27,136,0
+Banh-mi-2,Electrode-Box-A2,5,43,146,0
+Egg-tart-1,Electrode-Box-B2,0,7,56,1
+Egg-tart-1,Electrode-Box-B2,1,16,78,1
+Egg-tart-1,Electrode-Box-B2,2,23,103,1
+Egg-tart-1,Electrode-Box-B2,3,63,120,1
+Egg-tart-1,Electrode-Box-B2,4,69,132,1
+Egg-tart-1,Electrode-Box-B2,5,83,136,1
+Egg-tart-2,Electrode-Box-A2,0,9,32,0
+Egg-tart-2,Electrode-Box-A2,1,2,38,0
+Egg-tart-2,Electrode-Box-A2,2,31,93,0
+Egg-tart-2,Electrode-Box-A2,3,44,101,0
+Egg-tart-2,Electrode-Box-A2,4,54,131,0
+Egg-tart-2,Electrode-Box-A2,5,84,139,0
+Khoai-tay-1,Electrode-Box-A,0,14,62,1
+Khoai-tay-1,Electrode-Box-A,1,22,83,1
+Khoai-tay-1,Electrode-Box-A,2,6,79,1
+Khoai-tay-1,Electrode-Box-A,3,28,97,1
+Khoai-tay-1,Electrode-Box-A,4,60,134,1
+Khoai-tay-1,Electrode-Box-A,5,75,138,1
+Root-beer-1,Electrode-Box-B2,0,11,85,1
+Root-beer-1,Electrode-Box-B2,1,18,76,1
+Root-beer-1,Electrode-Box-B2,2,40,105,1
+Root-beer-1,Electrode-Box-B2,3,55,134,1
+Root-beer-1,Electrode-Box-B2,4,75,136,1
+Root-beer-1,Electrode-Box-B2,5,64,133,1
+Root-beer-2,Electrode-Box-A2,0,22,74,0
+Root-beer-2,Electrode-Box-A2,1,31,87,0
+Root-beer-2,Electrode-Box-A2,2,49,134,0
+Root-beer-2,Electrode-Box-A2,3,31,89,0
+Root-beer-2,Electrode-Box-A2,4,60,140,0
+Root-beer-2,Electrode-Box-A2,5,84,147,0
diff --git a/analysis/matlab/variations/boxa_b2_d0_5/result.txt b/analysis/matlab/variations/boxa_b2_d0_5/result.txt
new file mode 100644
index 0000000..c30cd85
--- /dev/null
+++ b/analysis/matlab/variations/boxa_b2_d0_5/result.txt
@@ -0,0 +1,68 @@
+==============================================================================
+VARIATION: boxa_b2_d0_5
+==============================================================================
+model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)
+day = training day within window (0 = first analyzed day)
+treatment (stim=1): Electrode-Box-A, Electrode-Box-B2
+control   (stim=0): Electrode-Box-A2
+N = 7 rats, 42 sessions   raw day coverage: treat 0..5, control 0..5
+(equal day coverage over this window)
+
+==============================================================================
+FULL MODEL SUMMARY -- fitlme
+==============================================================================
+
+Linear mixed-effects model fit by ML
+
+Model information:
+    Number of observations              42
+    Fixed effects coefficients           4
+    Random effects coefficients          7
+    Covariance parameters                2
+
+Formula:
+    behavior ~ 1 + day*stim + (1 | rat)
+
+Model fit statistics:
+    AIC       BIC       LogLikelihood    Deviance
+    341.62    352.04    -164.81          329.62  
+
+Fixed effects coefficients (95% CIs):
+    Name                   Estimate    SE        tStat       DF    pValue    
+    {'(Intercept)'}         12.524     5.8127      2.1546    38      0.037599
+    {'day'        }         9.8571     1.5596      6.3204    38    2.0711e-07
+    {'stim'       }        -4.7262     7.6895    -0.61463    38       0.54246
+    {'day:stim'   }         4.7071     2.0631      2.2816    38      0.028206
+
+
+    Lower      Upper 
+    0.75663    24.291
+     6.6999    13.014
+    -20.293     10.84
+    0.53057    8.8837
+
+Random effects covariance parameters (95% CIs):
+Group: rat (7 Levels)
+    Name1                  Name2                  Type           Estimate
+    {'(Intercept)'}        {'(Intercept)'}        {'std'}        5.8715  
+
+
+    Lower    Upper 
+    2.486    13.868
+
+Group: Error
+    Name               Estimate    Lower     Upper 
+    {'Res Std'}        11.3        8.9402    14.283
+
+
+effect                     t(df) / F(df1)     p (resid)    Satterthwaite: p (df)
+----------------------------------------------------------------------------
+stim x day (interaction)   t(38)=  2.28 F(1)=  5.206  p=0.02821    p=0.02871 (df=35)
+day (learning)             t(38)=  6.32 F(1)= 39.947  p=2.071e-07    p=2.928e-07 (df=35)
+stim (main, window start)  t(38)= -0.61 F(1)=  0.378  p=0.5425    p=0.5456 (df=20)
+interaction 95% CI: [+0.53, +8.88]
+HONEST LME (per-animal random slope, day|rat): interaction F(1,7.0)=2.92, p=0.131
+  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope
+   model above is the honest learning-rate test -- DF collapses toward the animal count.)
+INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.02821, slope diff=+4.71)
+Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
diff --git a/analysis/matlab/variations/boxa_b2_d6_10/analyze.m b/analysis/matlab/variations/boxa_b2_d6_10/analyze.m
new file mode 100644
index 0000000..82da5ed
--- /dev/null
+++ b/analysis/matlab/variations/boxa_b2_d6_10/analyze.m
@@ -0,0 +1,98 @@
+% Variation analysis -- the paper's linear mixed model on the successful-reach
+% COUNT, fit on this folder's curated data subset.
+%
+%   model: behavior ~ stim + day + stim:day + (1|rat)
+%     behavior = successful reaches (count per session)
+%     stim     = 1 for the treatment group(s), 0 for the control group(s)
+%     day      = training day within this window (0 = first analyzed day)
+%     rat      = subject (random intercept)
+%
+% Self-contained: reads data.csv beside this script and writes result.txt.
+% Run headless from this folder with:   matlab -batch "analyze"
+% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
+
+here = fileparts(mfilename('fullpath'));
+if isempty(here); here = pwd; end
+vname = regexprep(here, '.*[/\\]', '');            % folder name = variation id
+
+D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
+
+tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
+    'VariableNames', {'behavior', 'day', 'stim', 'rat'});
+m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
+C = m.Coefficients; A = anova(m); ci = coefCI(m);
+As = anova(m, 'DFMethod', 'satterthwaite');           % Satterthwaite denominator DF
+
+% Honest test: refit with a per-animal random SLOPE so the interaction DF
+% collapses toward the animal count (guarded -- may not converge in short windows).
+rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
+wst = warning('off', 'all');
+try
+    mr = fitlme(tbl, 'behavior ~ stim + day + stim:day + (day|rat)');
+    Ar = anova(mr, 'DFMethod', 'satterthwaite');
+    ri = strcmp(Ar.Term, 'day:stim');
+    rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
+catch
+end
+warning(wst);
+
+gi = @(t) find(strcmp(C.Name, t), 1);
+ga = @(t) find(strcmp(A.Term, t), 1);
+gs = @(t) find(strcmp(As.Term, t), 1);
+row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f  p=%.4g    p=%.4g (df=%.0f)\n', nm, ...
+    C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
+    As.pValue(gs(t)), As.DF2(gs(t)));
+
+maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
+maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
+if abs(maxT - maxC) > 2
+    cov = '** WARNING: unequal day coverage -- interaction may be confounded. **';
+else
+    cov = '(equal day coverage over this window)';
+end
+
+ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
+if pI >= 0.05
+    verdict = 'n.s. -- slopes parallel (no differential learning rate)';
+elseif eI > 0
+    verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
+else
+    verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
+end
+
+bar = repmat('=', 1, 78);
+raw = regexprep(evalc('disp(m)'), '</?strong>', '');
+s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
+s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)\n')];
+s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
+s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
+s = [s sprintf('control   (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
+s = [s sprintf('N = %d rats, %d sessions   raw day coverage: treat %d..%d, control %d..%d\n', ...
+    numel(unique(D.subject)), height(D), minT, maxT, minC, maxC)];
+s = [s sprintf('%s\n\n%s\nFULL MODEL SUMMARY -- fitlme\n%s\n%s\n', cov, bar, bar, raw)];
+s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (resid)', 'Satterthwaite: p (df)', repmat('-', 1, 76))];
+s = [s row('stim x day (interaction)', 'day:stim')];
+s = [s row('day (learning)', 'day')];
+s = [s row('stim (main, window start)', 'stim')];
+s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
+if rsOk
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
+else
+    s = [s sprintf('HONEST LME (per-animal random slope, day|rat): did not converge for this window.\n')];
+end
+s = [s sprintf(['  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
+    '   model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
+s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
+s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
+
+fprintf('%s', s);
+fid = fopen(fullfile(here, 'result.txt'), 'w');
+fprintf(fid, '%s', s);
+fclose(fid);
+
+% Machine-readable handoff for the summary table (see make_variations.m).
+VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
+    'nObs', height(D), 'interP', pI, 'interEst', eI, ...
+    'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
+    'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
+    'covEqual', abs(maxT - maxC) <= 2);
diff --git a/analysis/matlab/variations/boxa_b2_d6_10/data.csv b/analysis/matlab/variations/boxa_b2_d6_10/data.csv
new file mode 100644
index 0000000..238de8e
--- /dev/null
+++ b/analysis/matlab/variations/boxa_b2_d6_10/data.csv
@@ -0,0 +1,31 @@
+subject,group,day,success,total,stim
+Banh-mi-1,Electrode-Box-B2,6,90,131,1
+Banh-mi-1,Electrode-Box-B2,7,109,148,1
+Banh-mi-1,Electrode-Box-B2,8,104,137,1
+Banh-mi-1,Electrode-Box-B2,9,119,150,1
+Banh-mi-1,Electrode-Box-B2,10,121,158,1
+Banh-mi-2,Electrode-Box-A2,6,74,146,0
+Banh-mi-2,Electrode-Box-A2,7,70,148,0
+Banh-mi-2,Electrode-Box-A2,8,65,130,0
+Banh-mi-2,Electrode-Box-A2,9,79,151,0
+Banh-mi-2,Electrode-Box-A2,10,93,152,0
+Egg-tart-1,Electrode-Box-B2,6,71,142,1
+Egg-tart-1,Electrode-Box-B2,7,79,138,1
+Egg-tart-1,Electrode-Box-B2,8,98,142,1
+Egg-tart-1,Electrode-Box-B2,9,89,139,1
+Egg-tart-1,Electrode-Box-B2,10,96,143,1
+Egg-tart-2,Electrode-Box-A2,6,85,145,0
+Egg-tart-2,Electrode-Box-A2,7,79,143,0
+Egg-tart-2,Electrode-Box-A2,8,76,131,0
+Egg-tart-2,Electrode-Box-A2,9,88,149,0
+Egg-tart-2,Electrode-Box-A2,10,81,151,0
+Khoai-tay-1,Electrode-Box-A,6,81,137,1
+Khoai-tay-1,Electrode-Box-A,7,78,147,1
+Khoai-tay-1,Electrode-Box-A,8,91,132,1
+Khoai-tay-1,Electrode-Box-A,9,99,146,1
+Khoai-tay-1,Electrode-Box-A,10,94,143,1
+Root-beer-1,Electrode-Box-B2,6,104,139,1
+Root-beer-1,Electrode-Box-B2,7,98,148,1
+Root-beer-1,Electrode-Box-B2,8,81,145,1
+Root-beer-1,Electrode-Box-B2,9,89,156,1
+Root-beer-1,Electrode-Box-B2,10,105,158,1
diff --git a/analysis/matlab/variations/boxa_b2_d6_10/result.txt b/analysis/matlab/variations/boxa_b2_d6_10/result.txt
new file mode 100644
index 0000000..4d49b78
--- /dev/null
+++ b/analysis/matlab/variations/boxa_b2_d6_10/result.txt
@@ -0,0 +1,68 @@
+==============================================================================
+VARIATION: boxa_b2_d6_10
+==============================================================================
+model: behavior ~ stim + day + stim:day + (1|rat)   (behavior = success COUNT)
+day = training day within window (0 = first analyzed day)
+treatment (stim=1): Electrode-Box-A, Electrode-Box-B2
+control   (stim=0): Electrode-Box-A2
+N = 6 rats, 30 sessions   raw day coverage: treat 6..10, control 6..10
+(equal day coverage over this window)
+
+==============================================================================
+FULL MODEL SUMMARY -- fitlme
+==============================================================================
+
+Linear mixed-effects model fit by ML
+
+Model information:
+    Number of observations              30
+    Fixed effects coefficients           4
+    Random effects coefficients          6
+    Covariance parameters                2
+
+Formula:
+    behavior ~ 1 + day*stim + (1 | rat)
+
+Model fit statistics:
+    AIC       BIC       LogLikelihood    Deviance
+    231.05    239.46    -109.52          219.05  
+
+Fixed effects coefficients (95% CIs):
+    Name                   Estimate    SE        tStat      DF    pValue    
+    {'(Intercept)'}        74.2        6.2618      11.85    26    5.5401e-12
+    {'day'        }         2.4        1.8164     1.3213    26        0.1979
+    {'stim'       }          12        7.6691     1.5647    26       0.12974
+    {'day:stim'   }         1.9        2.2246    0.85409    26       0.40085
+
+
+    Lower      Upper 
+     61.329    87.071
+    -1.3336    6.1336
+     -3.764    27.764
+    -2.6727    6.4727
+
+Random effects covariance parameters (95% CIs):
+Group: rat (6 Levels)
+    Name1                  Name2                  Type           Estimate
+    {'(Intercept)'}        {'(Intercept)'}        {'std'}        6.2314  
+
+
+    Lower     Upper
+    2.9021    13.38
+
+Group: Error
+    Name               Estimate    Lower     Upper 
+    {'Res Std'}        8.123       6.1215    10.779
+
+
+effect                     t(df) / F(df1)     p (resid)    Satterthwaite: p (df)
+----------------------------------------------------------------------------
+stim x day (interaction)   t(26)=  0.85 F(1)=  0.729  p=0.4009    p=0.4015 (df=24)
+day (learning)             t(26)=  1.32 F(1)=  1.746  p=0.1979    p=0.1989 (df=24)
+stim (main, window start)  t(26)=  1.56 F(1)=  2.448  p=0.1297    p=0.142 (df=13)
+interaction 95% CI: [-2.67, +6.47]
+HONEST LME (per-animal random slope, day|rat): interaction F(1,6.0)=0.61, p=0.4631
+  (Satterthwaite DF ~= residual on this random-intercept model; the random-slope
+   model above is the honest learning-rate test -- DF collapses toward the animal count.)
+INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.4009, slope diff=+1.90)
+Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
diff --git a/analysis/matlab/variations/boxa_summary.csv b/analysis/matlab/variations/boxa_summary.csv
new file mode 100644
index 0000000..a1ff91c
--- /dev/null
+++ b/analysis/matlab/variations/boxa_summary.csv
@@ -0,0 +1,7 @@
+variation,grouping,window,nRats,nObs,interaction_p,interaction_p_satt,interaction_p_rs,interaction_est,stim_p
+boxa_a2_d0_5,boxa_a2,d0_5,7,42,0.0394100933497468,0.0387789570155384,0.134104933837531,4.67619047619049,0.934195418405735
+boxa_a2_d6_10,boxa_a2,d6_10,6,30,0.640948346908124,0.641274155570031,0.687193776070756,0.999999999999999,0.0959886716436125
+boxa_a2_d0_10,boxa_a2,d0_10,7,72,0.241397580095809,0.241169749082988,0.42454522177293,1.1984817751552,0.251749522863234
+boxa_b2_d0_5,boxa_b2,d0_5,7,42,0.0282058457878955,0.0287096264064777,0.131045386717085,4.70714285714286,0.542460948914174
+boxa_b2_d6_10,boxa_b2,d6_10,6,30,0.400853858752755,0.401497745167433,0.463104837605515,1.9,0.129739412462252
+boxa_b2_d0_10,boxa_b2,d0_10,7,72,0.108380977202764,0.108390578891594,0.166606323000689,1.60400880375615,0.794632800487035