Peaglet: A user-friendly probabilistic Kernel density estimation of intracranial cortical and subcortical stimulation sites

Background: Data on human brain function obtained with direct electrical stimulation (DES) in neurosurgical patients have been recently integrated and combined with modern neuroimaging techniques, allowing a connectome-based approach fed by intraoperative DES data. Within this framework is crucial to develop reliable methods for spatial localization of DES-derived information to be integrated within the neuroimaging workflow. New method: To this aim, we applied the Kernel Density Estimation for modelling the distribution of DES sites from different patients into the MNI space. The algorithm has been embedded in a MATLAB-based User Interface, Peaglet . It allows an accurate probabilistic weighted and unweighted estimation of DES sites location both at cortical level, by using shortest path calculation along the brain 3D geometric topology


Introduction
Building on the pioneering works of Penfield and colleagues, direct electrical stimulation (DES) in patients undergoing neurosurgery is widely recognized as a crucial tool for mapping human brain functions.Providing a functional description of brain regions is considered a longstanding goal in neuroscience, potentially offering crucial information necessary for developing new treatments for neurological disorders.However, DES alone does not have the potential to reveal the functional and structural patterns of the stimulated brain sector, hindering a network-based description.In this regard, the implementation of noninvasive functional neuroimaging techniques in the late twentieth century greatly extended the brain mapping potential, enabling whole brain connectomic analysis.As a side effect of this technological transition, functional neuroimaging can only provide correlational or indirect evidence that, compared to DES (or lesion-based approach), does not allow to investigate the causal link between a given behaviour or symptom and cortical-subcortical neuroanatomy (Siddiqi et al., 2022).This aspect is crucial for translating knowledge into therapeutic targets.
To overcome the intrinsic individual limitations of these different techniques, in the last decade the mapping of brain function based on DES in neurosurgical patients has been combined with modern neuroimaging techniques.The combination of DES-related probabilistic atlases within the neuroimaging environments carries the potential to highlight the functional and structural connectome of regions causally involved in specific brain functions, thus improving the knowledge about the brain architecture behind specific DES-related effects (Simone et al., 2021(Simone et al., , 2020;;Fornia et al., 2020b;Elmalem et al., 2023;Coletta et al., 2023;Rech et al., 2019).Moreover, the application of DES in patients undergoing awake neurosurgery for brain tumour resection further enriched this context by providing both cortical and subcortical white matter information regarding the localization of specific brain function (Viganò et al. 2022;Puglisi et al., 2019;Sarubbo et al., 2020), allowing a highly comprehensive mapping finally impacting on survival and quality of life of neuro-oncological patients (Rossi et al., 2021;Leonetti et al., 2021).However, the integration between DES results with neuroimaging techniques is still in its infancy.While functional and structural neuroimaging are subserved by consolidated and continuously evolving methods and available tools for spatial estimation of brain activity and dissection of white matter pathways, the spatial estimation of DESrelated information does not benefit from a comparable working framework.In this regard, to improve the specificity of the results obtained by combining modern neuroimaging techniques and intracranial stimulation, robust tools for a reliable spatial estimation of DES-related information are required.Although several groups have applied spatial estimation methods to quantify the localization of DES-derived information (Tate et al., 2014;Rech et al., 2019;Sarubbo et al., 2020;Elmalem et al., 2023), there are no available tools and common methodological approaches.
To overcome this limitation, in the present study we introduced a MATLAB-based tool, Peaglet, designed for comprehensive kernel density estimation of intraoperative stimulation sites in the MNI space.As a novelty, Peaglet probabilistic estimation can be applied on both cortical and subcortical DES sites taking into account the distinguishing intrinsic geometric features of the two tissues.Overall, the aim of Peaglet is to provide a robust probabilistic estimation of the cortical and subcortical distribution of DES sites going beyond a simple region of interest approach, allowing a user-friendly visualization and implementation of the results within the neuroimaging workflow to drive precise DESinformed connectomic analysis.

Materials and methods
Peaglet is a Matlab package of functions that displays a graphical user interface.It has been developed on MATLAB version 2019A and it requires Tools for NIfTI and ANALYZE image (https://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze -image).Peaglet is accessible through the following link https://github.com/mocalab-crew/peagletpub.The Matlab environment has been chosen since the availability of well-established and integrable MRI image processing libraries.

Cortical analysis
The cortical analysis is applied to DES cortical sites and it is based on three methodological aspects: 1) a surface representation of the human brain in the MNI space; 2) Kernel Density calculation to estimate the spatial probability of a group of MNI coordinates corresponding to each stimulation site; 3) Shortest Path calculation to quantify the distance among stimulated sites.

