Signal combination in flutter vibration perception.

Code ↔ Paper

6 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 6 matches

[1] § Results › Summation and suppression effects on flutter thresholds ↔ TestVersion.qmd, lines 1527–1646 · score 1.00 · approximately parallel handles, Greenhouse Geisser corrected, probability summation instead, 0.5–2 %, baseline stimuli vibrated, cumulative Gaussians
[2] § Results › Summation and suppression of neural responses ↔ TestVersion.qmd, lines 1527–1646 · score 0.99 · perfect linear summation, Greenhouse Geisser corrected, sub linear summation, EEG amplitudes increased, high baseline response, Responses increased monotonically
[3] § Results › Summation and suppression of neural responses ↔ vibrosummanuscript.qmd, lines 1470–1552 · score 0.98 · Greenhouse Geisser corrected, sub linear summation, EEG amplitudes increased, high baseline response, 32–64 %, increased monotonically
[4] § Results › Computational modelling results ↔ vibrosummanuscript.qmd, lines 1470–1552 · score 0.98 · directly tap mechanisms, sub linear summation, excellent account, suppression evident, neural population, sensitive subset
[5] § Materials and methods › Psychophysical procedures ↔ Experiment code/tactiledippers.m, lines 93–174 · score 0.58 · foot pedal, beep, pressing, headphones, staircases, interval
[6] § Materials and methods › EEG procedures ↔ Experiment code/tactileEEG.m, lines 1–60 · score 0.51 · repeated twice, blocks, SSSEP, EEG

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Quarto · 1,734 lines · 109 KB · no license · 2 matches