Surface representation
The MNI templates available are the ICBM 152 (306716 vertex) and the FSAverage (327684 vertex) extracted from Brainstorm (Tadel et al., 2011).The MNI templates are reconstructed as a network of points to faithfully capture its inherent three-dimensional topology.This representation is pivotal for maintaining the accuracy of subsequent analyses, such as kernel density estimation, as it avoids oversimplification and preserves the realistic topological features of the brain.

Kernel Density Estimation (KDE)
To facilitate KDE computation, a grid covering the spatial domain of interest is generated.For each point on this grid, the KDE is calculated by adding contributions from all data points.The kernel function, appropriately scaled by the specified bandwidth, ensures accurate representation of spatial relationships among DES sites.

Shortest path integration
Shortest paths between surface points are computed as an integral step in our methodology, aiming to preserve the inherent threedimensional topology of the brain throughout the KDE computation.This process involves determining the shortest routes between pairs of surface points on the reconstructed brain network.By doing so, we ensure that the spatial relationships between these points are accurately quantified, maintaining the structural integrity and complexity of the 3D brain topology in the KDE calculation.

Subcortical analysis
The subcortical analysis is based on brain volume and it has been developed for analysing stimulation coordinates obtained within the white matter.This is particularly important for spatial estimation of DES results obtained during subcortical disconnection of the brain tumour mass.In fact, as suggested by papers published in the last decade, subcortical networks play a critical role in brain functioning, partially in opposition with the previous concept of localizationism.Compared to cortical, for subcortical analyses a similar procedure has been applied, with the main difference in the calculation of probability density estimation through KDE.In this case, users can customise the volume of interest based on visual inspection of data distribution.

Peaglet utilisation
After setting the path, Peaglet can be launched by typing in the command window: PDE2.
Step 1. Load Data.Users start by loading a set of x,y,z MNI coordinates corresponding to stimulation sites provided as txt file.Each coordinate can be associated to a value (w) for weighted probability density estimation.If no "w" value is reported, non-weighted probability density estimation will be automatically performed (Figs.1-2).
Step 2. Brain Model Selection.The tool allows users to choose between different brain models depending on data pre-processing and normalisation procedure.This selection determines the subsequent analysis approach (Figs.1-2).
Step 3. Analysis.Cortical Analysis.In order to perform the cortical analysis, the tool automatically determines the neighbouring points in the brain network for each stimulation site based on user input data.This adaptive approach optimises sampling and computational efficiency (Fig. 1).Subcortical analysis.For subcortical analyses, users can customise the volume of interest and specify the sampling based on visual inspection of data distribution (Fig. 2).The tool provides a visual representation of the chosen sampling strategy (number of neighbours) on the screen, allowing users to inspect their selections before proceeding with the analysis.
Step 4. Parameters.Bandwidth Selection: users can flexibly choose the bandwidth parameter for the kernel density estimation, tailoring the smoothness of the estimated density according to specific criteria (i.e.amount of subcortical or cortical tissue possibly stimulated by DES).
Step 4.1.Weighting Options: Additionally, users can weigh each site with different strategies: linear, exponential or by using directly a multiplier, providing finer control over the contribution of individual A. Bellacicca et al. points to the estimation.The weighting strategy employed by the function is based on an oversampling approach aimed to shape the KDE calculation.Essentially, a coordinate with a higher weight is represented by a proportionally or exponentially greater number of "virtual points" at the same spatial location, finally impacting on density estimation.
Step 5. Estimate PD.When all parameters are fixed, the tool computes the kernel density estimation based on the selected parameters, producing a spatial probability density representation of the stimulated sites.The resulting estimation is displayed alongside the original sampling image for comparative analysis.
Step 6. Resample PDE.An important step of our analytical framework is the resampling procedure applied to Probability Density Estimates (PDEs).This process is crucial for adapting the granularity and range of PDEs.The function resample_pde_Callback presented herein encapsulates the resampling methodology, accommodating user-defined lower and upper limits (respectively "Lower Limit TextBox and Upper Limit TextBox").
Step 8. Save results.The user can save the resulting PDE to a NIfTI format.First, the user selects a NIfTI template file (1x1x1 "stamp.ni"file) through a user-friendly file selection dialog.The stimulation points and resulting probability density estimation, are then retrieved from the graphical interface and included in the NIfTI data.For visualisation purposes, the user can enhance the stimulation data using the "Voxel Filler" to increase the area shown.

Results
We provided some examples of Peaglet functions by analysing a set of cortical and subcortical DES coordinates extracted from seminal papers in the field of brain tumour resection (Sarubbo et al., 2019;Sarubbo et al., 2020).More specifically, we estimated the cortical localization of speech arrest sites and the subcortical localization of semantic language errors.Since these coordinates were not associated with any weighted value, to provide an example of the weighted probability density function we added for each coordinate a random value from 2 to 5 representing the usual DES intensity threshold range in mA required to obtain the desired effect.

Cortical parameters and results
This analysis was performed by using the following parameters: • Template.Since each coordinate was reconstructed and normalised following the ICBM 152 template, the same template was employed in further analysis.

Fig. 1.
Example of cortical analysis.1) STEP 1: load the txt file (x,y,z,w); 2) STEP 2: Template selection from the drop down menu; 3) STEP 3: select cortical analysis and play with neighbours to fit data distribution.Select "show sampling" to visualise the neighbours selection; 4) STEP 4: apply bandwidth; STEP 4.1: apply the weighting function if required; 5) STEP 5: select Estimate PD to run the analysis; 6) STEP 6: select the proper limits-voxel fillers: 7) STEP 7: select STEP 7 "resample PDE"; finally 8) STEP 8: save results as NIfTI file.Select Save PDE and then select stamp.niifile.Below three distinct outputs are shown: 9) unweighted, 10) weighted exponential and 11) linear.12) The file can be loaded into the Anatomy Toolbox for exploring their anatomical features.13) Furthermore, the files can be thresholded and saved as binary files to explore their functional connectivity with CONN.
• Analysis.A neighbour of 10 was chosen based on visual inspection of stimulation site distribution.• Parameters.A bandwidth of 8 was chosen.Ideally this value should be set based on the type of data.In our examples the DES was performed with a bipolar stimulation probe with a 5 mm distance between tips.In our experience a bandwidth parameter of 6-8 likely reflects the amount of tissue stimulated.However, we recommend setting this parameter following a trade off between the stimulated tissue and DES sites distribution.• Estimate PDE.After probability density estimation of DES sites, the results were resample using lower-upper limits "0-100" and a "0" Voxel filler (the latter can be increased in order to enlarge the results for a better visualisation).• Save NIfTI and Export Stim points.The results were saved by selecting save nifti and choosing stamp.niifile.Both unweighted and weighted probability density estimation are saved (the latter is labelled with "_w" as postfix).Single stimulation sites may also be exported as a NIfTI file.
As example, the results are easily visualised with SurfIce (https ://www.nitrc.org/projects/surfice/)and anatomically investigated with anatomy toolbox (https://github.com/inm7/jubrain-anatomy-toolbox).Additionally, the NIfTI file can be thresholded and used as DES-derived functional ROI to explore its pattern of functional connectivity.See Fig. S1A to visually inspect the functional connectivity results.

Sub-cortical parameters and results
This analysis was performed using the following parameters: • Template.Since each stimulation coordinate was reconstructed and normalised following the ICBM 152 template, the latter has been chosen for further analysis.• Analysis.After visual inspection of the site's distribution the analysis was performed by using a volume multiplier of 1.2 and a Sampling of 50.• Parameters.Compared to cortical analysis, a more conservative bandwidth of 10 was chosen considering the stimulation sites distribution and the different spread of current at subcortical level.• Estimate PD.After probability density estimation of DES sites, the results were resample using lower-upper limits "0-100" and "1" Voxel filler.• Save NIfTI and Export Stim points.The results were saved by selecting save nifti and choosing stamp.niifile.Both unweighted and weighted probability density estimation are saved (the latter is labeled with "_w" as postfix).
The results are easily visualised with MRCronGL, and compared with the atlas to investigate the white matter pathways mainly involved (https://www.nitrc.org/projects/mricrogl).Alternatively, the NIfTI file can be thresholded and used as DES-derived functional ROI to explore its disconnection map with BCBToolkit (Foulon et al., 2018).See Fig. S1B to Fig. 2. Example of subcortical analysis.For all the steps see Fig. 1.Since this analysis is related to subcortical white matter stimulation sites, the results can be explored with visualisation software such as MRIcroGL, or the maps can be thresholded and used as binary file for computing disconnectome analysis with BCBtoolkit.
visually inspect the disconnectome results.

Comparison with existing methods
The main distinction between Peaglet and other methods focused on spatial distribution of DES-related data (Sarubbo et al., 2020;Tate et al., 2014;Rech et al., 2019;Elmalem et al., 2023) lies firstly in the differentiation between cortical and subcortical analyses.In cortical analysis, point-to-point distances were computed considering the actual surface of the brain template, as opposed to a volumetric linear distance.This approach aims to faithfully reflect the complexity of cortical topology in the calculation of the distance among DES sites.Concerning subcortical analyses, a volumetric approach is applied, allowing users to flexibly customise parameters for volumetric analysis based on the nuances of their dataset.This flexible strategy allows for a better convergence between the specific requirements of each user and the variable distribution and numerosity of DES-derived information, which may also depend on the number of patients investigated.Furthermore, Peaglet can perform both weighted and unweighted probability density estimation based on an oversampling strategy allowing a fine-grained four-dimensional description (Avanzini et al., 2016) of DES-derived information.This function shapes the density calculation based on quantified parameters instead of considering the DES coordinates as an homogeneous dataset associated with a specific effect.
As a potential limitation, Peaglet does not apply clustering algorithm or statistical parametric analysis to refine the results (Tate et al., 2014;Elmalem et al., 2023).The implementation of these additional methodological aspects will be considered in the future.

Discussion
The advent of a multimodal brain mapping approach which combines and integrates modern neuroimaging with intracranial stimulation techniques allowed to extend DES-derived information at connectome level.In this regard, the development of systematic and consistent methods for spatial estimation of DES-derived information going beyond a region of interest approach (or average of coordinates) and easily integrated within the neuroimaging formats, is important to increase the anatomo-functional specificity of connectome results.With the aim of improving this new era of multimodal neuroimaging approaches and laying the foundation for spatial estimation of DES-related information, we have developed Peaglet.
The present MATLAB tool stands as a valuable contribution to estimate the probabilistic spatial localization of intracranial DES sites.By integrating advanced functionalities such as brain model selection, sampling, and customizable analysis parameters, the tool addresses the potential challenges posed by different, highly variable and clinically constrained, DES data in MNI coordinates.Researchers will benefit from a user-friendly environment that not only streamlines data processing, but also allows tailored analyses in both cortical and subcortical contexts fitting the geometrical features of the two neural substrates.
The outputs of Peaglet are NIfTI maps representing probabilistic estimation of a list of DES sites in MNI coordinates, possibly grouped based on the specific outcome evoked by stimulation (i.e.speech arrest or semantic errors).These maps can be used 1) to investigate their anatomical localization with common neuroimaging software and 2) as DES-derived functional regions of interest to inform structural and functional connectome analysis.While the former allows a better investigation of the anatomical counterpart of a specific brain function derived by DES, the latter allows investigating the networks subserving this function.In this regard, the greater causality and spatial resolution of intraoperative DES raises the question about how to integrate this benefit into the neuroimaging environment.The use of Peaglet, which estimates probabilistic cortical and subcortical DES-derived information, allows to increase the specificity of the functional and/or structural connectivity of DES-functionally defined regions which cannot be geometrically approximate to a single sphere representing average values among coordinates.This feature becomes crucial in the evaluation of the functional connectome that depends on the average spontaneous BOLD signal within the voxels composing the DES-derived seed.
Regarding limitations, as clearly acknowledged in the part dedicated to "comparison with other methods", Peaglet does not re-shape probability estimation of DES-derived information based on further statistical methods, such as clustering algorithm or statistical parametric analysis (Tate et al., 2014;Elmalem et al., 2023).These implementations will be embedded in a future release with other potential limitations for the users, including visual feedback indicating the status of the analysis, displaying the name of the loaded data file and an output file recording the parameters selected for the current analysis.

Conclusion
Peaglet has been developed and implemented in the field of DES for brain tumour resection.In this clinical context the key concept is the "onco-functional balance" i.e. to resect as much as tumour as possible while preserving the functional integrity of the patient.The DES allows surgeons to preserve eloquent brain structures and its efficacy in reaching an "onco-functional balance" has been demonstrated (De Wit et al., 2012).Furthermore, extensive application of intraoperative brain mapping strategies led to a switch from a surgery guided by the images to a surgery guided by the function, allowing to reach a supratotal resection with a greater impact on the oncological outcome (Rossi et al., 2021).In this framework, the methods implemented within Peaglet has been recently used to investigate DES-sites evoking specific motor and cognitive impairments (Fornia et al., 2020a;Fornia et al. 2020b;Fornia et al., 2022;Viganò et al. 2022;Fornia et al., 2024;Puglisi et al., 2019) providing evidence about the spatial localization of important eloquent functional borders to be preserved during the surgery.However, although implemented in this clinical field, Peaglet can be adapted to investigate a set of MNI coordinates regardless of their specific origins, extending its potentiality to spatial estimation of others intracranial information, such as electrophysiological responses recorded by stereo-EEG and ECoG electrodes.