---
title: "Signal combination in vibration perception"
author:
- name: Shasha Wei$^{1,2}$
- name: Alex R. Wade$^{1,3}$
- name: Catherine E.J. Preston$^{1}$
- name: \& Daniel H. Baker$^{1}$
format: pdf
# prefer-html: true # required if outputting to docx format
bibliography: references.bib
csl: pnas.csl
execute:
echo: false
output: false
---
$^1$Department of Psychology, University of York, UK, YO10 5DD\
$^2$Corresponding author, email: mfv507\@york.ac.uk\
$^3$York Biomedical Research Institute, University of York, UK, YO10 5DD\
**ORCID**:\
Shasha Wei: [https://orcid.org/0009-0002-8975-4214](https://orcid.org/0009-0002-8975-4214)\
Alex Wade: [https://orcid.org/0000-0003-4871-2747](https://orcid.org/0000-0003-4871-2747)\
Catherine Preston: [https://orcid.org/0000-0001-7158-5382](https://orcid.org/0000-0001-7158-5382)\
Daniel Baker: [https://orcid.org/0000-0002-0161-443X](https://orcid.org/0000-0002-0161-443X)\
**Corresponding author**\
Shasha Wei, Email: mfv507\@york.ac.uk
**Present address**:\
$^1$Department of Psychology, University of York, York, UK, YO10 5DD
**Author Contributions**\
D.H.B. designed research; S.W. performed research; S.W., A.R.W. and D.H.B. analysed data; and S.W., A.R.W., C.E.J.P. and D.H.B. wrote and revised the paper.
**Competing Interests**\
The authors declare no competing interest.
**Classification**: Social Science-Psychological and Cognitive Sciences\
**Keywords**: vibrotactile summation, suppression, somatosensory, computational modelling\
This manuscript was deposited as a preprint on PsyArXiv (DOI: 10.31234/osf.io/yjv6s)
under the CC-BY Attribution 4.0 International license.
```{r initialiseenvironment}
#| include: false
doanalysis <- 0
fitmodels <- 0
plotfigures <- 1
# install R packages
packagelist <- c('utils','knitr','reticulate','osfr','tinytex','rstatix','kableExtra')
missingpackages <- packagelist[!packagelist %in% installed.packages()[,1]]
if (length(missingpackages)>0){install.packages(missingpackages)}
toinstall <- packagelist[which(!packagelist %in% (.packages()))]
invisible(lapply(toinstall,library,character.only=TRUE))
use_python("/usr/bin/python3.10")
```
```{python setup}
#| include: false
# import packages etc.
import os
import mne
import numpy as np
import pandas as pd
import psignifit as ps
from matplotlib.lines import Line2D
import matplotlib.pyplot as plt
from mne import EpochsArray
from mne import viz
from mne.channels import make_standard_montage
import scipy.stats as stats
from scipy.fft import fft, fftfreq
from scipy.stats import norm
from scipy.optimize import minimize
from matplotlib.lines import Line2D
from matplotlib import cm, colors, colorbar
from statsmodels.stats.anova import AnovaRM
import pingouin as pg
from pingouin import pairwise_ttests
from mpl_toolkits.axes_grid1 import make_axes_locatable
```
```{r checkfordata}
#| include: false
# check if raw data are available and download if required
if (!dir.exists('local/')){dir.create('local/')}
osfproject <- osf_retrieve_node('m79d2')
osffiles <- osf_ls_files(osfproject,n_max=300)
if (!file.exists('local/8subjects.csv')){osf_download(osffiles[which(osffiles$name=='8subjects.csv'),],'local/',progress=TRUE)}
if (!file.exists('local/allthresh.npz')){osf_download(osffiles[which(osffiles$name=='allthresh.npz'),],'local/',progress=TRUE)}
if (!file.exists('local/allslope.npz')){osf_download(osffiles[which(osffiles$name=='allslope.npz'),],'local/',progress=TRUE)}
if (!file.exists('local/allF1.npy')){osf_download(osffiles[which(osffiles$name=='allF1.npy'),],'local/',progress=TRUE)}
if (!file.exists('local/allF2.npy')){osf_download(osffiles[which(osffiles$name=='allF2.npy'),],'local/',progress=TRUE)}
if (!file.exists('local/dippermodels.npz')){osf_download(osffiles[which(osffiles$name=='dippermodels.npz'),],'local/',progress=TRUE)}
if (!file.exists('local/EEGmodels.npz')){osf_download(osffiles[which(osffiles$name=='EEGmodels.npz'),],'local/',progress=TRUE)}
if (!file.exists('local/EEG_allresp.npy')){osf_download(osffiles[which(osffiles$name=='EEG_allresp.npy'),],'local/',progress=TRUE)}
if ((plotfigures+fitmodels)>0){
if (!file.exists('local/processedEEGdata.npz')){osf_download(osffiles[which(osffiles$name=='processedEEGdata.npz'),],'local/',progress=TRUE)}
if (!file.exists('local/montage-info.fif')){osf_download(osffiles[which(osffiles$name=='montage-info.fif'),],'local/',progress=TRUE)}
}
if (doanalysis==1){
osfproject <- osf_retrieve_node('z2fcm')
osffiles <- osf_ls_files(osfproject,n_max=300)
if (!dir.exists('local/rawEEGdata/')){dir.create('local/rawEEGdata/')
osf_download(osffiles[which(osffiles$name=='setfiles.zip'),],'local/',progress=TRUE,conflicts='skip')
unzip('local/setfiles.zip',exdir='local/rawEEGdata/')
file.remove('local/setfiles.zip')
d <- dir('local/rawEEGdata/setfiles/',full.names=TRUE)
file.copy(d,'local/rawEEGdata/')
file.remove(d)
for (i in 1:4){
osf_download(osffiles[which(osffiles$name==paste0('fdt',i,'.zip')),],'local/',progress=TRUE,conflicts='skip')
unzip(paste0('local/fdt',i,'.zip'),exdir='local/rawEEGdata/')
file.remove(paste0('local/fdt',i,'.zip'))
d <- dir(paste0('local/rawEEGdata/fdt',i),full.names=TRUE)
file.copy(d,'local/rawEEGdata/')
file.remove(d)
}
}
# for (s in 1:31){
# if (!file.exists(paste0('local/rawEEGdata/P',s,'_EEG.fdt'))){osf_download(osffiles[which(osffiles$name==paste0('P',s,'_EEG.fdt')),],'local/rawEEGdata/',progress=TRUE,conflicts='skip')}
# if (!file.exists(paste0('local/rawEEGdata/P',s,'_EEG.set'))){osf_download(osffiles[which(osffiles$name==paste0('P',s,'_EEG.set')),],'local/rawEEGdata/',progress=TRUE,conflicts='skip')}
# }
}
```
```{python analysedippers}
#| include: false
if r.doanalysis:
row = pd.read_csv('local/8subjects.csv')
ex1data = row[['Subject','Condition', 'PedestalContrast','TargetContrast','IsCorrect']].values
df = pd.DataFrame(ex1data, columns=['Subject','Condition', 'PedestalContrast','TargetContrast', 'IsCorrect'])
df['Condition'] = (df['Condition'] + 1) // 2
sublist = df['Subject'].unique()
#%% generate tofit data and fit psychometric functions
allthresh = np.zeros((len(sublist), 4, 8))
allslope = np.zeros((len(sublist), 4, 8))
# Loop over each subject
for i, subject in enumerate(sublist):
subdata = df[df['Subject'] == subject]
# Loop over conditions (1 to 4)
for cond in range(1, 5):
conddata = subdata[subdata['Condition'] == cond]
pedlevs = np.sort(conddata['PedestalContrast'].unique())
# Loop over pedestal levels (1 to 8)
for pedlevel in range(8):
if pedlevel < len(pedlevs):
blockdata = conddata[conddata['PedestalContrast'] == pedlevs[pedlevel]]
# Proceed only if there's data in blockdata
if not blockdata.empty:
targetcontrasts = np.sort(blockdata['TargetContrast'].unique())
targetcontrasts = pd.to_numeric(targetcontrasts, errors='coerce')
#targetcontrasts = targetcontrasts[targetcontrasts > 0]
ntrials = []
ncorrect = []
# Loop over each TargetContrast
for target in targetcontrasts:
temp = blockdata[blockdata['TargetContrast'] == target]
# Count ncorr and ntotal
ntrials.append(len(temp))
ncorrect.append(temp['IsCorrect'].sum())
level = np.round(20 * np.log10(targetcontrasts))
tofit = np.vstack((level, ncorrect, ntrials)).T
result_fit = ps.psignifit(tofit, experiment_type='2AFC')
# result_params = result_fit.get_parameters_estimate()
allthresh[i, cond - 1, pedlevel] = result_fit.parameter_estimate['threshold']
allslope[i, cond - 1, pedlevel] = 10.3 / (result_fit.parameter_estimate['width'] / (norm.ppf(1 - 0.05) - norm.ppf(0.05)))
#print("Thresholds:", allthresh)
#print("Slopes:", allslope)
allthresh[:, 2, 1:8] = allthresh[:, 2, 0:7]
allthresh[:, 3, 1:8] = allthresh[:, 3, 0:7]
allslope[:, 2, 1:8] = allslope[:, 2, 0:7]
allslope[:, 3, 1:8] = allslope[:, 3, 0:7]
allthresh[:, 2, 0] = allthresh[:, 0,0]
allthresh[:, 3, 0] = allthresh[:, 0, 0]
allslope[:, 2, 0] = allslope[:, 0,0]
allslope[:, 3, 0] = allslope[:, 0, 0]
np.savez('local/allthresh.npz', allthresh=allthresh)
np.savez('local/allslope.npz', allslope=allslope)
```
```{python plotdippers}
#| include: false
if r.plotfigures:
data = np.load('local/allthresh.npz')
allthresh = data['allthresh']
meanthresh=np.mean(allthresh, axis=0) #4*8
SEthresh = np.std(allthresh, axis=0, ddof=1) / np.sqrt(8)
x = np.arange(1, 9)
markers=['o','s','D','^']
lineColor = ['b', 'r', 'orange','g'] # Line colors for the four conditions
labels= ['Pentadactyl', 'Dekadactyl', 'Half-Dekadactyl','Dichodactyl']
markeredgecolors=['black','black','black','black']
legend_elements = []
#%% Plotting
plt.figure(figsize=(20, 10))
plt.subplot(1, 2, 1)
for m in range(meanthresh.shape[0]):#iterate over the columns of average_array
y=meanthresh[m, :]
y_error=SEthresh[m, :] #because SE array have the same number of columns as average_array
y_upper = y + y_error
y_lower = y - y_error
plt.plot(x, y, label=labels[m], color=lineColor[m], marker=markers[m], markersize=10, markeredgecolor=markeredgecolors[m])
plt.fill_between(x, y_lower, y_upper,alpha=0.15, color=lineColor[m])
legend_elements.append(Line2D([0], [0], color=lineColor[m], linestyle='-', label=labels[m], markersize=10, marker=markers[m],markeredgecolor=markeredgecolors[m]))
plt.ylim(-15,30)
plt.yticks([-12,-6,0,6,12,18,24,30],[0.25, 0.5, 1, 2, 4, 8, 16, 32],fontsize = 25)
plt.xticks(x,[0, 0.5, 1, 2, 4, 8, 16, 32], fontsize = 25 )
plt.xlabel('Baseline intensity level (%)', fontsize=30)
plt.ylabel('Threshold (%)', fontsize=30)
plt.legend(handles=legend_elements)
leg = plt.legend(frameon=False, loc='lower right', fontsize=24)
for text in leg.get_texts():
text.set_fontsize(18)
fig = plt.gcf()
plt.text(0.085, 0.92, '(a)', fontsize=40, va='top', ha='left', transform=fig.transFigure)
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(2.5)
ax.tick_params(axis='both', which='major', width=2)
################ plot slope
data1= np.load('local/allslope.npz')
allslope = data1['allslope']
allslope = 20 * np.log10(allslope)
meanslope = np.mean(allslope, axis=0) #4*8
SEslope = np.std(allslope, axis=0, ddof=1) / np.sqrt(8)
plt.subplot(1, 2, 2)
for m in range(meanslope.shape[0]):#iterate over the columns of average_array
y=meanslope[m, :]
y_error=SEslope[m,:] #because SE array have the same number of columns as average_array
y_upper = y + y_error
y_lower = y - y_error
plt.plot(x, y, label=labels[m], color=lineColor[m], marker=markers[m], markersize=10,markeredgecolor=markeredgecolors[m])
plt.fill_between(x, y_lower, y_upper,alpha=0.15, color=lineColor[m])
plt.ylim(-7,19)
plt.yticks([-6, 0, 6, 12, 18], ['0.5','1', '2', '4', '8'],fontsize=25)
plt.minorticks_off()
plt.xticks(x, [0, 0.5, 1, 2, 4, 8, 16, 32], fontsize=25)
plt.xlabel('Baseline intensity level (%)', fontsize=30)
plt.ylabel('Weibull ' + r"$\mathrm{\beta}$", fontsize=30)
plt.subplots_adjust(wspace=0.6)
fig = plt.gcf()
plt.text(0.58, 0.92, '(b)', fontsize=40, va='top', ha='left', transform=fig.transFigure)
ax = plt.gca()
for spine in ax.spines.values():
spine.set_linewidth(2.5) # Adjust the thickness of the frame
ax.tick_params(axis='both', which='major', width=2)
plt.tight_layout()
plt.savefig('Figures/Figure1.pdf')
```
```{python analyseEEG}
#| include: false
if r.doanalysis:
raw_path = "local/rawEEGdata/" # path containing raw data
nsubjs = 31
electrodenames = np.array(['F1', 'F2', 'Fz', 'FC1', 'FC2', 'FCz']) # electrodes of interest
allsubjsF1 = np.zeros((nsubjs,61,6,5), dtype='complex')
allsubjsF2 = np.zeros((nsubjs,61,6,5), dtype='complex')
selectSpec = np.zeros((nsubjs,len(electrodenames),6,5,10001), dtype='complex')
F1 = 26
F2 = 23
F1index = F1*10
F2index = F2*10
tmin = 1.0
tmax = 11.0
fmin = 1.0
fmax = 50.0
reject_criteria = None #dict(eeg=200e-6)
legaltriggers = np.concatenate([range(11,16),range(21,26),range(31,36),range(41,46),range(51,56),range(61,66),range(71,76),range(81,86),range(91,96),range(101,106),range(111,116),range(121,126)])
eventList = ['Mon/L/C1', 'Mon/L/C2', 'Mon/L/C3', 'Mon/L/C4', 'Mon/L/C5',
'Mon/R/C1', 'Mon/R/C2', 'Mon/R/C3', 'Mon/R/C4', 'Mon/R/C5',
'Bin/LR/C1', 'Bin/LR/C2', 'Bin/LR/C3', 'Bin/LR/C4', 'Bin/LR/C5',
'Bin/RL/C1', 'Bin/RL/C2', 'Bin/RL/C3', 'Bin/RL/C4', 'Bin/RL/C5',
'Dich/L/C1', 'Dich/L/C2', 'Dich/L/C3', 'Dich/L/C4', 'Dich/L/C5',
'Dich/R/C1', 'Dich/R/C2', 'Dich/R/C3', 'Dich/R/C4', 'Dich/R/C5',
'XMon/L/C1', 'XMon/L/C2', 'XMon/L/C3', 'XMon/L/C4', 'XMon/L/C5',
'XMon/R/C1', 'XMon/R/C2', 'XMon/R/C3', 'XMon/R/C4', 'XMon/R/C5',
'XBin/LR/C1', 'XBin/LR/C2', 'XBin/LR/C3', 'XBin/LR/C4', 'XBin/LR/C5',
'XBin/RL/C1', 'XBin/RL/C2', 'XBin/RL/C3', 'XBin/RL/C4', 'XBin/RL/C5',
'XDich/L/C1', 'XDich/L/C2', 'XDich/L/C3', 'XDich/L/C4', 'XDich/L/C5',
'XDich/R/C1', 'XDich/R/C2', 'XDich/R/C3', 'XDich/R/C4', 'XDich/R/C5']
levellist = ["C1","C2","C3","C4","C5"]
condlist = ['Mon','Bin','Dich','XMon','XBin','XDich']
# loop through each individual participant, calculate the spectra and amplitude for 6 conditions and 5 intensity levels
for p in range(0,nsubjs):
filename = raw_path + "P" + str(p+1) + '_EEG.set'
raw = mne.io.read_raw_eeglab(filename,preload=True)
raw.drop_channels(["HEOG", "VEOG", "M1", "M2","Fpz"])
ANT_montage = mne.channels.make_standard_montage("standard_1020")
raw.set_montage(ANT_montage)
info = raw.info
events, event_id = mne.events_from_annotations(raw)
ch_names = np.array(raw.ch_names)
selectE = np.where(np.isin(ch_names,electrodenames))[0] # find indices of electrodes of interest
extracted_values = []
for trigger in legaltriggers:
trigger_str = str(trigger) # convert trigger to string for comparison
if trigger_str in event_id:
extracted_values.append(event_id[trigger_str])
event_dict = dict(zip(eventList, extracted_values))
allepochs = mne.Epochs(raw, events, event_id=event_dict, tmin=tmin, tmax=tmax, baseline=(tmin,tmax), reject=reject_criteria, preload=True)
allblocksF1 = np.zeros((61,len(condlist),len(levellist)), dtype='complex')
allblocksF2 = np.zeros((61,len(condlist),len(levellist)), dtype='complex')
spectra = np.zeros((len(selectE), len(condlist), len(levellist), 10001))
for c in range(len(condlist)):
for l in range(len(levellist)):
condstr = str(condlist[c]) + '/' + str(levellist[l])
temp = allepochs[condstr]
conditionmean = temp.average()
d = 1000000*conditionmean.data # rescale to microvolts (from volts)
s = d.shape
for electrode in range(s[0]):
spec = fft(d[electrode,:])/10001
# Check if the electrode is in the selected electrodes list
if electrode in selectE:
idx = np.where(np.isin(selectE,electrode))[0]
spectra[idx, c, l, :] = spec # store all frequencies from 6 electrodes
allblocksF1[electrode,c,l] = spec[F1index]
allblocksF2[electrode,c,l] = spec[F2index]
selectSpec[p,:,:,:,:] = spectra # store spectra: participant X electrode X condition X level X frequency
allsubjsF1[p,:,:,:] = allblocksF1 ##amplitudes at 26Hz: participant X electrode X condition X level
allsubjsF2[p,:,:,:] = allblocksF2 ##amplitude form 23Hz: participant X electrode X condition X level
np.savez('local/processedEEGdata.npz', allsubjsF1=allsubjsF1, allsubjsF2=allsubjsF2, selectSpec=selectSpec, selectE=selectE, info=info)
info.save('local/montage-info.fif') # save montage info in correct format
```
```{python plotEEGdata}
#| include: false
if r.plotfigures:
excludelist = 13
mask = np.ones(31, dtype=bool)
mask[excludelist] = False
freqs = np.linspace(0,49,491)
x = [12,18,24,30,36]
eegdata = np.load('local/processedEEGdata.npz',allow_pickle=True)
allsubjsF1 = eegdata['allsubjsF1']
allsubjsF2 = eegdata['allsubjsF2']
selectSpec = eegdata['selectSpec']
selectE = eegdata['selectE']
info = mne.io.read_info('local/montage-info.fif')
levellist = ["C1","C2","C3","C4","C5"]
condlist = ['Mon','Bin','Dich','XMon','XBin','XDich']
# Row data for plotting spectra
Deka_spec = np.mean(np.abs(selectSpec[mask, :, 1, 4,:]), axis=1) #### average across 6 electrodes
Penta_spec = np.mean(np.abs(selectSpec[mask, :, 0, 4,:]), axis=1)
CrossP_spec = np.mean(np.abs(selectSpec[mask, :, 3, 4,:]), axis=1)
CrossDeka_spec = np.mean(np.abs(selectSpec[mask, :, 4, 4,:]), axis=1)
#### Do outlier rejection, calculate the data for plotting topmap
topomapF1 = np.zeros((61,len(condlist),len(levellist)))
topomapF2 = np.zeros((61,len(condlist),len(levellist)))
for el in range(61):
for cond in range(6):
for lev in range(5):
temp = np.abs(allsubjsF1[mask,el,cond,lev])
topomapF1[el,cond,lev] = np.mean(temp)
temp = np.abs(allsubjsF2[mask,el,cond,lev])
topomapF2[el,cond,lev] = np.mean(temp)
##### Remove one outlier, the rest of 30 participants are used to calculate mean amplitude across conditions, levels
subjsF1 = np.zeros((sum(mask), 61, len(condlist), len(levellist)))
subjsF2 = np.zeros((sum(mask), 61, len(condlist), len(levellist)))
for el in range(61):
for cond in range(6):
for lev in range(5):
# F1 data
subjsF1[:, el, cond, lev] = np.abs(allsubjsF1[mask, el, cond, lev])
# F2 data
subjsF2[:, el, cond, lev] = np.abs(allsubjsF2[mask, el, cond, lev])
######################## filtered F1, F2 contrast response data
EsubjsF1 = subjsF1[:,selectE, :, :]
dataF1 = np.nanmean(EsubjsF1, axis=1) #average across 6 electrodes
meanF1 = np.mean(dataF1,axis=0)
SEF1 = np.std(dataF1,axis=0)/np.sqrt(sum(mask))
EsubjsF2 = subjsF2[:,selectE, :, :]
dataF2 = np.nanmean(EsubjsF2, axis=1) #average across 6 electrodes
meanF2 = np.mean(dataF2,axis=0)
SEF2 = np.std(dataF2,axis=0)/np.sqrt(sum(mask))
np.save('local/allF1.npy', dataF1)
np.save('local/allF2.npy', dataF2)
####################################### plot spectrum + topomaps
fig=plt.figure(figsize=(30,45))
ax1=fig.add_subplot(3,2,1)
ax2=fig.add_subplot(3,2,2)
ax3=fig.add_subplot(3,2,3)
ax4=fig.add_subplot(3,2,4)
ax5=fig.add_subplot(3,2,5)
ax6=fig.add_subplot(3,2,6)
pos = ax1.get_position() # set the position of each subplot manually
new_pos = [pos.x0, pos.y0+0.05, pos.width, pos.height] ###adjust position
ax1.set_position(new_pos)
pos = ax2.get_position() # set the position of each subplot manually
new_pos = [pos.x0 + 0.05, pos.y0+0.05, pos.width, pos.height] ###adjust position
ax2.set_position(new_pos)
pos = ax3.get_position()
new_pos = [pos.x0, pos.y0, pos.width, pos.height]
ax3.set_position(new_pos)
pos = ax4.get_position()
new_pos = [pos.x0 + 0.05, pos.y0, pos.width, pos.height]
ax4.set_position(new_pos)
pos = ax5.get_position()
new_pos = [pos.x0, pos.y0-0.05, pos.width, pos.height]
ax5.set_position(new_pos)
pos = ax6.get_position()
new_pos = [pos.x0 + 0.05, pos.y0-0.05, pos.width, pos.height]
ax6.set_position(new_pos)
axes = [ax1, ax2, ax3, ax4, ax5, ax6]
line_weight = 2
for ax in axes:
for spine in ax.spines.values():
spine.set_linewidth(line_weight)
freq_range = range(11, 471)
def plot_PSD_with_topomap(spectrum, label, topomap_data, info, ax):
spectrum = spectrum[:,list(freq_range)]
nsubjs = spectrum.shape[0]
psds_mean = spectrum.mean(axis=0)
psds_std = spectrum.std(axis=0)
psds_sr = psds_std / np.sqrt(nsubjs)
ax.plot(freqs[list(freq_range)], psds_mean, label=label)
ax.fill_between(freqs[list(freq_range)], psds_mean - psds_sr, psds_mean + psds_sr, color="b", linewidth=0, alpha=0.2)
ax.set_title(f'{label}', fontsize=35)
ax.set_ylabel('Amplitude (μV)', fontsize=35)
ax.set_xlabel('Frequency (Hz)', fontsize=35)
ax.set_ylim(0, 0.5)
ax.tick_params(axis='both', which='major', labelsize=25)
# Plot topomap
topomap_ax = ax.inset_axes([0.55, 0.55, 0.38, 0.38])
img, _ = mne.viz.plot_topomap(topomap_data, info, axes=topomap_ax, show=False, cmap='viridis', vlim=(0, 0.5))
divider = make_axes_locatable(topomap_ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbar = plt.colorbar(img, cax=cax, ticks=[0, 0.25, 0.5])
cbar.ax.tick_params(labelsize=20)
cbar.ax.yaxis.set_label_position('left')
cbar.ax.yaxis.set_ticks_position('right')
# add the unit at the top center
cbar.ax.set_title('µV', fontsize=25, pad=2, loc='center')
# Create colorbar for topomap
####Add text (a), (b)......to the plot
ax1.text(0.03, 0.87+0.05, '(a)', fontsize=80, transform=fig.transFigure)
ax2.text(0.51, 0.87+0.05, '(b)', fontsize=80, transform=fig.transFigure)
ax3.text(0.03, 0.55+0.05, '(c)', fontsize=80, transform=fig.transFigure)
ax4.text(0.51, 0.55+0.05, '(d)', fontsize=80, transform=fig.transFigure)
ax5.text(0.03, 0.23+0.05, '(e)', fontsize=80, transform=fig.transFigure)
ax6.text(0.51, 0.23+0.05, '(f)', fontsize=80, transform=fig.transFigure)
plot_PSD_with_topomap(Penta_spec, 'Pentadactyl', topomapF1[:61, 0, 4], info, ax1)
plot_PSD_with_topomap(Deka_spec, 'Dekadactyl', topomapF1[:61, 1, 4], info, ax2)
plot_PSD_with_topomap(CrossP_spec, 'Cross-Pentadactyl', topomapF2[:61, 3, 4], info, ax3)
plot_PSD_with_topomap(CrossDeka_spec, 'Cross-Dekadactyl', topomapF1[:61, 4, 4], info, ax4)
################################################## Plot contrast response function at 26Hz
fillx = np.concatenate([x, x[::-1]])
x_ticks = [4, 8, 16, 32, 64]
labels= ['Pentadactyl', 'Dekadactyl', 'Dichodactyl','Cross-Pentadactyl','Cross-Dekadactyl','Cross-Dichodactyl']
condlist = ['Mon','Bin','Dich','XMon','XBin','XDich']
lineColor = ['b', 'r', 'g', 'y', 'm', 'k']
markers = ['o', 's', '^', '<', 'p', 'd']
###########################
legend_elements = []
for c, marker in zip(range(len(condlist)), markers):
a = meanF1[c, 0:5] + SEF1[c, 0:5]
b = meanF1[c, (4, 3, 2, 1, 0)] - SEF1[c, (4, 3, 2, 1, 0)]
filly = np.concatenate([a, b])
ax5.fill(fillx, filly, facecolor=lineColor[c], alpha=0.2)
ax5.plot(x, meanF1[c, 0:5], color=lineColor[c], marker=marker,markersize=10, label=condlist[c])
legend_elements.append(Line2D([0], [0], color=lineColor[c], linestyle='-', marker=marker, markersize=15, label=labels[c]))
ax5.set_ylim(0, 0.5)
y_ticks = np.arange(0, 0.6, 0.1)
formatted_y_ticks = [f'{y:.1f}' for y in y_ticks]
ax5.set_yticks(y_ticks)
ax5.set_yticklabels(formatted_y_ticks, fontsize=25)
ax5.set_xticks(x)
ax5.set_xticklabels(x_ticks, fontsize=25)
ax5.set_ylabel('Response at 26 Hz (µV)', fontsize=35)
ax5.set_xlabel('Intensity level (%)', fontsize=35)
################################################## Plot contrast response function at 23 Hz
for c, marker in zip(range(len(condlist)) , markers):
a2 = meanF2[c,0:5]+SEF2[c,0:5]
b2 = meanF2[c,(4,3,2,1,0)]-SEF2[c,(4,3,2,1,0)]
filly2 = np.concatenate([a2,b2])
ax6.fill(fillx,filly2,facecolor=lineColor[c],alpha=0.2)
ax6.plot(x, meanF2[c,0:5], color=lineColor[c], marker=marker,markersize=10, label=condlist[c])
#legend_elements.append(Line2D([0], [0], color=lineColor[c], linestyle='-', marker=marker, label=labels[c]))
leg1 = ax6.legend(
handles=legend_elements,
loc='upper left',
bbox_to_anchor=(0.01, 1),
ncol=2, # 3 columns per row
fontsize=22,
frameon=False
)
y_ticks = np.arange(0, 0.6, 0.1)
formatted_y_ticks = [f'{y:.1f}' for y in y_ticks]
ax6.set_yticks(y_ticks)
ax6.set_yticklabels(formatted_y_ticks, fontsize=25)
ax6.set_ylim(0, 0.5)
ax6.set_xticks(x)
ax6.set_xticklabels(x_ticks, fontsize=25)
ax6.set_xlabel('Intensity level (%)', fontsize=35)
ax6.set_ylabel('Response at 23 Hz (µV)', fontsize=35)
axes = [ax1, ax2, ax3, ax4, ax5, ax6]
# set the spine linewidth
for ax in axes:
for spine in ax.spines.values():
spine.set_linewidth(2.5)
ax.tick_params(axis='both', which='major', width=2)
plt.savefig('Figures/Figure2.pdf')
```
```{python fitdippers}
#| include: false
if r.fitmodels:
fitpath = "local/fits1/"
if not os.path.exists(fitpath):
os.makedirs(fitpath)
data = np.load('local/allthresh.npz')
allthresh = data['allthresh']
meanthresh = np.mean(allthresh, axis=0)
SEthresh = np.std(allthresh, axis=0, ddof=1) / np.sqrt(8)
datatofit = meanthresh
ntotalsimplexfits = 100
pedcontrasts = np.arange(-12,31,6)
pedlist = pedcontrasts
def getmodelresp(p, L, R):
Lresp = (L**p[2]) / (p[3] + L + p[5] * R)
Rresp = (R**p[2]) / (p[3] + R + p[5] * L)
bs = (Lresp**p[7] + Rresp**p[7])**(1/p[7])
resp = (bs**p[0]) / (p[4] + bs**p[1])
return resp
def discriminate(p, pedC, cond):
if cond == 1:
baseline = getmodelresp(p, pedC, 0)
elif cond == 2:
baseline = getmodelresp(p, pedC, pedC)
elif cond == 3:
baseline = getmodelresp(p, pedC, pedC)
elif cond == 4:
baseline = getmodelresp(p, pedC, 0)
modelresp = -10
contrastinc = 0
if baseline > -999:
while (modelresp - baseline) < p[6]:
contrastinc += 0.1
if cond == 1:
modelresp = getmodelresp(p, pedC + contrastinc, 0)
elif cond == 2:
modelresp = getmodelresp(p, pedC + contrastinc, pedC + contrastinc)
elif cond == 3:
modelresp = getmodelresp(p, pedC + contrastinc, pedC)
elif cond == 4:
modelresp = getmodelresp(p, pedC, contrastinc)
if contrastinc > 100:
modelresp = 999
if modelresp < 999:
while (modelresp - baseline) > p[6]:
contrastinc -= 0.001
if cond == 1:
modelresp = getmodelresp(p, pedC + contrastinc, 0)
elif cond == 2:
modelresp = getmodelresp(p, pedC + contrastinc, pedC + contrastinc)
elif cond == 3:
modelresp = getmodelresp(p, pedC + contrastinc, pedC)
elif cond == 4:
modelresp = getmodelresp(p, pedC, contrastinc)
else:
contrastinc = 99
return contrastinc
def getslope(p, pedC, cond):
targetlevelsdB = np.array(range(-18,36,3))
targetlevelsC = 10 ** (targetlevelsdB/20)
if cond == 1:
baseline = getmodelresp(p, pedC, 0)
elif cond == 2:
baseline = getmodelresp(p, pedC, pedC)
elif cond == 3:
baseline = getmodelresp(p, pedC, pedC)
elif cond == 4:
baseline = getmodelresp(p, pedC, 0)
propcorr = targetlevelsC*0
for t in range(len(targetlevelsC)):
if cond == 1:
modelresp = getmodelresp(p, pedC + targetlevelsC[t], 0)
elif cond == 2:
modelresp = getmodelresp(p, pedC + targetlevelsC[t], pedC + targetlevelsC[t])
elif cond == 3:
modelresp = getmodelresp(p, pedC + targetlevelsC[t], pedC)
elif cond == 4:
modelresp = getmodelresp(p, pedC, targetlevelsC[t])
dprime = (modelresp - baseline)/(p[6]/0.9538726)
propcorr[t] = stats.norm.cdf(dprime/np.sqrt(2), loc=0, scale=1)
level = targetlevelsdB
ncorrect = np.round(propcorr * 100)
ntrials = (targetlevelsdB*0) + 100
tofit = np.vstack((level, ncorrect, ntrials)).T
result_fit = ps.psignifit(tofit, experiment_type='2AFC')
threshval = result_fit.parameter_estimate['threshold']
slopeval = 10.3 / (result_fit.parameter_estimate['width'] / (norm.ppf(1 - 0.05) - norm.ppf(0.05)))
return slopeval
def errorfit(p):
p = 10**np.array(p)
p[1] = p[1] + 1
p[1] = min(p[1], 16)
p[0] = 1 + p[0] + p[1]
p[0] = min(p[0], 20)
p[2] = 1 + p[2]
if len(p) == 7:
p = np.append(p, 1) # append p[7]=1 Minkowski exponent of 1 in the summation model
pedlevelsC = 10**(pedlist / 20)
pedlevelsC[0] = 0
allpred = np.zeros((4, len(pedlist)))
for cond in range(1, 5):
for pedlev in range(len(pedlist)):
allpred[cond - 1, pedlev] = discriminate(p, pedlevelsC[pedlev], cond)
allpreddB = 20 * np.log10(allpred)
rms = np.sqrt(np.mean((allpreddB - datatofit) ** 2)) # return RMS error
if np.isnan(rms):
rms = 999
return rms
#%% fit linear summation model
allp = None
allrms = None
nfits = ntotalsimplexfits
###################################
def runfitting():
initial_params = np.random.normal(0, 0.1, 7 ) + np.log10([0.5, 5.5, 0.3, 1, 0.01, 1, 0.2])
sout = minimize(errorfit, initial_params, method='Nelder-Mead', tol= 1e-10)
if sout.success:
print("Optimization successful", sout.x)
allout = np.concatenate(([errorfit(sout.x)], sout.x, [0]))
filecount = 1
while True:
filepath = f"{fitpath}M1f{filecount}.npz"
if not os.path.exists(filepath):
np.savez(filepath, allout=allout)
break
filecount += 1
return allout
else:
raise ValueError("Optimization failed")
allparams = np.zeros((ntotalsimplexfits, 9))
nfits = ntotalsimplexfits
for n in range(1, ntotalsimplexfits + 1):
if os.path.exists(os.path.join(fitpath, f'M1f{n}.npz')):
nfits -= 1
if nfits > 0:
for i in range(nfits):
result = runfitting()
allparams[i] = result # Store each result in the allparams array
####################################
finalout = np.zeros((ntotalsimplexfits, 9))
for n in range(1, ntotalsimplexfits + 1):
filepath = f"{fitpath}M1f{n}.npz"
if os.path.exists(filepath):
data = np.load(filepath)
allout = data['allout']
finalout[n - 1, :] = allout
i = np.argmin(finalout[:, 0]) #Finds the row index i of the minimum value in the first column of finalout
p = 10 ** finalout[i, 1:9] #take the corresponding parameters
allrms = finalout[i, 0]
p[1] += 1
p[1] = min(p[1], 16)
p[0] = 1 + p[0] + p[1]
p[0] = min(p[0], 20)
p[2] += 1
allp = p
allparams = finalout
model1params = allp
model1rms = allrms
####Generate prediction data
pedlist = np.arange(-12, 40, 0.1)
pedlevelsC = 10 ** (pedlist / 20)
pedlevelsC[0] = 0
pedcontrasts = np.arange(-12,31,6)
pedlevelsC2 = 10 ** (pedcontrasts / 20)
pedlevelsC2[0] = 0
allpred = np.zeros((4, len(pedlist)))
allslopepred = np.zeros((4, len(pedlevelsC2)))
for cond in range(1,5):
for pedlev in range(len(pedlist)):
allpred[cond-1, pedlev] = discriminate(allp, pedlevelsC[pedlev], cond)
for cond in range(1,5):
for pedlev in range(len(pedlevelsC2)):
allslopepred[cond-1, pedlev] = getslope(allp, pedlevelsC2[pedlev], cond)
allpreddB = 20 * np.log10(allpred)
model1rms = round(model1rms,2)
print(model1rms)
print(model1params)
#%%################################ fit Minkowski summation model
allp2 = None
allrms2 = None
pedlist = np.arange(-12,31,6)
datatofit = meanthresh
nfits = ntotalsimplexfits
###################################
def runfitting2():
initial_params = np.random.normal(0, 0.1, 8 ) + np.log10([0.5, 5.5, 0.3, 1, 0.01, 1, 0.2, 4])
sout = minimize(errorfit, initial_params, method='Nelder-Mead', tol= 1e-10)
if sout.success:
print("Optimization successful", sout.x)
allout = np.concatenate(([errorfit(sout.x)], sout.x))
# Save to a unique file
filecount = 1
while True:
filepath = f"{fitpath}M2f{filecount}.npz"
if not os.path.exists(filepath):
np.savez(filepath, allout=allout)
break
filecount += 1
return allout
else:
raise ValueError("Optimization failed")
allparams2 = np.zeros((ntotalsimplexfits, 9))
nfits = ntotalsimplexfits
for n in range(1, ntotalsimplexfits + 1):
if os.path.exists(os.path.join(fitpath, f'M2f{n}.npz')):
nfits -= 1
if nfits > 0:
for i in range(nfits):
result2 = runfitting2()
allparams2[i] = result2 # Store each result in the allparams array
finalout2 = np.zeros((ntotalsimplexfits, 9))
for n in range(1, ntotalsimplexfits + 1):
filepath = f"{fitpath}M2f{n}.npz"
if os.path.exists(filepath):
data = np.load(filepath)
allout = data['allout']
finalout2[n - 1, :] = allout
i = np.argmin(finalout2[:, 0]) #Finds the row index i of the minimum value in the first column of finalout
p = 10 ** finalout2[i, 1:9] #take the corresponding parameters
allrms2 = finalout2[i, 0]
p[1] += 1
p[1] = min(p[1], 16)
p[0] = 1 + p[0] + p[1]
p[0] = min(p[0], 20)
p[2] += 1
allp2 = p
allparams2 = finalout2
#np.save(os.path.join(fitpath, 'simplexfitsM2.npy'), {'allp': allp, 'allrms': allrms, 'allparams': allparams, 'ntotalsimplexfits': ntotalsimplexfits})
model2params = allp2
model2rms = allrms2
pedlist = np.arange(-12, 40, 0.1)
pedlevelsC = 10 ** (pedlist / 20)
pedlevelsC[0] = 0
pedcontrasts = np.arange(-12,31,6)
pedlevelsC2 = 10 ** (pedcontrasts / 20)
pedlevelsC2[0] = 0
#### generate prediction data
allpred2 = np.zeros((4, len(pedlist)))
allslopepred2 = np.zeros((4, len(pedlevelsC2)))
for cond in range(1,5):
for pedlev in range(len(pedlist)):
allpred2[cond-1, pedlev] = discriminate(allp2, pedlevelsC[pedlev], cond)
for cond in range(1,5):
for pedlev in range(len(pedlevelsC2)):
allslopepred2[cond-1, pedlev] = getslope(allp2, pedlevelsC2[pedlev], cond)
allpreddB2 = 20 * np.log10(allpred2)
model2rms=round(model2rms,2)
print(model2rms)
print(model2params)
np.savez('local/dippermodels.npz', model1params=model1params, allpreddB=allpreddB, model1rms=model1rms, model2params=model2params, allpreddB2=allpreddB2, model2rms=model2rms, allslopepred=allslopepred, allslopepred2=allslopepred2)
```
```{python fitEEG}
#| include: false
if r.fitmodels:
fitpath = "local/fits2/"
if not os.path.exists(fitpath):
os.makedirs(fitpath)
EEGF1 = np.load("local/allF1.npy") # 30*6*5
meanF1 = np.mean(EEGF1, axis=0) #6*5
EEGF2 = np.load("local/allF2.npy") # 30*6*5
meanF2 = np.mean(EEGF2, axis=0) #6*5
mean = np.concatenate((meanF1, meanF2) , axis=0) #12*5 , mean data of 26 and 23Hz
def get_twoStage_model1(A, B, C, D, p): # [A, B]= 26Hz, [C, D]= 23Hz
p = np.power(10, p)
respAL = (A ** p[0]) / (p[2] + A + p[3] * B + p[3] * D)
respAR = (B ** p[0]) / (p[2] + B + p[3] * A + p[3] * C)
bsA = p[4]*(respAL + respAR)
respA = bsA + p[1]
return respA
def get_twoStage_model2(A, B, C, D, p): # [A, B]= 26Hz, [C, D]= 23Hz
p = np.power(10, p)
respBL = (C ** p[0]) / (p[2] + C + p[3] * D + p[3] * B)
respBR = (D ** p[0]) / (p[2] + D + p[3] * C + p[3] * A)
bsB = p[4]*(respBL + respBR)
respB = bsB + p[1]
return respB
def get_twoStage_modelresp(p, contrastsdB):
contrastsC = 10**(contrastsdB / 20)
responses = np.array([
# 26 Hz
get_twoStage_model1(contrastsC, 0, 0, 0, p), # Pentedactyl_resp
get_twoStage_model1(contrastsC, contrastsC, 0, 0, p), # Dekadactyl_resp
get_twoStage_model1(contrastsC, 32, 0, 0, p), # Dichodactyl_resp
get_twoStage_model1(0, 0, contrastsC, 0, p), # crossPentedactyl_resp
get_twoStage_model1(contrastsC, 0, 0, contrastsC, p), # crossDekadactyl_resp
get_twoStage_model1(contrastsC, 0, 0, 32, p), # crossDichodactyl_resp
## 23 Hz
get_twoStage_model2(contrastsC, 0, 0, 0, p), # Pentedactyl_resp
get_twoStage_model2(contrastsC, contrastsC, 0, 0, p), # Dekadactyl_resp
get_twoStage_model2(contrastsC, 32, 0, 0, p), # Dichodactyl_resp
get_twoStage_model2(0, 0, contrastsC, 0, p), # crossPentedactyl_resp
get_twoStage_model2(contrastsC, 0, 0, contrastsC, p), # crossDekadactyl_resp
get_twoStage_model2(contrastsC, 0, 0, 32, p), # crossDichodactyl_resp
])
return responses #12*5
def geterror(p):
contrastsdB = np.arange(12,37,6)
modelpredictions = get_twoStage_modelresp(p, contrastsdB) # 6*5
rms = np.sqrt(np.mean((modelpredictions - mean)**2))
return rms
allp = None
allrms = None
ntotalsimplexfits = 100
nfits = ntotalsimplexfits
def runfitting():
initial_params = [1.5, 0.05, 20, 0.5, 0.2]+ 0.01 * np.random.randn(5)
sout = minimize(geterror, np.log10(initial_params), method='Nelder-Mead', tol= 1e-10)
if sout.success:
print("Optimization successful", sout.x)
allout = np.concatenate(([geterror(sout.x)], sout.x))
filecount = 1
while True:
filepath = f"{fitpath}M1f{filecount}.npz"
if not os.path.exists(filepath):
np.savez(filepath, allout=allout)
break
filecount += 1
return allout
else:
print("Optimization did not converge.")
#raise ValueError("Optimization failed")
allparams = np.zeros((ntotalsimplexfits, 6))
nfits = ntotalsimplexfits
for n in range(1, ntotalsimplexfits + 1):
if os.path.exists(os.path.join(fitpath, f'M1f{n}.npz')):
nfits -= 1
if nfits > 0:
for i in range(nfits):
result = runfitting()
allparams[i] = result
finalout = np.zeros((ntotalsimplexfits, 6))
for n in range(1, ntotalsimplexfits + 1):
filepath = f"{fitpath}M1f{n}.npz"
if os.path.exists(filepath):
data = np.load(filepath)
allout = data['allout']
finalout[n - 1, :] = allout
i = np.argmin(finalout[:, 0])
p = finalout[i, 1:6]
contrastsfine = np.arange(12,37,0.1)
allresp = get_twoStage_modelresp(p, contrastsfine)
allrms = finalout[i, 0]
allp = 10**p
allparams = finalout
modelparams = allp
print("model 1 optimized p:", modelparams)
modelRMS = allrms
modelRMS = round(modelRMS,2)
print("model 1 rms:", modelRMS)
np.savez('local/EEGmodels.npz', allresp=allresp, modelparams=modelparams, modelRMS=modelRMS)
```
```{python plotmodelling}
#| include: false
if r.plotfigures:
pedcontrasts = np.arange(-12,31,6)
pedlist = np.arange(-12, 40, 0.1)
data = np.load('local/allthresh.npz')
allthresh = data['allthresh']
meanthresh = np.mean(allthresh, axis=0) #4*8
SEthresh = np.std(allthresh, axis=0, ddof=1) / np.sqrt(8)
data = np.load('local/dippermodels.npz')
model1params = data['model1params']
allpreddB = data['allpreddB']
model1rms = data['model1rms']
model2params = data['model2params']
allpreddB2 = data['allpreddB2']
model2rms = data['model2rms']
allslopepred = data['allslopepred']
allslopepred2 = data['allslopepred2']
fig=plt.figure(figsize=(20,28))
ax1=fig.add_subplot(3,2,1)
ax2=fig.add_subplot(3,2,2)
ax3=fig.add_subplot(3,2,3)
ax4=fig.add_subplot(3,2,4)
ax5=fig.add_subplot(3,2,5)
ax6=fig.add_subplot(3,2,6)
pos = ax1.get_position()
new_pos = [pos.x0, pos.y0+0.05, pos.width, pos.height] ###adjust position
ax1.set_position(new_pos)
pos = ax2.get_position()
new_pos = [pos.x0 + 0.08, pos.y0+0.05, pos.width, pos.height] ###adjust position
ax2.set_position(new_pos)
pos = ax3.get_position()
new_pos = [pos.x0, pos.y0, pos.width, pos.height]
ax3.set_position(new_pos)
pos = ax4.get_position()
new_pos = [pos.x0 + 0.08, pos.y0, pos.width, pos.height]
ax4.set_position(new_pos)
pos = ax5.get_position()
new_pos = [pos.x0, pos.y0-0.05, pos.width, pos.height]
ax5.set_position(new_pos)
pos = ax6.get_position()
new_pos = [pos.x0 + 0.08, pos.y0-0.05, pos.width, pos.height]
ax6.set_position(new_pos)
axes = [ax1, ax2, ax3, ax4, ax5, ax6]
line_weight = 2
for ax in axes:
for spine in ax.spines.values(): spine.set_linewidth(line_weight)
ax1.text(0.05, 0.950, '(a)', fontsize=40, transform=fig.transFigure)
ax2.text(0.55, 0.950, '(b)', fontsize=40, transform=fig.transFigure)
ax3.text(0.05, 0.625, '(c)', fontsize=40, transform=fig.transFigure)
ax4.text(0.55, 0.625, '(d)', fontsize=40, transform=fig.transFigure)
ax5.text(0.05, 0.300, '(e)', fontsize=40, transform=fig.transFigure)
ax6.text(0.55, 0.300, '(f)', fontsize=40, transform=fig.transFigure)
#%% Plot fitting result of psychophysical data
condition = [ 'Pentadactyl', 'Dekadactyl', 'Half-Dekadactyl', 'Dichodactyl']
collist = [ 'blue', 'red', 'orange', 'darkgreen']
markers=['o','s','D','^']
for cond, marker in zip(range(4), markers):
ax1.scatter(pedcontrasts, meanthresh[cond, :], color=collist[cond], marker=marker, label=condition[cond], s=70)
ax1.plot(pedlist,allpreddB[cond,:], color=collist[cond], lw=2)
ax1.set_xlabel('Baseline intensity level (%)',fontsize=30 )
ax1.set_ylabel('Threshold (%)', fontsize=30)
ax1.set_title('Linear summation model', fontsize=30)
ax1.set_xlim(-13,31)
ax1.set_ylim(-13,31)
ax1.set_yticks([-12,-6,0,6,12,18,24,30],[0.25, 0.5, 1, 2, 4, 8, 16, 32],fontsize = 25)
ax1.set_xticks([-12,-6,0,6,12,18,24,30],[0, 0.5, 1, 2, 4, 8, 16, 32] ,fontsize = 25)
ax1.text(10,-11,"RMSE = " + str(np.round(model1rms,2)) + "dB", fontsize=30)
ax1.plot([-12, 30], [-12, 30], color='black', linestyle="--")
legend_elements2 = []
for cond, marker in zip(range(4) , markers):
ax2.scatter(pedcontrasts,meanthresh[cond, :], color=collist[cond], marker=marker, label=condition[cond], s=70)
ax2.plot(pedlist,allpreddB2[cond], color=collist[cond], lw=2 )
legend_elements2.append(Line2D([0], [0], color=collist[cond], linestyle='-', label=condition[cond], markersize=12, marker=marker))
ax2.set_xlabel('Baseline intensity level (%)',fontsize=30 )
ax2.set_ylabel('Threshold (%)', fontsize=30)
ax2.set_title('Minkowski summation model', fontsize=30)
ax2.set_xlim(-13,31)
ax2.set_ylim(-13,31)
ax2.set_yticks([-12,-6,0,6,12,18,24,30],[0.25, 0.5, 1, 2, 4, 8, 16, 32],fontsize = 25)
ax2.set_xticks([-12,-6,0,6,12,18,24,30],[0, 0.5, 1, 2, 4, 8, 16, 32] ,fontsize = 25 )
leg2 = ax2.legend(
handles=legend_elements2,
loc='upper left',
bbox_to_anchor=(0.01, 1),
fontsize=20,
frameon=False
)
ax2.text(10,-11,"RMSE = " + str(np.round(model2rms,2)) + "dB", fontsize=30)
ax2.plot([-12, 30], [-12, 30], color='black', linestyle="--")
################ plot slope
data1 = np.load('local/allslope.npz')
allslope = data1['allslope']
allslope = 20 * np.log10(allslope)
allslopepred = 20 * np.log10(allslopepred)
allslopepred2 = 20 * np.log10(allslopepred2)
meanslope = np.mean(allslope, axis=0) #4*8
SEslope = np.std(allslope, axis=0, ddof=1) / np.sqrt(8)
for cond, marker in zip(range(4), markers):
ax3.scatter(pedcontrasts, meanslope[cond, :], color=collist[cond], marker=marker, label=condition[cond], s=70)
ax3.plot(pedcontrasts,allslopepred[cond,:], color=collist[cond], lw=2)
ax3.set_xlabel('Baseline intensity level (%)',fontsize=30 )
ax3.set_ylabel('Weibull ' + r"$\mathrm{\beta}$", fontsize=30)
ax3.set_xlim(-13,31)
ax3.set_ylim(-7,19)
ax3.set_xticks([-12,-6,0,6,12,18,24,30],[0, 0.5, 1, 2, 4, 8, 16, 32],fontsize = 25)
ax3.set_yticks([-6,0,6,12,18],[0.5, 1, 2, 4, 8] ,fontsize = 25)
#ax3.text(10,-11,"RMSE = " + str(np.round(model1rms,2)) + "dB", fontsize=30)
ax3.plot([-12, 30], [20*np.log10(1.3), 20*np.log10(1.3)], color='black', linestyle="--")
for cond, marker in zip(range(4), markers):
ax4.scatter(pedcontrasts, meanslope[cond, :], color=collist[cond], marker=marker, label=condition[cond], s=70)
ax4.plot(pedcontrasts,allslopepred2[cond,:], color=collist[cond], lw=2)
ax4.set_xlabel('Baseline intensity level (%)',fontsize=30 )
ax4.set_ylabel('Weibull ' + r"$\mathrm{\beta}$", fontsize=30)
ax4.set_xlim(-13,31)
ax4.set_ylim(-7,19)
ax4.set_xticks([-12,-6,0,6,12,18,24,30],[0, 0.5, 1, 2, 4, 8, 16, 32],fontsize = 25)
ax4.set_yticks([-6,0,6,12,18],[0.5, 1, 2, 4, 8] ,fontsize = 25)
#ax4.text(10,-11,"RMSE = " + str(np.round(model1rms,2)) + "dB", fontsize=30)
ax4.plot([-12, 30], [20*np.log10(1.3), 20*np.log10(1.3)], color='black', linestyle="--")
#%% Plot fitting result of EEG data
EEGF1 = np.load("local/allF1.npy") # 30*6*5
meanF1 = np.mean(EEGF1, axis=0) #6*5
EEGF2 = np.load("local/allF2.npy") # 30*6*5
meanF2 = np.mean(EEGF2, axis=0) #6*5
data = np.load('local/EEGmodels.npz')
allresp = data['allresp']
modelparams = data['modelparams']
modelRMS = data['modelRMS']
contrastsdB = np.arange(12, 37, 6)
contrastsC = 10**(contrastsdB / 20)
contrastsfine = np.arange(12, 37, 0.1)
conditions = ['Pentadactyl', 'Dekadactyl', 'Dichodactyl','Cross-Pentadactyl', 'Cross-Dekadactyl', 'Cross-Dichodactyl', 'Pentadactyl', 'Dekadactyl', 'Dichodactyl','Cross-Pentadactyl', 'Cross-Dekadactyl', 'Cross-Dichodactyl']
collist = ['b', 'r', 'g','y','m', 'k', 'b', 'r', 'g','y','m', 'k']
markers = ['o', 's', '^', '<', 'p', 'd' , 'o', 's', '^', '<', 'p', 'd']
allresp = np.load('local/EEG_allresp.npy')
# Plot 26 Hz
legend_elements = []
for i in range(6):
idx = i % len(markers)
ax5.plot(contrastsfine, allresp[i, :], color=collist[idx], lw=2)
ax5.scatter(contrastsdB, meanF1[i, :], marker=markers[idx], color=collist[idx], s=70)
legend_elements.append(Line2D([0], [0], marker=markers[idx], color=collist[idx], linestyle='-', label=conditions[i], markersize=12))
ax5.set_title("26 Hz", fontsize=30)
ax5.set_xlabel('Intensity level (%)', fontsize=30)
ax5.set_ylabel('Amplitude (μV)', fontsize=30)
ax5.set_xticks(contrastsdB)
ax5.set_xticklabels(['4', '8', '16', '32', '64'], fontsize=25)
ax5.set_yticks([0,0.1,0.2,0.3,0.4,0.5], [0,0.1,0.2,0.3,0.4,0.5], fontsize = 25)
ax5.text(12,0.48,"RMSE = " + str(np.round(modelRMS,3)) + "μV", fontsize=30)
# Plot 23 Hz
legend_elements = []
for i in range(6):
idx = i % len(markers)
ax6.plot(contrastsfine, allresp[i+6, :], color=collist[idx], lw=2)
ax6.scatter(contrastsdB, meanF2[i, :], marker=markers[idx],color=collist[idx], s=70)
legend_elements.append(Line2D([0], [0], marker=markers[idx], color=collist[idx], linestyle='-', label=conditions[i], markersize=12))
ax6.set_title("23 Hz", fontsize=30)
ax6.set_xlabel('Intensity level (%)', fontsize=30)
ax6.set_ylabel('Amplitude (μV)', fontsize=30)
ax6.set_xticks(contrastsdB)
ax6.set_xticklabels(['4', '8', '16', '32', '64'], fontsize=25)
ax6.set_yticks([0,0.1,0.2,0.3,0.4,0.5], [0,0.1,0.2,0.3,0.4,0.5], fontsize = 25)
# Add legend to ax4
leg6 = ax6.legend(
handles=legend_elements,
loc='upper left',
bbox_to_anchor=(0.005, 1),
ncol=2,
fontsize=20,
frameon=False
)
plt.savefig('Figures/Figure3.pdf')
```
```{python dostats}
#| include: false
row = pd.read_csv('local/8subjects.csv')
ex1data = row[['Subject','Condition', 'PedestalContrast','TargetContrast','IsCorrect']].values
df = pd.DataFrame(ex1data, columns=['Subject','Condition', 'PedestalContrast','TargetContrast', 'IsCorrect'])
df['Condition'] = (df['Condition'] + 1) // 2
r.ntotaltrials = ex1data.shape[0]
r.trialspersubj = int(np.round(ex1data.shape[0]/8))
data = np.load('local/allthresh.npz')
allthresh = data['allthresh'] # 8*4*8
meanthresh = np.mean(allthresh, axis=0) #4*8
r.monbinsupradB = np.round(np.mean(meanthresh[0,1:7]-meanthresh[1,1:7]),2)
r.monbinsupra = np.round(10**(r.monbinsupradB/20),2)
r.dichthdelevdB = np.round(meanthresh[3,7] - meanthresh[3,0], 2)
r.dichthdelev = np.round(10**(r.dichthdelevdB/20), 2)
merged_array = allthresh.reshape(256, 1)
# Create index arrays
subs = np.repeat(np.arange(1, 9), 32).reshape(-1, 1) # 8 participants
conditions = np.tile(np.repeat(np.arange(1, 5), 8), 8).reshape(-1, 1) # 4 conditions
levels = np.tile(np.arange(1, 9), 32).reshape(-1, 1) # 8 levels, repeated 32 times
# Combine all data into one array
alldata_array = np.hstack((subs, conditions, levels, merged_array))
df = pd.DataFrame(alldata_array, columns=['subs', 'condition', 'level', 'threshold'])
df['subs'] = df['subs'].astype(int)
df['condition'] = df['condition'].astype(int)
df['level'] = df['level'].astype(int)
anova_results = pg.rm_anova(dv='threshold', within=['condition', 'level'], subject='subs', data=df, detailed=True)
pd.set_option('display.max_columns', None) # Show all columns
pd.set_option('display.width', None) # Adjust the display width to show the full DataFrame
print(anova_results)
dipperstats = np.array(anova_results)
dipperstats2dp = np.array(anova_results.round(2))
dipperstats3dp = np.array(anova_results.round(3))
posthoc_condition = pg.pairwise_tests(dv='threshold', within='condition', subject='subs', data=df, padjust='bonf')
print("Pairwise comparisons for condition:")
print(posthoc_condition)
t_value = posthoc_condition.iloc[0, 5]
r.t_value = t_value.round(2)
r.dipperF1 = dipperstats2dp[0,5]
r.dipperF2 = dipperstats2dp[1,5]
r.dipperF3 = dipperstats2dp[2,5]
r.dipperng1 = dipperstats3dp[0,8]
r.dipperng2 = dipperstats3dp[1,8]
r.dipperng3 = dipperstats3dp[2,8]
posthoc = pairwise_ttests(dv='threshold', within=['level','condition'], padjust='bonf', data=df, subject='subs')
print(posthoc)
r.dipperposthocs = posthoc
EEGF1 = np.load("local/allF1.npy") # 30*6*5
EEGF2 = np.load("local/allF2.npy") # 30*6*5
meanF1 = np.mean(EEGF1, axis=0) #6*5
monbinratios = meanF1[1,:]/meanF1[0,:]
monbinratiosdB = 20*np.log10(monbinratios)
meanratiodB = np.round(np.mean(monbinratiosdB), 2)
r.meanratio = np.round(10**(meanratiodB/20), 2)
merged_array1 = EEGF1.reshape(np.prod(EEGF1.shape), 1)
subs = np.repeat(np.arange(1, 31), 30).reshape(-1, 1) # 30 participants
conditions = np.tile(np.repeat(np.arange(1, 7), 5), 30).reshape(-1, 1)
levels = np.tile(np.arange(1, 6), 180).reshape(-1, 1)
# Combine all data into one array
alldata_array1 = np.hstack((subs, conditions, levels, merged_array1))
df = pd.DataFrame(alldata_array1, columns=['subs', 'condition', 'level', 'amplitudes'])
df['subs'] = df['subs'].astype(int)
df['condition'] = df['condition'].astype(int)
df['level'] = df['level'].astype(int)
###ANOVA F1
aov_results1 = pg.rm_anova(dv='amplitudes', within=['condition', 'level'], subject='subs', data=df, detailed=True)
pd.set_option('display.max_columns', None) # Show all columns
pd.set_option('display.width', None) # Adjust the display width to show the full DataFrame
print(aov_results1)
#########post-hoc test
posthocF1= pairwise_ttests(dv='amplitudes', within=['level','condition'], padjust='bonf', data=df, subject='subs')
#posthocF1.to_csv('C:/documents/meanEEGresp/posthocF1.csv', index=False)
print(posthocF1)
EEGF1stats = np.array(aov_results1)
EEGF1stats2dp = np.array(aov_results1.round(2))
EEGF1stats3dp = np.array(aov_results1.round(3))
r.EEGF1_F1 = EEGF1stats2dp[0,5]
r.EEGF1_F2 = EEGF1stats2dp[1,5]
r.EEGF1_F3 = EEGF1stats2dp[2,5]
r.EEGF1_ng1 = EEGF1stats3dp[0,8]
r.EEGF1_ng2 = EEGF1stats3dp[1,8]
r.EEGF1_ng3 = EEGF1stats3dp[2,8]
merged_array2 = EEGF2.reshape(np.prod(EEGF2.shape), 1)
subs = np.repeat(np.arange(1, 31), 30).reshape(-1, 1) # 30 participants
conditions = np.tile(np.repeat(np.arange(1, 7), 5), 30).reshape(-1, 1)
levels = np.tile(np.arange(1, 6), 180).reshape(-1, 1)
# Combine all data into one array
alldata_array2 = np.hstack((subs, conditions, levels, merged_array2))
df = pd.DataFrame(alldata_array2, columns=['subs', 'condition', 'level', 'amplitudes'])
df['subs'] = df['subs'].astype(int)
df['condition'] = df['condition'].astype(int)
df['level'] = df['level'].astype(int)
###ANOVA F2
aov_results2 = pg.rm_anova(dv='amplitudes', within=['condition', 'level'], subject='subs', data=df, detailed=True)
pd.set_option('display.max_columns', None) # Show all columns
pd.set_option('display.width', None) # Adjust the display width to show the full DataFrame
print(aov_results2)
#########post-hoc test
posthocF2= pairwise_ttests(dv='amplitudes', within=['level','condition'], padjust='bonf', data=df, subject='subs')
#posthocF2.to_csv('C:/documents/meanEEGresp/posthocF2.csv', index=False)
print(posthocF2)
EEGF2stats = np.array(aov_results2)
EEGF2stats2dp = np.array(aov_results2.round(2))
EEGF2stats3dp = np.array(aov_results2.round(3))
r.EEGF2_F1 = EEGF2stats2dp[0,5]
r.EEGF2_F2 = EEGF2stats2dp[1,5]
r.EEGF2_F3 = EEGF2stats2dp[2,5]
r.EEGF2_ng1 = EEGF2stats3dp[0,8]
r.EEGF2_ng2 = EEGF2stats3dp[1,8]
r.EEGF2_ng3 = EEGF2stats3dp[2,8]
### one-way ANOVA on Crossdichodactyl condition
Crossdichodactyl_dataF2 = EEGF2[:,5,:]
merged_array3 = Crossdichodactyl_dataF2.reshape(150, 1)
subs = np.repeat(np.arange(1, 31), 5).reshape(-1, 1)
levels = np.tile(np.arange(1, 6), 30).reshape(-1, 1)
alldata_array3 = np.hstack((subs, levels, merged_array3))
df = pd.DataFrame(alldata_array3, columns=['subs', 'level', 'amplitudes'])
df['subs'] = df['subs'].astype(int)
df['level'] = df['level'].astype(int)
aov_results3 = pg.rm_anova(dv='amplitudes', within=['level'], subject='subs', data=df, detailed=True, effsize='np2')
pd.set_option('display.max_columns', None) # Show all columns
pd.set_option('display.width', None) # Adjust the display width to show the full DataFrame
print('Crossdichodactyl:', aov_results3)
EEGF2Bstats = np.array(aov_results3)
EEGF2Bstats2dp = np.array(aov_results3.round(2))
EEGF2Bstats3dp = np.array(aov_results3.round(3))
r.EEGF2B_F1 = EEGF2Bstats2dp[0,4]
r.EEGF2B_ng1 = EEGF2Bstats3dp[0,7]
data = np.load('local/EEGmodels.npz')
temp = data['modelparams']
r.EEGparams_m = np.round(temp[0],3)
r.EEGparams_S = np.round(temp[2],3)
r.EEGparams_w = np.round(temp[3],3)
r.EEGparams_k = np.round(temp[4],3)
r.EEGparams_R = np.round(temp[1],3)
temp = data['modelRMS']
r.EEGRMS = np.round(temp,3)
data = np.load('local/dippermodels.npz')
temp = data['model1params']
r.model1params_p = np.round(temp[0],3)
r.model1params_q = np.round(temp[1],3)
r.model1params_m = np.round(temp[2],3)
r.model1params_S = np.round(temp[3],3)
r.model1params_Z = np.round(temp[4],3)
r.model1params_w = np.round(temp[5],3)
r.model1params_k = np.round(temp[6],3)
r.model1params_y = np.round(temp[7],3)
temp = data['model1rms']
r.model1rms = np.round(temp,2)
temp = data['model2params']
r.model2params_p = np.round(temp[0],3)
r.model2params_q = np.round(temp[1],3)
r.model2params_m = np.round(temp[2],3)
r.model2params_S = np.round(temp[3],3)
r.model2params_Z = np.round(temp[4],3)
r.model2params_w = np.round(temp[5],3)
r.model2params_k = np.round(temp[6],3)
r.model2params_y = np.round(temp[7],3)
temp = data['model2rms']
r.model2rms = np.round(temp,2)
```
# Abstract
While the brain’s integration of auditory and visual inputs has been extensively investigated, the mechanisms underlying the combination of somatosensory signals remain less explored. Here, we investigated vibrotactile summation across fingers using psychophysical and electrophysiological methods. In Experiment 1, discrimination thresholds for 26 Hz vibrations were measured using a two-interval forced-choice (2IFC) task. Thresholds exhibited a 'dipper' pattern when plotted against baseline intensity level. Detection thresholds decreased by approximately 1dB when all ten fingers were stimulated ('dekadactyl' condition) simultaneously compared to when each alternate finger was stimulated ('pentadactyl' condition), suggesting a process of probability summation. When a target stimulus was presented to five digits and a baseline stimulus to the remaining digits (‘dichodactyl’ condition) simultaneously, thresholds increased, consistent with suppression between digit representations. In Experiment 2, steady-state somatosensory evoked potential (SSSEP) signals showed approximately a 1.4-fold amplitude increase in dekadactyl compared to pentadactyl conditions, indicating a summation effect. Conversely, SSSEP amplitudes decreased when targets and masks were presented simultaneously at different frequencies, providing additional evidence of suppression. These results are consistent with a model featuring inhibition between digits and reveal that the weight of suppression is intermediate between that observed in binocular vision and binaural hearing.
## Keywords:
*vibrotactile summation*, *suppression*, *somatosensory*, *computational modelling*
# Significance
Suppression is known to occur both between the eyes and the ears, but how does the brain combine tactile inputs from multiple fingers? Using psychophysical threshold measurements and EEG recordings, we found that doubling the inputs led to a weak improvement in threshold and an increase in EEG amplitude, though less than double. This suggests that the combination of vibration signals involves both probability summation and suppressive processes. Computational modelling further suggests a unique mechanism of tactile integration, distinct from those observed in vision and audition.
# Introduction
Signal combination is crucial for human perception and our interaction with the environment. Our sensory systems comprise a variety of receptors and neural pathways that detect and process a wide range of stimuli, including sight, hearing, touch, and smell. To construct a coherent and comprehensive representation of our surroundings and our own body, the brain must filter out overlapping or redundant information and avoid excessively strong signals to prevent sensory overload. Therefore, in addition to additive processes, suppressive processes are also involved in signal combination. This suppression during signal combination has been investigated both within [e.g., vision, hearing, touch, @Baker2017; @Baker2020; @Biermann1998; @Gandevia1983] and between [e.g., visuo-tactile, audio-visual, @Ide2013; @Hidaka2015] modalities.
In visual and auditory perception, psychophysical studies have shown that detection performance for binocular or binaural presentation is typically between a factor of $\sqrt{2}$ and 2 better than for monocular or monaural presentation [@Baker2018; @Baker2020; @Campbell1965]. This summation at threshold implies the existence of physiological mechanisms that combine signals across eyes or ears. Above threshold, detection performance improves when a weak fixed intensity (‘baseline’) stimulus is added, a phenomenon known as facilitation. At higher baseline intensities performance worsens, producing a masking effect. When plotted against baseline intensity, the discrimination thresholds exhibit a "dipper" shape (see [@fig-methods]b). The facilitation and masking effects arise from the brain transducing physical signals into neural responses in a non-linear manner, and are observed for many sensory stimuli [reviewed by @Solomon2009].
To complement psychophysical work, we can directly measure the response to stimuli of different intensities by recording brain activity. One convenient method is the steady state technique, in which periodic stimulus oscillations are reflected in electromagnetic neural responses at the same frequency, which can be recorded using EEG or MEG. For instance, some EEG studies have investigated the signal combination process in visual and auditory modalities by recording steady-state signals [@Baker2017; @Baker2020]. These studies found that responses increased when input channels (eyes or ears) were doubled, but by less than a factor of two. Additionally, when a mask was added to one channel instead of a signal input (i.e. oscillating at a different frequency), the signal response decreased due to suppression [@Busse2009].
The process by which the brain combines multiple signals has attracted great interest in recent years, leading to the development of several computational models. The two stage gain control model of signal combination [@Meese2006] successfully accounts for both binocular and binaural perception, positing that the signal from one channel (left or right) is inhibited by the other channel before being summed. The equations describing this model are defined as:
$$Stage1_L = \frac{I^m_L}{S + I_L + \omega I_R}$$ {#eq-stage1L}
$$Stage1_R = \frac{I^m_R}{S + I_R + \omega I_L}$$ {#eq-stage1R}
$$sum = Stage1_L + Stage1_R$$ {#eq-binsum}
$$resp = \frac{sum^p}{Z + sum^q}$$ {#eq-stage2}
In this model, sensory inputs from adjacent channels are processed in two stages. At the first stage, each input undergoes gain control, where the response from one channel is suppressed by the input from the other channel. This is described by Equations 1 and 2. For example, the left channel response ($Stage1_L$) depends not only on its own stimulus intensity ($I_L$) but also on the intensity of the right channel ($I_R$), modulated by the suppression weight $\omega$. The exponent $m$ and saturation constant $S$ are free parameters and determine how the input is transformed nonlinearly. The outputs of the left and right channels are then summed to form $sum$ (Equation 3), which represents the combined signal strength after inter-channel suppression. This combined signal is then passed through a nonlinear transducer function (Equation 4), which produces the final response ($resp$) based on the free parameters $p$, $q$ and $Z$. These parameters govern the gain and saturation characteristics of the system. In vision, the weight of suppression is approximately $\omega = 1$, whereas in auditory perception, suppression between the ears is dramatically weaker, with $\omega$ close to 0. Therefore, while suppression can be key to signal combination, its strength varies across different sensory modalities.
In tactile perception, some neuroimaging findings suggest there exist summation and suppression of vibration signals also. For instance, studies have recorded brain responses to vibration stimuli delivered to two fingers separately and simultaneously. Results indicated that at low levels of stimulation near the perceptual threshold, there was a facilitation or summation effect between two fingers. However, at higher levels of stimulation, a suppression effect was observed [@Gandevia1983]. Furthermore, the brain's response to simultaneous vibration of two fingers is less than (approximately 50\% of) the sum of the responses to individual finger stimuli [@Biermann1998; @Gandevia1983]. In addition, the extent of suppression is dependent on the spatial distance between the fingers [@Gandevia1983; @Hoechstetter2001]. In some psychophysical studies, researchers measured vibration detection thresholds using various contactor sizes [@Gescheider2005; @Gu2013]. Results showed that the threshold decreased as the contactor size increased, a phenomenon known as area summation. Specifically, when the contactor size was doubled, the threshold decreased by approximately 3 dB (a factor of 1.4 @Gescheider2005). These findings focused on the summation rather than suppression between inputs, because suppressive processes are minimal at low intensities near threshold. In summary, although substantial evidence suggests the existence of suppression in vibration combination, there has been relatively little work measuring psychophysical threshold, and a corresponding computational model has yet to be developed.
Here, we combine psychophysical measurements (Experiment 1) with EEG recordings (Experiment 2) to assess human vibration perception, and fit a computational model to interpret the processes of vibration signal combination and suppression. In Experiment 1, we measure how thresholds vary across different conditions and baseline levels using a two-interval forced-choice (2IFC) task. Our aim is to investigate two processes: the summation effect (by doubling the number of stimulated digits), and the masking effect (by introducing an additional mask stimulus) (see *Psychophysical procedures*). Thresholds decreased by approximately 1 dB when the number of stimulated digits was doubled, a finding consistent with probability summation rather than physiological summation. In psychophysical tasks, probability summation refers to improved detection performance resulting from multiple independent detection opportunities at the decision level [@Tyler2000], whereas physiological summation reflects the neural integration of sensory signals across populations of neurons. Our computational modelling indicated that this effect could potentially be explained without requiring suppression between digits. This is because psychophysical thresholds are determined primarily by the most sensitive subset of neurons relevant to the task, and suppression does not necessarily impact this 2IFC task. In comparison, EEG recordings reflect the overall activity of the entire neuronal population, where suppression between neurons is an important characteristic of sensory perception. Therefore, in our EEG Experiment, we tested these neural mechanisms directly by measuring steady-state somatosensory evoked potentials (SSSEPs). We examined summation by stimulating adjacent digits at the same frequency, and suppression by introducing a second frequency simultaneously (see *EEG procedures*). These different frequencies allowed us to isolate suppressive interactions between digits without the complicating factor of summation. The EEG modelling results indicated an inhibitory weight of suppression of $\omega \sim 0.5$, intermediate between the weights observed for vision ($\omega \sim 1$) and audition ($\omega \sim 0$).
# Results
## Summation and suppression effects on vibration thresholds
In Experiment 1, participants were required to report which interval included the target stimulus. Individual thresholds were measured using a 3-down-1-up staircase procedure. These data were used to fit psychometric functions, estimating thresholds at 75% correct and slope parameters via cumulative Gaussians. The average thresholds (geometric means) across 8 participants for each condition and baseline intensity level (equivalent to the pedestal contrast in studies of visual contrast discrimination) are shown in [@fig-dippers]a. A 4 (condition) $\times$ 8 (baseline intensity level) repeated measures ANOVA was used to assess statistical differences between factors. We found significant main effects of condition (F = `r dipperF1`, $p$ < 0.001, $\eta_G^2$ = `r dipperng1`, Greenhouse-Geisser corrected) and baseline intensity level (F = `r dipperF2`, $p$ < 0.001, $\eta_G^2$ = `r dipperng2`, Greenhouse-Geisser corrected), as well as a significant interaction between the two factors (F = `r dipperF3`, $p$ < 0.001, $\eta_G^2$ = `r dipperng3`, Greenhouse-Geisser corrected). When plotting the thresholds against baseline intensity level, except for the dichodactyl condition (green diamonds), the other three conditions exhibited a dipper shape. Specifically, at low baseline intensity level (0.5 to 2\%), thresholds decreased with increasing baseline intensity level, indicating a facilitation effect. When the baseline intensity level exceeded 2\%, thresholds increased due to a masking effect, resulting in approximately parallel 'handles' across the pentadactyl (blue circles), dekadactyl (red squares) and half-dekadactyl (orange triangles) conditions. At detection threshold (baseline intensity level = 0\%), the threshold for the dekadactyl condition was around 1 dB lower than for the pentadactyl condition, with no statistically significant difference. This weak summation between digits is lower than that typically attributed to physiological summation (~3-6 dB), suggesting a process of probability summation instead [@Quick1974; @Tyler2000].
::: {.content-visible unless-format="docx"}
![Thresholds (a) and psychometric slope values (b) from Experiment 1. Data are averaged across N=8 participants, with shaded regions indicating ±1SE across participants.](Figures/Figure1.pdf){#fig-dippers}
:::
Above detection threshold, thresholds in the dekadactyl condition were significantly lower (t$_7$ = `r t_value`, $p$ < 0.01, Bonferroni-corrected) than those in the pentadactyl condition. On average, the thresholds decreased by a factor of approximately `r monbinsupra` (`r monbinsupradB` dB). Additionally, thresholds in the half-dekadactyl condition, where the baseline stimulus vibrated ten digits, but the target stimulus only vibrated five digits, were higher than those in the dekadactyl condition. These findings could be interpreted as a summation effect when the target inputs were doubled. However, the psychometric function is linearised from the bottom of the dip onwards, which might increase the effects of probability summation [@Tyler2000]. We return to this point in the modelling section. In the dichodactyl condition, thresholds increased across all baseline intensity levels, with thresholds elevated by a factor of `r dichthdelev` (`r dichthdelevdB` dB) at the highest baseline intensity level (32\%). This result is consistent with a suppression effect when adding a mask between digits, where the target stimulus could only be detected when it approached the baseline intensity.
The slopes of the psychometric function for each condition are illustrated in [@fig-dippers]b. At detection threshold, all conditions exhibited relatively steep slopes, around $\beta = 4$. At low baseline intensity levels (0 to 4\%), except for in the dichodactyl condition (green triangles), the slopes decreased with increasing baseline intensity level, approaching $\beta = 1$ at a baseline intensity level of 4\%. At higher baseline intensity levels (4\% to 32\%), the slopes for these three conditions remained shallow, consistent with previous studies using the same paradigm in other senses [@Foley1981; @Meese2006]. In contrast, the slopes in the dichodactyl condition remained steep across all baseline intensity levels, with slight variations in the range $2 < \beta < 5$. This pattern differs from that observed in dichoptic masking in vision [@Meese2006; @Baker2024], where slopes were steep at detection threshold, became shallow at lower contrast levels, and then became very steep again at higher contrast levels. The very steep slopes ($\beta \sim 6$) in dichoptic masking are attributed to mandatory physiological summation between the eyes [@Baker2013]. Therefore, the approximately constant slopes observed in the dichodactyl condition suggest that mandatory physiological summation between digits may not occur here.
## Summation and suppression of neural responses
The average amplitude spectra and scalp distributions for four conditions of Experiment 2 are shown in [@fig-EEG]a-d. The steady-state EEG signals were strongest at fronto-central electrodes for both frequencies, aligning with findings from previous studies involving finger vibration [@Porcu2014; @Timora2018]. Therefore, we averaged the EEG responses across six electrodes ($F1$, $F2$, $Fz$, $FC1$, $FC2$, $FCz$) to calculate intensity-response functions at both 23 Hz and 26 Hz.
::: {.content-visible unless-format="docx"}
![Results of the EEG experiment. Panels (a-d) show Fourier spectra and inset scalp topographies for a subset of four conditions at the highest intensity. Panels (e,f) show intensity-response functions for all conditions at 26Hz (e) and 23Hz (f). Data are averaged across N=30 participants, and shaded regions indicate ±1SE across participants.](Figures/Figure2.pdf){#fig-EEG}
:::
In the pentadactyl and dekadactyl conditions, where a 26 Hz stimulus was used, we observed a peak at 26 Hz in both [@fig-EEG]a and b, with a stronger activity at the fronto-central region in the dekadactyl condition. A similar pattern was observed in the cross-pentadactyl condition, where a 23 Hz stimulus resulted in a peak at 23 Hz ([@fig-EEG]c). In the cross-dekadactyl condition ([@fig-EEG]d), when two different frequencies were simultaneously presented, peaks were observed at the original frequencies (26 Hz (F1) and 23 Hz (F2)) as well as at specific intermodulation frequencies (20 Hz ($2 \times F2 - F1$) and 29 Hz ($2 \times F1 - F2$)). Additionally, the distribution of activity across the scalp was similar in the cross-dichodactyl and cross-dekadactyl conditions.
Responses increased monotonically as a function of vibration intensity for most conditions ([@fig-EEG]e,f). Two 6 (condition) $\times$ 5 (intensity) repeated measures ANOVAs were conducted on the EEG amplitudes at each frequency separately. At 26 Hz, significant main effects of condition (F = `r EEGF1_F1`, $p$ < 0.001, $\eta_G^2$ = `r EEGF1_ng1`, Greenhouse-Geisser corrected), intensity level (F = `r EEGF1_F2`, $p$ < 0.001, $\eta_G^2$ = `r EEGF1_ng2`, Greenhouse-Geisser corrected), and a significant interaction effect (F = `r EEGF1_F3`, $p$ < 0.001, $\eta_G^2$ = `r EEGF1_ng3`) were observed. In the pentadactyl condition (blue circles), EEG amplitudes increased monotonically with increasing intensity level ([@fig-EEG]e). In the dekadactyl condition (red squares), doubling the inputs led to stronger responses. Compared to the pentadactyl condition, average response amplitudes were higher by a factor of `r meanratio` across all intensity levels—reflected as a multiplicative scaling of the red line. This suggests a sub-linear summation, where the combined response is greater than in the pentadactyl condition, but less than would be expected from perfect linear summation. If adjacent digit inputs were combined linearly without suppression, we would expect a doubling of response amplitude (i.e., a factor of ~2). In the dichodactyl condition (green triangles), the target stimulus was presented to Set A, while a fixed vibration at 32\% of maximum intensity was presented to Set B, resulting in a high baseline response, and responses that increased with intensity level. At higher intensity levels (32 to 64\%), the responses in the dichodactyl condition approximated those in the dekadactyl condition. In the cross-dekadactyl (pink pentagons) and cross-dichodactyl (black diamonds) conditions, two different frequencies were presented simultaneously. Instead of summing, they suppressed each other, leading to weaker responses compared to the pentadactyl condition.
At 23 Hz ([@fig-EEG]f), we also found significant main effects of condition (F = `r EEGF2_F1`, $p$ < 0.001, $\eta_G^2$ = `r EEGF2_ng1`, Greenhouse-Geisser corrected), intensity level (F = `r EEGF2_F2`, $p$ < 0.001, $\eta_G^2$ = `r EEGF2_ng2`), as well as a significant interaction effect (F = `r EEGF2_F3`, $p$ < 0.001, $\eta_G^2$ = `r EEGF2_ng3`, Greenhouse-Geisser corrected). Further evidence of suppression was observed in the cross-dekadactyl condition (pink pentagons), where simultaneous 23 Hz and 26 Hz stimuli led to reduced responses compared to the cross-pentadactyl condition (yellow triangles). In the cross-dichodactyl condition (black diamonds), we observed a strong initial response to the 23 Hz mask at 4\% target intensity level. However, as the intensity of the 26 Hz target increased, the overall response decreased (F = `r EEGF2B_F1`, $p$ < 0.001, $\eta_G^2$ = `r EEGF2B_ng1`, Greenhouse-Geisser corrected). This negative trend suggests that the stronger the 26 Hz target, the more it suppressed the response to the 23 Hz mask [@Busse2009]. The resulting pattern reflects a masking effect and inter-digit suppression, rather than summation, leading to weaker responses as intensity levels increase. As expected, the three conditions involving only a 26 Hz vibration did not show any measurable responses at 23 Hz ([@fig-EEG]f).
## Computational modelling results
We began by fitting the data of Experiment 1 using the model of Meese et al. [@Meese2006], with seven free parameters including the weight of suppression between channels ($\omega$ in @eq-stage1L and @eq-stage1R). Because the combination across channels is linear (following an initial transducer stage), we refer to this as the Linear summation model. The model gave a reasonable fit, as shown in [@fig-models]a, with an RMS error of `r model1rms`dB (see @tbl-parametertable for parameter values). However, two systematic shortcomings are apparent. At detection threshold, the model overestimates the amount of summation (notice the red curve is below the data points for the first three baseline intensity levels). Also, suppression between channels is relatively strong ($\omega$ = `r model1params_w`), but this results in the slope of the dichodactyl condition (green curve) being steeper than in the empirical data (green triangles).
To address these shortcomings, we added a further parameter to the model - an exponent at the summation stage (see @eq-binsum2), which is the defining feature of the Minkowski summation model. The Minkowski exponent determines how multiple sensory inputs are combined at the perceptual level. At high values, it serves as an approximation of probability summation, which has traditionally been modelled using a MAX operator within a 2IFC signal detection framework [@Tyler2000]. As the Minkowski exponent increases, the model shifts from linear integration toward a MAX-like operation, where the strongest input dominates the response. This provides a flexible way to capture both linear and nonlinear integration, depending on the value of the exponent. The Minkowski summation model provides a much more satisfactory fit (see [@fig-models]b), with an RMS error of `r model2rms`dB, and no obvious systematic issues with the fit. The Minkowski exponent was estimated at $\gamma$ = `r model2params_y` - substantially greater than the implicit value of $\gamma$ = 1 in the Linear summation model, and consistent with the low levels of summation at threshold in the data. Several other parameter values also differ between the models (see @tbl-parametertable), most notably the weight of suppression ($\omega$) reduces to near-zero. With this architecture, threshold elevation in the dichodactyl condition is caused by the MAX-like operation from the high Minkowski exponent (a true MAX operator would have an exponent of $\gamma = \infty$), meaning that explicit suppression is not required.
To predict the psychometric slopes, we used the estimated parameters from the threshold fitting for each of the two models separately. For each model, we simulated model responses across a range of target stimulus levels. For each target level and condition, we calculated d-prime $(d^{\prime})$ as the signal difference divided by the internal noise, and obtained the predicted probability of a correct response by passing $d^{\prime}$ divided by $\sqrt{2}$ through the cumulative normal distribution. The resulting predicted psychometric function was fitted with a cumulative Gaussian using $psignifit$ 4, and the predicted slope parameters were extracted from the fitted functions. The results of psychometric slope estimates are shown in [@fig-models]c and 3d. Both models capture the general features of the slope data. Slopes start steep ($\beta = 4$) at detection threshold, and are linearised to around $\beta = 1.3$ at higher baseline intensities in all but the dichodactyl condition. The dichodactyl slopes remain steep ($\beta > 2$) across the full range of baseline intensities. We note that the additional nonlinearity in the Minkowski model does cause the slope values around detection threshold to be slightly overestimated, but the general trend is correct. Because these slope values are predicted with no additional free parameters, based on the fits at detection threshold (see also [@Meese2006]), this increases our confidence in the accuracy of the models.
We also performed maximum likelihood fits to the raw proportion correct data, that constrained the model to the full psychometric function (i.e. capturing both threshold and slope). These fits, summarised in the Supporting Information Appendix, resulted in a consistent pattern of values for the key parameters ($\omega$ and $\gamma$), and qualitatively similar behavior of the two models.
The absence of suppression in the Minkowski model does not necessarily mean that there is no suppression between digits - since psychophysical responses are assumed to be based on the most sensitive subset of neurons relevant to a given task, it is possible that they directly tap mechanisms that are suppression-free. Our EEG data measure responses from the whole neural population, and so we next model these results to obtain a more general estimate of suppression. The model used to fit the EEG data omits the output nonlinearity (@eq-stage2), and adds a scaling parameter ($R_{max}$), resulting in 5 free parameters. It produced an excellent account of the results at both temporal frequencies, and across all 6 stimulus conditions (see [@fig-models]e,f), with an RMS error of `r EEGRMS`$\mu$V. The model correctly captures the sub-linear summation between channels (i.e. the red curve in [@fig-models]e is less than a factor of two greater than the blue curve), and the suppression evident in the cross-dekadactyl and cross-dichodactyl conditions. The weight of suppression ($\omega$ = `r EEGparams_w`) was intermediate between the two dipper models, as well as being in between previously published values for vision ($\omega$ = 1; [@Baker2017]) and hearing ($\omega$ = 0; [@Baker2020]).
```{r}
#| label: tbl-parametertable
#| output: true
#| tbl-cap: 'Summary of fitted model parameters for three computational models. Brackets indicate fixed values, and dashes denote parameters not included in a given model.'
tabledata <- matrix(0,nrow=3,ncol=11)
tabledata <- data.frame(tabledata)
tabledata[1,3] <- as.character(model1params_p)
tabledata[1,4] <- as.character(model1params_q)
tabledata[1,5] <- as.character(model1params_m)
tabledata[1,6] <- as.character(model1params_S)
tabledata[1,7] <- as.character(model1params_Z)
tabledata[1,8] <- as.character(model1params_w)
tabledata[1,9] <- as.character(model1params_k)
tabledata[1,10] <- paste0('(',as.character(model1params_y),')')
tabledata[2,3] <- as.character(model2params_p)
tabledata[2,4] <- as.character(model2params_q)
tabledata[2,5] <- as.character(model2params_m)
tabledata[2,6] <- as.character(model2params_S)
tabledata[2,7] <- as.character(model2params_Z)
tabledata[2,8] <- as.character(model2params_w)
tabledata[2,9] <- as.character(model2params_k)
tabledata[2,10] <- as.character(model2params_y)
tabledata[3,2] <- as.character(EEGparams_R)
tabledata[3,5] <- as.character(EEGparams_m)
tabledata[3,6] <- as.character(EEGparams_S)
tabledata[3,8] <- as.character(EEGparams_w)
tabledata[3,9] <- as.character(EEGparams_k)
tabledata[,11] <- c(as.character(model1rms), as.character(model2rms), as.character(EEGRMS))
tabledata[,1] <- c('Linear summation', 'Minkowski summation', 'EEG model')
tabledata[1:2,2] <- '-'
tabledata[3,c(3,4,7,10)] <- '-'
colnames(tabledata) <- c('Model',"$R_{max}$",'$p$','$q$','$m$','$S$','$Z$',"$\\omega$",'$k$',"$\\gamma$",'RMSE')
kableExtra::kbl(tabledata, align='lcccccccccc',booktabs=TRUE,linesep='',escape = FALSE) %>% kable_styling(font_size = 8)
# %>% kableExtra::kable_styling(latex_options = "HOLD_position")
```
::: {.content-visible unless-format="docx"}
![Summary of model fitting. Panel (a) shows the fit of a model in which summation is linear (following an early transducer), and panel (b) shows the fit of the Minkowski summation model. Panels (c) and (d) show the slope estimates for the Linear summation mode (c) and Minkowski summation model (d). Panels (e) and (f) show the best model fit to the EEG data at 26Hz (c) and 23Hz (d).](Figures/Figure3.pdf){#fig-models}
:::
# Discussion
In this study, we present evidence detailing the processes of signal summation and suppression in vibrotactile signal combination. In the psychophysical experiment, we observed that doubling the inputs resulted in a slight reduction in detection and discrimination thresholds, indicating a weak summation effect. Additionally, adding a masking stimulus led to increased discrimination thresholds, indicating a suppression process. Furthermore, the EEG experiment demonstrated that adding inputs at the same frequency increases brain activity, whereas adding inouts at a different frequency decreases it. These findings allow us to quantify the summation and suppression effect between digits. Our computational model fitting results show that at a population level, responses to two vibrotactile inputs suppress each other, though our psychophysical data are more consistent with probability summation rather than neural/physiological summation. The weight of suppression between digits is around $\omega=0.5$, which is intermediate between the corresponding values observed in visual [@Baker2017] and auditory [@Baker2020] signal combination. Overall, our results suggest the presence of suppression in vibrotactile signal combination, with a suppression effect distinct from that observed in visual and auditory signal combination. In the remainder of this Discussion, we consider the underlying mechanisms of vibrotactile signal combination and the differences across sensory modalities.
Our results show that the detection threshold in the dekadactyl condition was \sim 1 dB lower than in the pentadactyl condition, although this difference did not reach statistical significance. Nevertheless, the trend is consistent with previous studies [@Gescheider2002; @Gescheider2005; @Verrillo1963] suggesting that doubling inputs can reduce detection thresholds. However, the summation effect observed in our study was weaker than the \sim 3 dB reductions reported in those studies, which used higher-frequency vibrations (>40 Hz) and doubled the contactor size. This discrepancy may be attributed to the different vibration frequencies used in the studies, and our choice to stimulate alternating fingers rather than to vary stimulus area. Some psychophysical studies have proposed that four distinct mechanoreceptive channels mediate vibration perception: Pacinian (P), Non-Pacinian I (NP I), Non-Pacinian II (NP II), and Non-Pacinian III (NP III) [@Bolanowski1988; @Verrillo1963]. Each channel has unique characteristics and responds to specific frequency ranges. For instance, the P channel is more sensitive to higher frequencies (e.g., greater than 40 Hz) and is the only channel that possesses the property of spatial summation [@Bolanowski1988; @Gescheider2002; @Verrillo1965]. In contrast, vibrations between 2 and 40 Hz are mediated by the NP I channel, which has been shown not to exhibit spatial summation [@Bolanowski1988; @Gescheider1994; @Verrillo1963]. Consequently, the weak summation effect observed in our study may be attributed to the 26 Hz vibration, which is determined by the NP I channel that lacks spatial summation. Additionally, unlike the study by @Gescheider2005, which used varying contactor sizes to stimulate random areas within a test region, thereby preventing participants from knowing the exact location of stimulation, our study specifically involved the stimulation of individual digits, with participants aware of which digits were being stimulated in some conditions. This methodological difference likely contributed to the observed discrepancy between the studies.
Our EEG experiment showed a significant increase in neural responses when doubling the number of stimulated digits at the same frequency, although the increase was less than the twofold change expected if the digits were completely independent. These findings correspond with previous studies, demonstrating that brain responses to concurrent stimuli at the same frequencies are stronger than those elicited by individual stimuli, suggesting a summation effect between inputs [@Hoechstetter2001; @Gandevia1983; @Ishibashi2000]. However, the observed response is less than the anticipated summed response, generally falling between 10\% and 50\% of the expected value, signifying a suppression effect between inputs [@Gandevia1983; @Ching-Liang1995]. In contrast, adding a masking stimulus at a different frequency resulted in a reduction in brain activity. These results provide further evidence of suppression between digits and are consistent with previous neuroimaging studies investigating the interaction between different frequencies [@Severens2010; @Pang2015]. For instance, @Severens2010 examined the suppression effect by vibrating two fingers at the same (18 Hz) or different frequencies (18 vs. 22 Hz, 18 vs. 26 Hz). Their findings revealed that the event-related potential (ERP) and SSSEP responses to simultaneous stimulation of both fingers were significantly lower than the linear sum of individual responses, with no difference in suppression based on frequency variation. This ruled out the hypothesis that suppression is due to neuronal occlusion from overlapping cortical areas (which support adjacent areas of skin), as stronger suppression would be expected at the same frequency if this were the case [@Gandevia1983; @Krause2001]. Instead, the results support the alternative explanation of lateral inhibition, where activation of one cortical neuron suppresses neighbouring neurons’ activity [@Laskin1979; @Brumberg1996; @DiCarlo1998; @Dykes1984; @Mirabella2001]. This inhibitory mechanism has been extensively observed in both animal and human studies [@Iwamura1978; @Iwamura1985; @Biermann1998; @Gandevia1983; @Hoechstetter2001].
Another interesting finding is that we observed an intermediate level of suppression between vibrotactile stimuli compared with vision and audition. In visual perception, the forward-facing positioning of the eyes results in substantially overlapping visual fields. To merge these overlapping inputs from each eye (monocular vision) into a cohesive single percept (binocular single vision), neural signals from both eyes must exhibit strong mutual inhibition to achieve ‘ocularity invariance’, ensuring the constancy of perception through one or both eyes [@Baker2007]. This is consistent with recent findings indicating that some neurons in the primary visual cortex are monocularly excitable, responding exclusively to the dominant eye, while binocular stimulation can suppress their activity [@Dougherty2019]. In contrast, the lateral placement of the ears results in minimal overlap between auditory inputs, reducing the necessity for strong interaural inhibition to integrate these signals into a unified percept, as is required in visual processing. Instead, binaural perception benefits from the comparison and integration of these disparate inputs to localise sound sources and discern subtle differences in timing and intensity [@Brugge1973; @Benson1976]. Although prior evidence supports the existence of binaural suppression [@Tiihonen1989; @Gransier2017; @Baker2020], the extent of this suppression may differ depending on hemispheric laterality [@Kaneko2003; @Fujiki2002].
Tactile perception involves the stimulation of each of the ten fingers individually, creating a more complex signal integration compared to binocular or binaural perception. The primary somatosensory cortex (S1) is divided into four distinct Brodmann areas (BA 3a, 3b, 1, and 2), each responsible for mapping different parts of the body surface. Evidence suggests that receptive field of each finger is arranged in somatotopic order in BA 3b, while in BA 1/2, they are arranged in clusters with large overlap [@Kurth2000; @Krause2001; @Ishibashi2000]. Although SSSEPs reflect activity across the whole scalp, EEG [@Pang2015; @Arslanova2022] and MEG [@Ishibashi2000; @Hoechstetter2001; @Tame2015] studies also support overlapping finger representations in S1. Notably, the overlap between adjacent fingers is greater than between non-adjacent ones, indicating that suppression between fingers depends on their spatial proximity. Additionally, hand posture plays a significant role in signal summation and suppression. For instance, a study [@Hamada2003] found that a "CLOSE" posture (thumb and index finger positioned as if to pick something up) produces stronger suppression than an "OPEN" posture (hand fully open). In the present study, finger posture and inter-finger distance were controlled by placing the fingertips in a fixed position on solenoids, suggesting that the amount of suppression may vary depending on finger arrangements or hand posture. Furthermore, studies have shown that blind individuals exhibit stronger somatosensory evoked potentials compared to sighted individuals, potentially indicating an expanded S1 region [@Giriyappa2009; @Burton2004]. This suggests that use-dependent cortical reorganisation may occur across different body regions, and the amount of suppression may vary depending on the body part involved.
## Conclusion
We demonstrate that vibrotactile signal combination involves both probability summation and suppression, as evidenced by psychophysical thresholds and EEG activity. A computational model further characterizes this process, indicating that the strength of suppression between digits is weaker than in vision but stronger than in audition. These findings establish a foundational framework for understanding vibrotactile integration and provide a basis for future research to explore whether the weight of suppression varies across different body parts, offering deeper insights into somatosensory processing and clinical conditions in which tactile processing or body experience is affected (e.g., chronic hand pain).
# Materials and Methods
## Participants
Eight adult subjects participated in the psychophysics experiment, and thirty-one adult subjects participated in the EEG experiment. All participants self-reported as healthy, with no diagnosed neurological disorders, and no history of exposure to severe hand-transmitted vibration. Both experiments were approved by the Ethics Committee of the Department of Psychology at the University of York (application IDs 2277 and 2303). Written informed consent was obtained from all participants prior to conducting the experiments.
## Apparatus & stimuli
Vibration stimuli were generated by a specially constructed board with ten fixed solenoids ('tactor' devices, from Dancer Design Ltd., Yorkshire, UK) controlled by a computer ([@fig-methods]c). Each solenoid could be independently driven by a pair of 6-channel USB sound cards to vibrate each finger. The outputs of the sound cards were amplified by a 10-channel linear amplifier, which produced a maximum output modulation of $\pm 7.5$V, for a maximum vibration amplitude of $\pm 0.375$N. All the stimuli were generated and presented using MATLAB and Psychtoolbox 3 [@Kleiner2007; @Brainard1997]. Any sounds produced by the vibrations were rendered inaudible by playing a 440 Hz tone on the remaining two sound card outputs, which was delivered to participants over a pair of headphones.
::: {.content-visible unless-format="docx"}
![Overview of methods details. Panel (a) illustrates the arrangements of baseline (blue circles) and target (red squares) stimuli in the two intervals of a trial in Experiment 1, for four stimulus arrangements (rows). Panel (b) shows an example 'dipper' function. Panel (c) shows the board containing the ten tactors. Panel (d) illustrates the allocation of digits into two sets. Panel (e) illustrates the stimulus conditions used in Experiment 2.](Figures/Figure4.pdf){#fig-methods}
:::
EEG signals were recorded using a 64-channel electrode cap and an ANT Neuroscan (ANT Neuro, Netherlands) amplifier sampling at 1kHz. Electrodes were arranged according to the 10-20 system, and impedances were kept below $5k \Omega$. Digital triggers were sent to the EEG amplifier using a USB TTL module (Black Box Toolkit Ltd., UK), signifying the start of the trial. The whole head average was used as a reference for the EEG data, and the ground electrode was located at position $AFz$.
## Psychophysical procedures
In Experiment 1, a two-interval forced-choice (2IFC) task was used to measure detection and discrimination thresholds. During the experiment, participants were instructed to place their ten digits on the corresponding solenoids, and a series of 26Hz vibration stimuli were delivered to their hands. The digits on each hand were coded 1 to 5 from the thumb to the little finger in that order. The digits 1, 3 and 5 on the left hand and 2 and 4 on the right hand are denoted "Set A", and the digits 2 and 4 on the left hand, and 1, 3 and 5 on the right hand are denoted "Set B" (illustrated in [@fig-methods]d).
Each trial consisted of two 500 ms intervals: one containing the baseline stimulus (equivalent to a pedestal stimulus in studies of visual contrast discrimination) and the other containing the baseline stimulus plus a target increment, separated by a 400 ms interstimulus interval. The next trial began 200 ms after the participant responded. To mask any sound produced by the solenoids and indicate the stimulus intervals, a 440 Hz beep sound was delivered through headphones simultaneously with the vibration stimuli. The order of the two intervals was randomised, and participants were required to determine which interval contained the target increment by pressing a foot pedal. Feedback was provided by a coloured square on the computer screen, with green indicating a correct response and red indicating an incorrect one. The amplitude of the target increment was determined using a pair of 3-down-1-up staircases, with a step size of 3 dB (where dB units are defined as $20 \times log_{10}(100 \times I)$, and $I$ is the stimulus intensity expressed as a proportion of the maximum system output), aiming to distribute trials around the detection threshold at 75\% correct. Threshold measurement was terminated after either 70 trials or 12 reversals, whichever occurred first.
The baseline and target stimuli were presented under four conditions (illustrated in [@fig-methods]a) and at eight baseline intensity levels (0, 1, 2, 4, 8, 16, 32, 64\%), with the baseline intensity level expressed as a percentage of the maximum vibration (0.375N). In the "pentadactyl" condition (meaning five-fingered in Greek), the baseline and target stimuli vibrated only "Set A" (analogous to the stimulation of one eye or one ear). In the "dekadactyl" (ten-fingered) condition, the baseline and target stimuli vibrated both "Set A" and "Set B" (analogous to the stimulation of both eyes or both ears). Comparing the "pentadactyl" and "dekadactyl" conditions reveals the summation effect of doubling the number of inputs. In the "dichodactyl" condition (named for consistency with dichoptic conditions in vision experiments, where different stimuli are shown to the two eyes), the baseline stimulus vibrated "Set A" and the target stimulus vibrated "Set B", permitting the measurement of masking effects across digits. The "half-dekadactyl" condition involved vibrating all ten digits for the baseline stimulus, while only "Set A" vibrated for the target stimulus, enabling measurement of summation (by comparison with the dekadactyl condition) while controlling the number of digits receiving the baseline stimulus [see also the half-binocular condition of Meese et al. @Meese2006]. In all conditions, the vibration frequency for both baseline and target stimuli was 26 Hz. "Set A" and "Set B" were counterbalanced across trials. Each baseline intensity level was tested in a single block lasting approximately 15 minutes and repeated three times. The entire experiment took around 7 hours per participant, resulting in a total of `r ntotaltrials` trials (`r trialspersubj` per participant on average), and was completed over multiple days.
## EEG procedures
After EEG cap set-up, participants in Experiment 2 were exposed to a series of vibrations under six different conditions (illustrated in [@fig-methods]e) at five intensity levels (4, 8, 16, 32, 64\%). Steady-state somatosensory evoked potential (SSSEP) signals were recorded throughout the experiment. During each trial, participants received an 11-s vibration with a 3-s interstimulus interval, and were required to keep their hands still until designated break periods. The "Set A" and "Set B" allocations were the same as those used in the psychophysical experiment (see [@fig-methods]d) and were counterbalanced across trials. The order of conditions was randomised, and each condition was repeated twice for each set, resulting in a total of 120 trials. The entire experiment lasted approximately 30 minutes, split into four 7-minute blocks with rest breaks between blocks.
Two different frequencies were used: 26 Hz (F1) and 23 Hz (F2). F1 is considered the primary frequency, as SSSEP responses are greatest around 26Hz when vibration is delivered to the hands [@Snyder1992]. Using two distinct frequencies enables measurement of both the separate responses to each frequency and suppression between them, without confounding effects from summation. In the "pentadactyl" condition, F1 only vibrated "Set A"; in the "dekadactyl" condition, F1 vibrated both "Set A" and "Set B"; in the "dichodactyl" condition, F1 vibrated "Set A" while a mask at 26 Hz with an intensity level of 32\% vibrated "Set B". These three conditions contributed to the observation of summation and suppression effects at the same frequency.
In the remaining three conditions, a second frequency (F2) was introduced to investigate interactions between different frequencies. In the "cross-pentadactyl" condition, F2 vibrated "Set A", providing a comparison with the "pentadactyl" condition and serving as a baseline for other cross-frequency conditions. In the "cross-dekadactyl" condition, F1 vibrated "Set A" and F2 vibrated "Set B", permitting measurement of suppression effects between different frequencies. In the "cross-dichodactyl" condition, F1 vibrated "Set A" while the F2 mask (intensity level of 32\%) vibrated "Set B", again to measure suppression between digits (see [@fig-methods]e for a diagram of all conditions).
## Data analysis
For both experiments, off-line analysis and statistical testing was performed in Python 3. For psychophysical data, the $psignifit$ 4 package [@Schutt2016] was used to estimate thresholds (at 75\% correct) and slope parameters of the psychometric functions by fitting a cumulative Gaussian, parameterized as:
$$p = \gamma + (1 - \lambda - \gamma)\Phi(C\frac{x-m}{w})$$ {#eq-cumgaus}
where $\lambda$ and $\gamma$ define the upper and lower asymptotes (the lapse rate and baseline accuracy), $x$ is the stimulus intensity (in logarithmic units), $m$ is the threshold parameter, and $w$ determines the slope. (Note that $\gamma$ here is unrelated to $\gamma$ in @eq-binsum2, but we retain this notation for consistency with that used in @Schutt2016.) $\Phi$ denotes the cumulative standard normal distribution, and $C$ is a constant ($C = \Phi^{-1}(0.95) - \Phi^{-1}(0.05) =$ `r qnorm(0.95)-qnorm(0.05)`) (for details, see @Schutt2016). We then converted the slope parameter to equivalent Weibull $\beta$ values using the approximation $\beta = 10.3/\sigma$, where $\sigma = w/C$. We performed this conversion because Weibull $\beta$ is a more commonly used measure of slope, which makes for easier comparison with previous studies, and because steep slopes intuitively correspond to large $\beta$ values. Psychometric function fitting was conducted independently for each participant and condition, and we calculated geometric means across participants of the threshold and slope parameters (implemented as the arithmetic mean of logarithmic values).
For EEG data, all preprocessing was conducted using MNE-Python [@Gramfort2013]. For each trial, the initial 1000 ms post-stimulus presentation was discarded to eliminate onset transients. The remaining 10 s were Fourier transformed and the amplitudes were averaged across repetitions and participants. After removing one outlier, data from 30 participants remained for calculating the intensity-response functions and performing statistical analysis. The amplitudes from six electrodes centred on the area of greatest response ($F1$, $F2$, $Fz$, $FC1$, $FC2$, $FCz$) were then averaged to plot the amplitude spectra and intensity-response functions. The primary dependent variables were the Fourier amplitudes at 26 Hz and 23 Hz.
## Computational modelling
The psychophysical data were fitted using the two stage model of Meese et al [@Meese2006], as outlined in equations 1-4, with 7 free parameters ($m$, $S$, $\omega$, $p$, $q$, $Z$ and $k$). The final parameter ($k$) determines the threshold criteria, such that the model response in the target interval must exceed that in the null interval by $k$ for threshold to be reached. It is proportional to additive noise in the model. We also considered a modified model, in which the summation process described by @eq-binsum involved an additional exponent:
$$sum = (Stage1_L^\gamma + Stage1_R^\gamma)^{1/\gamma}$$ {#eq-binsum2}
The Minkowski exponent ($\gamma$) is implicitly 1 in @eq-binsum, but was free to vary in our second model.
To fit the EEG data, we used a simpler model described by Equations 1-3, but omitting the output nonlinearity described by @eq-stage2, which serves primarily to determine the shape of the dipper function in psychophysical studies, but is not required for EEG data [@Baker2017]. The output was instead defined as:
$$resp = R_{max} \times sum + k$$ {#eq-stage2b}
where $R_{max}$ is a scaling parameter on the output, and $k$ again represents additive noise. The simpler model had 5 free parameters ($m$, $S$, $\omega$, $R_{max}$ and $k$).
We used a downhill simplex algorithm to minimize the squared error between model and data. Fitting for each model was repeated 100 times from random starting vectors, and the parameters that gave the best numerical fit were selected. For the psychophysical data, we first fitted the model to the averaged threshold data by determining the target contrast for each condition that was required to increase the model response by the criterion ($k$). At the suggestion of a reviewer, we additionally performed maximum likelihood fits to the raw proportion correct data. These fits are summarised in the Supporting Information Appendix, and lead to similar conclusions to the modelling reported in the body of the manuscript.
## Data and code availability
All experiment code, raw data, processed data and analysis code are available at the project repository: https://doi.org/10.17605/OSF.IO/M79D2. The linked GitHub repository also contains a fully computationally reproducible version of this manuscript in Quarto format.
## Acknowledgements
We thank Mark Green for constructing the vibration board, and Kirralise Hansford for her support during EEG data collection. This work was supported by the China Scholarship Council [202308510049], and by BBSRC grant BB/V007580/1 awarded to DHB and ARW.The funder had no role in the design of the study, data collection, analysis, or decision to publish.
# References

TestVersion.qmd, no license · at the source

Overview

Authors: Shasha Wei¹, Alex R Wade^1,2, Catherine EJ Preston¹, Daniel H Baker¹

ORCID iDs: Shasha Wei, Alex R Wade, Daniel H Baker

Department of Psychology, University of York, York, United Kingdom
York Biomedical Research Institute, University of York, York, United Kingdom

Institutions: University of York (United Kingdom)

Journal: PloS one, volume 21, issue 6, article e0350140

Dates: received 17 November 2025; accepted 8 May 2026; published online 11 June 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1371/journal.pone.0350140 · PMID 42275386 · PMCID PMC13258021 · OpenAlex W7164346268

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: EEG (modality), human (organism), cognitive (subfield)

MeSH: Evoked Potentials, Somatosensory*, Touch Perception*, Vibration*, Adult, Electroencephalography, Female, Fingers, Humans, Male, Sensory Thresholds, Young Adult (* major topic)

Topic: Multisensory perception and integration (Experimental and Cognitive Psychology, Psychology), according to OpenAlex

Funding: Biotechnology and Biological Sciences Research Council (BB/V007580/1)

Citations: not cited yet (Europe PMC); 65 references in the paper

Abstract

While the brain’s integration of auditory and visual inputs has been extensively investigated, the mechanisms underlying somatosensory signal combination remain less explored. Here, we combine psychophysical thresholds with steady-state somatosensory evoked potentials (SSSEPs) to investigate how vibrotactile inputs are combined across fingers. We find that doubling the number of stimulated digits leads to a weak improvement in detection threshold, consistent with probability summation, whereas introducing a masking stimulus to interleaved digits induces inter-digit suppression. Correspondingly, EEG recordings reveal a ~ 1.4-fold increase in SSSEP amplitude when doubling the number of digits stimulated at the same frequency, reflecting a summation effect. In contrast, SSSEP amplitudes decrease when digits are vibrated at two different frequencies, further supporting the presence of suppression. These results are consistent with a model featuring inhibition between digits and reveal that the weight of suppression is intermediate between that observed in binocular vision and binaural hearing.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 6 matches between paragraphs and lines of code.

OSF m79d2

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Languages: Quarto (3), MATLAB (2)

Size: 50 files, 5 scripts

Software Heritage: not checked

Found in: “Data Availability”

Holds: environment (docker-compose.yml, docker/Dockerfile), continuous integration, 3 notebooks

Not found: README, license file, CITATION.cff, tests, documentation

Tools: Matplotlib (2 files), MNE-Python (2 files), NumPy (2 files), pandas (2 files), Pingouin (2 files), psignifit (2 files), Psychtoolbox (2 files), SciPy (2 files), statsmodels (2 files)

Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)

27 September 2026: the link answers (HTTP 200)

5 files

Experiment code/tactileEEG.m, MATLAB, 323 lines, 1 match
Experiment code/tactiledippers.m, MATLAB, 440 lines, 1 match
SupportingInfo.qmd, Quarto, 661 lines
TestVersion.qmd, Quarto, 1,734 lines, 2 matches
vibrosummanuscript.qmd, Quarto, 1,578 lines, 2 matches

At the source: doi.org/10.17605/OSF.IO/M79D2

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
5 scripts, each with its path and the digest of its content;
6 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

All experiment code, raw data, processed data and analysis code are available at the project repository: https://doi.org/10.17605/OSF.IO/M79D2. The linked GitHub repository also contains a fully computationally reproducible version of this paper in Quarto format.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 11 MeSH terms, 1 funder, 64 references.

Cite

This paper

Wei, S., Wade, A. R., Preston, C. E., & Baker, D. H. (2026). Signal combination in flutter vibration perception. PloS one, 21(6), e0350140. https://doi.org/10.1371/journal.pone.0350140

BibTeX

@article{wei2026signal, author = {Wei, Shasha and Wade, Alex R and Preston, Catherine EJ and Baker, Daniel H}, title = {{Signal combination in flutter vibration perception}}, journal = {PloS one}, year = {2026}, month = jun, volume = {21}, number = {6}, pages = {e0350140}, publisher = {PLOS}, issn = {1932-6203}, doi = {10.1371/journal.pone.0350140}, url = {https://doi.org/10.1371/journal.pone.0350140}, pmid = {42275386}, pmcid = {PMC13258021} }

RIS

TY - JOUR AU - Wei, Shasha AU - Wade, Alex R AU - Preston, Catherine EJ AU - Baker, Daniel H TI - Signal combination in flutter vibration perception T2 - PloS one J2 - PLoS One PY - 2026 DA - 2026/06/11 VL - 21 IS - 6 SP - e0350140 SN - 1932-6203 PB - PLOS DO - 10.1371/journal.pone.0350140 UR - https://doi.org/10.1371/journal.pone.0350140 LA - en ER -

CSL-JSON

{ "id": "10.1371/journal.pone.0350140", "type": "article-journal", "title": "Signal combination in flutter vibration perception", "container-title": "PloS one", "author": [ { "family": "Wei", "given": "Shasha" }, { "family": "Wade", "given": "Alex R" }, { "family": "Preston", "given": "Catherine EJ" }, { "family": "Baker", "given": "Daniel H" } ], "container-title-short": "PLoS One", "volume": "21", "issue": "6", "page": "e0350140", "DOI": "10.1371/journal.pone.0350140", "PMID": "42275386", "PMCID": "PMC13258021", "ISSN": "1932-6203", "publisher": "PLOS", "URL": "https://doi.org/10.1371/journal.pone.0350140", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 11 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.