# Deep learning-based brain segmentation model performance validation with
clinical radiotherapy CT

Selena Huisman<EMAIL_ADDRESS>Matteo Maspero  Marielle
Philippens  Joost Verhoeff  Szabolcs David  Department of Radiation Oncology,
Amsterdam UMC, De Boelelaan 1117, 1081 HV Amsterdam, The Netherlands
Department of Radiation Oncology, UMC Utrecht, Heidelberglaan 100, 3508 GA
Utrecht, The Netherlands.  Computational Imaging Group for MR Diagnostics $\&$
Therapy, UMC Utrecht, Heidelberglaan 100, 3508 GA Utrecht, The Netherlands.

###### Abstract

Introduction Manual segmentation of medical images is labor intensive and
especially challenging for images with poor contrast or resolution. The
presence of disease exacerbates this further, increasing the need for an
automated solution. To this extent, SynthSeg is a robust deep learning model
designed for automatic brain segmentation across various contrasts and
resolutions. This study validates the SynthSeg robust brain segmentation model
on computed tomography (CT), using a multi-center dataset.
Methods An open access dataset of 260 paired CT and magnetic resonance imaging
(MRI) from radiotherapy patients treated in 5 centers was collected. Brain
segmentations from CT and MRI were obtained with SynthSeg model, a component
of the Freesurfer imaging suite. These segmentations were compared and
evaluated using Dice scores and Hausdorff 95 distance (HD95), treating MRI-
based segmentations as the ground truth. Brain regions that failed to meet
performance criteria were excluded based on automated quality control (QC)
scores.
Results Dice scores indicate a median overlap of 0.76 (IQR: 0.65-0.83). The
mean volume difference is 7.79% (CI: 6.41%–9.18%), with CT segmentations
typically smaller than MRI-based. The median HD95 is 2.95 mm (IQR: 1.73-5.39).
QC score based thresholding improves median dice by 0.1 and median HD95 by
0.05mm. Morphological differences related to sex and age, as detected by MRI,
were also replicated with CT, with an approximate 17% difference between the
CT and MRI results for sex and 10% difference between the results for age.
Conclusion SynthSeg can be utilized for CT-based automatic brain segmentation,
but only in applications where precision is not essential. CT performance is
lower than MRI based on the integrated QC scores, but low-quality
segmentations can be excluded with QC-based thresholding. Additionally,
performing CT-based neuroanatomical studies is encouraged, as the results show
correlations in sex- and age-based analyses similar to those found with MRI.

###### keywords:

Deep learning, Validation, Clinical Brain MRI, Clinical Brain CT

##  1 Introduction

Segmentation of the human brain can be used to measure regional volumetric
differences between healthy subjects and patients, aiding in visualizing and
quantifying brain structures [1]. Manual segmentation requires extensive labor
while providing variable performance between observers, with median Dice
scores falling between 0.46 and 0.89 and median HD95 between 2.8 and 3.3
depending on the brain region [2]. Consequently, automatic and reproducible
segmentation models have been developed, which usually require a T1-weighted
magnetic resonance imaging (MRI). However, when images are only available from
computed tomography (CT), the task of tissue segmentation becomes challenging
due to the inferior soft tissue contrast of CT compared to MRI. To address
shortcomings, contrast-agnostic tools have been trained on artificially
downsampled MRIs. Such models enable automatic and reproducible segmentation
of brain images regardless of contrast and resolution, effectively enabling
CT-based soft-tissue segmentation. Recently, SynthSeg has been proposed [3] as
a promising contrast and resolution agnostic, whole-brain segmentation model.
This model, trained through exclusively synthetic images with randomized
contrasts and resolutions, exemplifies this robust approach, offering
satisfactory CT segmentation performance [4]. Direct comparative validation of
SynthSeg and similar models is challenging, due to the need for paired CT and
MRI, in which image pair availability is limited in both healthy volunteers
and patients. In radiotherapy (RT), however, patients with brain tumors
routinely undergo both CT and MRI for radiation planning, making this setting
ideal for validation. Usually, a large part of the patient’s brain is not
affected by tumors, allowing us to compare the performance of CT and MRI-based
segmentation of such regions. Additionally, SynthSeg incorporates automatic
quality control (QC) scores to assess the accuracy of segmented brain regions,
enabling the exclusion of diseased areas from the analysis.

We have assembled a benchmark dataset from three open sources, comprising 260
paired oncological brain CT and MRI images sourced from a total of five
medical centers. This diverse dataset offers a comprehensive resource for
validating the SynthSeg model. While quantifying CT and MRI segmentation
performance is valuable, tissue segmentations are ultimately tools in
exploring various medical or scientific questions, for example the analysis of
volumetric differences based on sex or age, which have been extensively
explored in the human brain. [5, 6, 7] .

In this work, we compared the segmentation performance of SynthSeg on CT and
MRI of patients with brain tumors. Additionally, we examined the volumes of
various brain regions between males and females using both imaging techniques.
Furthermore, we assessed how the brain volumes of patients change with age and
investigated whether these results obtained from MRI segmentation are
reproducible using CT-based segmentation. Although our study does not
primarily investigate these well-documented differences, we employ both MRI
and CT-based segmentation to replicate such findings conceptually. This
approach allows us to test the viability of CT-based segmentation in
neuroscientific studies, which have predominantly used MRI. This comparison
aims to expand the methodologies available for such research and validate the
effectiveness of CT in capturing relevant brain differences, when MRI is not
available.

##  2 Methods

Data was sourced from three openly accessible datasets, containing patients
with primary or metastatic brain tumor for a total of 260 paired CT and MRI.
The first dataset was the SynthRAD2023 Grand Challenge
(https://SynthRAD2023.grand-challenge.org/) training dataset111Retrieved on
08-09-2023 from https://zenodo.org/doi/10.5281/zenodo.7260704., comprising a
subset of 180 brain CT and T1w 3D gradient echo MRI and are available for
public access [8]. These images were collected from three Dutch hospitals: UMC
Utrecht, UMC Groningen, and Radboud Nijmegen. Furthermore, we utilized similar
CT-T1w paired subsets from the GLIS-RT [9] and Burdenko-GBM-Progression
projects [10], from the TCIA database [11], which are also available openly.
The GLIS-RT data was collected at Massachusetts General Hospital, while the
Burdenko data is from the Burdenko National Medical Research Center of
Neurosurgery. All CT and T1 images have been registered and rescaled to 1mm
isotropic resolution, no additional processing was applied to the data. An
overview of the datasets and demographics can be found in table 1. GLIS-RT
provides contrast-enhanced 3D-T1 weighed MRI and radiotherapy planning CT.
Meanwhile, Burdenko-GBM-Progression includes T1-weighed MRI and topometric CT,
in which the MRIs were obtained from four vendors with varying scanning
protocols and the CTs were obtained with a single scanning protocol. In
principle, our analysis could be replicated by many radiotherapy departments
with access to similar data, the open nature of our sources and the mixture of
contributing centers enhancing the transparency and reproducibility of our
results.

Next, we used the SynthSeg tool within the Freesurfer software suite, version
7.4.1 [4] to perform automatic brain segmentation based on 95 brain regions.
The model was run in robust mode, with the integrated QC score and volume
options enabled. The QC scores evaluate the reliability of each segmentation
by comparing the segmented output against a learned model of what constitutes
a high-quality segmentation. These scores are generated by assessing features
of the segmentation against the predicted values from the ’regressor R’ of the
SynthSeg model, providing a numerical indication of segmentation quality. For
CTs, the corresponding option was enabled to facilitate CT-specific
processing. The QC scores are divided into eight categories: general white
matter, general grey matter, general cerebrospinal fluid (csf), cerebellum,
brainstem, thalamus, putamen & pallidum, and hippocampus & amygdala. The
resulting data was examined using Rstudio version 2023.9.1.494 [12], applying
a manual region-specific threshold for QC score exclusion instead of the
static 0.65 threshold suggested by the SynthSeg manuscript [4]. This manual
thresholding approach was chosen to account for variations in QC score
averages across different brain regions. In some cases the 0.65 threshold is
too lenient resulting in unacceptable segmentations being included. Our manual
thresholds correspond with the range of 0.55-0.75 found in the SynthSeg
manuscript. Brain regions with inadequate QC scores were selectively removed
from each image, preserving the remaining regions for subsequent analysis. The
analysis included calculating the percentage mean volume difference for each
brain region, visualized using bar plots and percentage-based Bland-Altmann
plots, also known as Giavarina plots [13]. Dice scores and Hausdorff 95
distances (HD95) for each brain region were computed using the
segmentation_metrics python package, based on the segmented outputs [14]. Non-
overlapping regions were excluded from the analyses as they resulted in 0 Dice
score and physically implausible HD95. Also, when a label was missing for a
patient from either the CT or the MRI, that particular label was excluded from
analyses. Note, that only CT-based segmentations resulted in any missing
regions. Upon visual inspection, all missing labels in the CT or the non-
overlapping regions only occurred in the tumor affected areas, for which
regions SynthSeg is not meant to produce brain segmentation labels. The
segmentations were visually inspected using FSLeyes from FSL version 6.0 [15].
For visualization purposes, the brainstem, left thalamus, left hippocampus and
left cerebellum white matter were highlighted in the figures.

Table 1: An overview of the data used.

A sex and age-based analysis was conducted using the 60 patients from UMC
Utrecht. These analysis considered eight regions namely the left-thalamus,
-insula cortex, -superior frontal cortex, -hippocampus, -putamen, -pallidum,
-lateral ventricle, and the brainstem. The demographic details for these
patients are not publicly available due to privacy regulations and were
retrieved in accordance with the institutional review boards (nWMO research
number 22-474). Figure 1 provides a full overview of the analysis pipeline.

Figure 1: Pipeline of data processing and analyses. From left to right: The
patient’s CT and MRI are processed by SynthSeg, which outputs QC scores and
segmented volumes for both image modalities. The resulting segmentations are
then filtered by a QC-based threshold, after which statistical analyses is
performed.

##  3 Results

To visualize the difference between a CT that received a low QC versus a high
QC score, figure 2 and 3 showcase CT and MRI segmentations in two example
patients. Figure 2 shows a patient with a low mean CT QC score of 0.364 and a
high MRI QC score of 0.813. While the segmentation in both image modalities
fails around the tumor and treatment affected area of the brain, the CT
segmentation also fails to capture accurately non-affected areas, for example
the brainstem, the cerebellum and even the contralateral cortical areas
relative to the tumor cavity.

Figure 2: Segmentation of patient 1BB177 of the SynthSeg cohort. The MRI and
CT segmentations are shown at the top and bottom of the image, respectively.

Figure 3 shows an example segmentation with high QC scores for both MRI and
CT, with a mean of 0.851 and 0.839, respectively. Both segmentations capture
nearly all normal-appearing parts of the brain.

Figure 3: An example segmentation of patient 1BA014. The MRI and CT
segmentations are shown at the top and bottom of the image, respectively.

MRI QC scores show low variability (SD = 0.025, 0.019 excluding CSF), for CSF,
which exhibits a higher SD of 0.072, while CT QC scores show high variability
(SD = 0.139). Figure 4 shows the comparison of CT-based and MRI-based QC
scores from all 260 patients. Brain regions with CT QC scores below the
corresponding threshold were subsequently excluded from the analysis, which
results in a SD of 0.036. Most data points lie above the unity line (in blue)
in the figures, meaning that MRI-based segmentations yield higher QC scores
than those from CT-based segmentations. In all categories, MRI QC scores are
significantly higher than CT QC scores. Furthermore, no brain regions were
excluded based on low MRI-based QC scores if they were not already excluded
due to low CT QC scores. A detailed breakdown of filtered and unfiltered mean
QC scores for each brain region is available in the supplementary materials,
table 1. The remaining four QC score plots can be found in supplementary
materials, figure 1.

Figure 4: Four plots showing the CT QC scores against the MRI QC scores per
patient. The x-axis shows the CT QC scores, while the y-axis shows the MRI QC
scores. The blue line has a slope of one to indicate where the scores would be
equal. The red line is the threshold selected for the exclusion of
segmentations.

The distribution of volumes per brain region is shown in figure 5. On average,
the CT volumes are significantly smaller than the MRI volumes, with an average
of 7.79% (CI: 6.41%–9.18%) absolute difference, although the shape of the
distributions is similar. 73 regions have significantly lower volume for CT
compared to MRI, 16 show no statistical difference, while 6 regions
(parahippocampal cortex in both hemispheres, both thalami; posteriorcingulate
and paracentral in the left cortex) have a significantly larger volume for CT
compared to MRI.

Detailed data on each region’s numerical percentage differences and total
volumes are provided in supplementary table 2. If the confidence interval of
the difference between regions includes 0, it indicates no significant
difference. A visualization of these percentage differences can be found in
supplementary figure 1.

Figure 5: Four violin plots show the volume distribution of between the MRI
and CT segmentation volumes for four brain regions. On the x-axis the volume
in mm3 is found, while the imaging type is listed on the y-axis. CT is colored
red, while MRI is colored blue.

To illustrate the bias in the distribution of the volume differences between
image modalities, figure 6 shows the percentage differences of volumes plotted
against the mean volume for four of the 95 brain regions. The distribution
tends to vary more for regions with smaller mean volumes. Additionally, CT
segmentation volumes are generally smaller than those from MRI, except in the
thalamus and ventricles.

Figure 6: Four Giavarina plots for four different brain regions. The x-axis
has the mean volume in mm3, while the y-axis shows the percentage difference
for each segmentation. The bias and standard deviations are shown within each
plot.

A median Dice score of 0.76 (IQR: 0.65-0.83) was found from all regions,
calculated on the overlap of each region between MRI and CT. Before QC
thresholding the median overall Dice score was 0.66 (IQR: 0.50-0.78). Figure
7/A shows the distribution of Dice scores for four brain regions. White matter
regions generally have higher Dice scores than regions of grey matter. Figure
7/B shows the HD95 distances with a median HD95 distance of 2.95 mm (IQR:
1.73-6.16) across all regions. Before QC thresholding the median HD95 was 3.
The Dice scores and HD95 distances for all regions are available in
supplementary table 2.

Figure 7: Two sets of boxplots show the Dice scores (A) and HD95 distance (B)
for four brain regions on the y axis, respectively. The selected four brain
regions are listed on the x axis. HD95 was cut off at 8mm for visual clarity,
which removed some outliers for the brainstem.

After considering the effects of sex on brain region volume, density plots of
the brainstem and the hippocampus volume were obtained, shown in figure 8. The
distributions are separate for the MRI and CT volumes, while the x axes are
scaled identically within the regions. The differences in the distributions
between sexes established in MRI are similar to those based-on CT. The
quantitative overview of the sex-based analysis is located in supplementary
table 3. For all sex based analyses, the direction of the volume differences
are identical with both CT and MRI based comparisons, with a mean 17%
difference overall. Of the selected regions, the brainstem shows the most
considerable percentage difference of 52% between the MRI and CT.

Figure 8: Four density plots comparing the volume distributions in the
brainstem and hippocampus across sexes for both CT and MRI.

The results for the age-based analysis are shown in figure 9 . Similarly to
the sex-based analysis, most of the MR-based results and relationships are
preserved in the CT-based analyses. In four selected regions, the volume of
the brain areas are plotted against patient age. A linear regression slope
with confidence intervals is used to visualize the effect of aging on brain
region volumes, separately for CT and MRI-segmented volumes. All regions,
except the brainstem, show the same direction in the volume change-age
relationship. For these regions, the mean absolute difference between the
regression slopes is 10%. The brainstem show non-significant regressions in
both CTs and MRIs with different slope directions. However, these differences
are negligible when adjusted for the brainstem’s size. The quantitative
summary of the age-based analysis is also located in supplementary table 3.

Figure 9: The figure presents four plots, each illustrating the relationship
between brain region volume and age in four different brain regions. Orange
points and lines denote CT data, whereas blue points and lines represent MRI
data. The points mark the volumes from individual segmentation, and the lines
indicate the linear regression trajectories for each imaging modality,
complete with their respective confidence intervals.

##  4 Discussion

In this work, we conducted a comprehensive analysis to evaluate the
performance of the SynthSeg model for CT brain segmentations against MRI-based
ones. The quantitative analysis suggests that SynthSeg delivers CT
segmentation performance similar to the range of inter observer variability in
manual segmentation, particularly in scenarios with minimal anomalies or data
censoring. Should segmentation failure occur for any reason, it is reflected
in the QC score, allowing for the rejection of unreliable segmentations.
Nonetheless, CT segmentations typically differ from MRI-based baseline, as
evidenced by the Dice score. Generally, CT QC scores are lower than those for
MRI, implying superior performance in MRI segmentations. The difference is
marginal in certain brain regions, such as the brainstem or the thalamus
(0.037 and 0.022, respectively). Still, the disparity is more pronounced in
others, such as general grey and white matter, with mean filtered QC score
differences of 0.098 and 0.077, respectively. However, these regions occupy
large portion of the brain and are nearly always affected by lesions. This
observation aligns with the expectation of lower soft tissue contrast in CT.

The application of QC scores significantly improved the Dice score and HD95
metrics by filtering out regions unsuitable for this analysis. Typically,
these are areas with abnormal tissue, such as those affected by tumors or
surgical cavities. Additionally, regions like the lateral ventricles and grey
matter exhibit highly variable Dice scores and regions that do not overlap
despite adequate QC scores. This may suggest that QC scores are not always
reliable indicators of segmentation quality in some brain regions, though they
are generally effective. However, the HD95 distance indicates that segmented
regions are always relatively close to each other, unless a certain region is
not segmented at all. Correlations between demographic factors such as sex or
age and brain volume remain consistent in MRI and CT segmentations, with a 13%
mean absolute difference for age-based analyses. This suggests that SynthSeg’s
robustness could improve the applicability of CT databases for research. The
ability to reliably segment CT comparably to MRIs opens the door to revisiting
archived clinical datasets.

Our study, while comprehensive, has certain limitations. It primarily utilizes
260 paired images from radiotherapy patients. Despite this, the ability to
exclude diseased and censored brain regions using QC scores suggests that our
findings could have broader applicability in settings requiring brain region
segmentation. It is anticipated that segmentation performance would be
enhanced in the absence of diseased brain areas. One limitation regarding
using MRI as a basis for comparison is that MRI contain distortions, which are
not present in CT. This could cause a discrepancy in the Dice scores and HD95
distances between MRI and CT. Keeping in mind that current investigations did
not correct for MRI distortions, our expectation is that if we correct for MRI
distortions beforehand, segmentation performance would improve and HD95 and
Dice would be lower. However, since the Dice scores and HD95 still fall within
a similar range as the interobserver variation, we believe that CT performance
will still be similar to the interobserver variability with these distortions
taken into account.

This work confirms that SynthSeg provides a valuable tool for utilizing CT
data in research applications. Additionally, SynthSeg and related tools could
be used for adaptive RT applications, possibly supporting these treatment
methods within RT, where MRI is not available for soft tissue contrast [16,
17, 18]. Due to the multi-center, international nature of the data used in the
study, the results should be generalizable to other datasets. Except for the
sex and age for the UMC Utrecht cohort, the rest of the data is publicly
available, allowing for full transparency over the majority of the analysis.
For future work, validating SynthSeg’s robustness with lower quality CT, such
as cone-beam CT, could determine if the model can effectively extract data
from a wide-range of applications and even archived datasets. Additionally,
the model’s applicability could be further generalized by acquiring paired
CT/MR data from healthy individuals. Our findings suggest that SynthSeg can be
effectively utilized to obtain segmentation data from CT, enhancing research
capabilities.

##  5 Conclusion

The robust version of SynthSeg demonstrates segmentation of brain CT
performance comparable to ranges found in interobserver variation. While CT-
based segmentation typically underperforms compared to MRI-based segmentation,
the model’s automated provision of QC scores is a valuable feature. It allows
for the exclusion of regions with substandard quality, whether due to low
image quality or abnormal tissue, such as tumors. The relationships between
sex, age, and brain region volumes are preserved when performed using CT,
suggesting that the segmentations offer comparable insights. Therefore,
SynthSeg could be effectively employed for segmenting CTs in scenarios where
absolute precision is not a critical requirement.

## References

  * Despotović et al. [2015] Ivana Despotović, Bart Goossens, and Wilfried Philips.  Mri segmentation of the human brain: Challenges, methods, and applications.  _Computational and Mathematical Methods in Medicine_ , 2015:1–23, 2015.  doi: 10.1155/2015/450341.
  * Lorenzen et al. [2021] Ebbe Laugaard Lorenzen, Jesper Folsted Kallehauge, Camilla Skinnerup Byskov, Rikke Hedegaard Dahlrot, Charlotte Aaquist Haslund, Trine Lignell Guldberg, Yasmin Lassen-Ramshad, Slávka Lukacova, Aida Muhic, Petra Witt Nyström, and et al.  A national study on the inter-observer variability in the delineation of organs at risk in the brain.  _Acta Oncologica_ , 60(11):1548–1554, Oct 2021.  doi: 10.1080/0284186x.2021.1975813.
  * Billot et al. [2023a] Benjamin Billot, Douglas N. Greve, Oula Puonti, Axel Thielscher, Koen Van Leemput, Bruce Fischl, Adrian V. Dalca, and Juan Eugenio Iglesias.  Synthseg: Segmentation of brain mri scans of any contrast and resolution without retraining.  _Medical Image Analysis_ , 86:102789, 2023a.  doi: 10.1016/j.media.2023.102789.
  * Billot et al. [2023b] Benjamin Billot, Colin Magdamo, You Cheng, Steven E. Arnold, Sudeshna Das, and Juan Eugenio Iglesias.  Robust machine learning segmentation for large-scale analysis of heterogeneous clinical brain mri datasets.  _Proceedings of the National Academy of Sciences_ , 120(9), 2023b.  doi: 10.1073/pnas.2216399120.
  * Armstrong et al. [2019] N.M. Armstrong, Y. An, L. Beason-Held, J. Doshi, G. Erus, L. Ferrucci, C. Davatzikos, and S.M. Resnick.  Sex differences in brain aging and predictors of neurodegeneration in cognitively healthy older adults, Jun 2019.  URL https://www.sciencedirect.com/science/article/abs/pii/S0197458019301769.
  * Ruigrok et al. [2014] Amber N.V. Ruigrok, Gholamreza Salimi-Khorshidi, Meng-Chuan Lai, Simon Baron-Cohen, Michael V. Lombardo, Roger J. Tait, and John Suckling.  A meta-analysis of sex differences in human brain structure.  _Neuroscience and Biobehavioral Reviews_ , 39:34–50, Feb 2014.  doi: 10.1016/j.neubiorev.2013.12.004.
  * Peters [2006] R Peters.  Ageing and the brain.  _Postgraduate Medical Journal_ , 82(964):84–88, Feb 2006.  doi: 10.1136/pgmj.2005.036665.
  * Thummerer et al. [2023] Adrian Thummerer, Erik van der Bijl, Arthur Galapon Jr, Joost J. C. Verhoeff, Johannes A. Langendijk, Stefan Both, Cornelis (Nico) A. T. van den Berg, and Matteo Maspero.  Synthrad2023 grand challenge dataset: Generating synthetic ct for radiotherapy.  _Medical Physics_ , 50(7):4664–4674, 2023.  doi: https://doi.org/10.1002/mp.16529.  URL https://aapm.onlinelibrary.wiley.com/doi/abs/10.1002/mp.16529.
  * Shusharina and Bortfeld [2021] N. Shusharina and T. Bortfeld.  Glioma image segmentation for radiotherapy: Rt targets, barriers to cancer spread, and organs at risk (glis-rt), 2021.  URL https://wiki.cancerimagingarchive.net/pages/viewpage.action?pageId=95224486.
  * Zolotova et al. [2023] S.V. Zolotova, A.V. Golanov, I.N. Pronin, A.V. Dalechina, A.A. Nikolaeva, A.S. Belyashova, D.Y. Usachev, E.A. Kondrateva, P.V. Druzhinina, B.N. Shirokikh, and et al.  Burdenko-gbm-progression, Apr 2023.  URL https://doi.org/10.7937/E1QP-D183.
  * Clark et al. [2013] Kenneth Clark, Bruce Vendt, Kirk Smith, John Freymann, Justin Kirby, Paul Koppel, Stephen Moore, Stanley Phillips, David Maffitt, Michael Pringle, and et al.  The cancer imaging archive (tcia): Maintaining and operating a public information repository.  _Journal of Digital Imaging_ , 26(6):1045–1057, Jul 2013.  doi: 10.1007/s10278-013-9622-7.
  * Posit team [2023] Posit team.  _RStudio: Integrated Development Environment for R_.  Posit Software, PBC, Boston, MA, 2023.  URL http://www.posit.co/.
  * Giavarina [2015] Davide Giavarina.  Understanding bland altman analysis.  _Biochemia Medica_ , 25(2):141–151, 2015.  doi: 10.11613/bm.2015.015.
  * Jia et al. [2024] Jingnan Jia, Marius Staring, and Berend C Stoel.  seg-metrics: a python package to compute segmentation metrics.  _medRxiv_ , pages 2024–02, 2024.
  * McCarthy [2023] Paul McCarthy.  Fsleyes, 9 2023.  URL https://zenodo.org/record/7038115.
  * Tseng et al. [2024] Chia-Lin Tseng, K Liang Zeng, Eric A Mellon, Scott G Soltys, Mark Ruschin, Angus Z Lau, Natalia S Lutsik, Rachel W Chan, Jay Detsky, James Stewart, and et al.  Evolving concepts in margin strategies and adaptive radiotherapy for glioblastoma: A new future is on the horizon.  _Neuro-Oncology_ , 26, Mar 2024.  doi: 10.1093/neuonc/noad258.
  * Matsuyama et al. [2022] Tomohiko Matsuyama, Yoshiyuki Fukugawa, Junichiro Kuroda, Ryo Toya, Takahiro Watakabe, Tadashi Matsumoto, and Natsuo Oya.  A prospective comparison of adaptive and fixed boost plans in radiotherapy for glioblastoma.  _Radiation Oncology_ , 17(1), Feb 2022.  doi: 10.1186/s13014-022-02007-4.
  * Guevara et al. [2023] Beatriz Guevara, Kaylie Cullison, Danilo Maziero, Gregory A. Azzam, Macarena I. De La Fuente, Karen Brown, Alessandro Valderrama, Jessica Meshman, Adrian Breto, John Chetley Ford, and et al.  Simulated adaptive radiotherapy for shrinking glioblastoma resection cavities on a hybrid mri–linear accelerator.  _Cancers_ , 15(5):1555, Mar 2023.  doi: 10.3390/cancers15051555.


# Sparse-firing regularization methods for spiking neural networks with time-
to-first spike coding

Yusuke Sakemi  Research Center for Mathematical Engineering, Chiba Institute
of Technology, Narashino, Japan  International Research Center for
Neurointelligence (WPI-IRCN), The University of Tokyo, Tokyo, Japan  Kakei
Yamamoto  Center for Advanced Intelligence Project, RIKEN, Tokyo, Japan  Takeo
Hosomi  NEC Corporation, Kawasaki, Japan  Kazuyuki Aihara  Research Center for
Mathematical Engineering, Chiba Institute of Technology, Narashino, Japan
International Research Center for Neurointelligence (WPI-IRCN), The University
of Tokyo, Tokyo, Japan

###### Abstract

The training of multilayer spiking neural networks (SNNs) using the error
backpropagation algorithm has made significant progress in recent years. Among
the various training schemes, the error backpropagation method that directly
uses the firing time of neurons has attracted considerable attention because
it can realize ideal temporal coding. This method uses time-to-first spike
(TTFS) coding, in which each neuron fires at most once, and this restriction
on the number of firings enables information to be processed at a very low
firing frequency. This low firing frequency increases the energy efficiency of
information processing in SNNs, which is important not only because of its
similarity with information processing in the brain, but also from an
engineering point of view. However, only an upper limit has been provided for
TTFS-coded SNNs, and the information-processing capability of SNNs at lower
firing frequencies has not been fully investigated. In this paper, we propose
two spike timing-based sparse-firing (SSR) regularization methods to further
reduce the firing frequency of TTFS-coded SNNs. The first is the membrane
potential-aware SSR (M-SSR) method, which has been derived as an extreme form
of the loss function of the membrane potential value. The second is the firing
condition-aware SSR (F-SSR) method, which is a regularization function
obtained from the firing conditions. Both methods are characterized by the
fact that they only require information about the firing timing and associated
weights. The effects of these regularization methods were investigated on the
MNIST, Fashion-MNIST, and CIFAR-10 datasets using multilayer perceptron
networks and convolutional neural network structures.

## Introduction

Spiking neural networks (SNNs) can process information in the form of spikes
in a manner similar to the way information is processed in the brain. SNNs are
thereby expected to be able to achieve both high computational functionality
and energy efficiency [1]. The spikes are represented as all-or-none binary
values, and how information is represented by spikes is closely related to the
information-processing mechanism in SNNs. The spike-based information
representation methods are divided into two major categories, rate coding and
temporal coding [2, 3]. In rate coding, information is contained in the
average number of spikes generated by a neuron. In this case, the firing
frequency can take approximately continuous values as a function of the input
intensities; therefore, the resulting SNNs can be treated as differentiable
models similar to an artificial neural network (ANN). Using rate coding, ANNs
can be converted to SNNs, and the high learning ability of ANNs has been
successfully transferred to SNNs [4, 5, 6]. However, when rate coding is used,
information processing in the SNNs is just an approximation of that in ANNs.
Furthermore, the precise approximation of an ANN requires many spikes, which
reduces energy efficiency when implemented in neuromorphic hardware [7]. It
has been experimentally shown that physiologically, neurons in certain brain
regions or specific neuron types exhibit extremely sparse firing
characteristics [8], and it is thought that temporal coding using not only the
firing frequency but also the firing time is realized in at least some brain
regions [9, 10, 11, 12]. Therefore, to achieve brain-like high-capacity,
energy-efficient information processing capabilities in SNNs, it is important
to use temporal coding that also considers spike timing information.

Because in temporal coding, the precise timing of spikes carries information,
it is necessary to train the SNNs directly instead of using a converted ANN.
In recent years, by incorporating deep learning techniques, it has become
possible to directly train SNNs using the backpropagation algorithm [13, 14,
15, 16]. Among the various methods proposed, methods that focus on the
displacement of the membrane potential and those that focus on the
displacement of the spike time have attracted particular attention because of
their high learning performance. In membrane potential displacement methods,
the derivative of the output spike with respect to the membrane potential is
almost always zero because the spike is a binary value. However, it is
possible to approximate this derivative using a surrogate function [17]. This
method has been proposed in various forms by various groups [18, 19, 20, 21].
It has been proven to be very flexible, works with various surrogate functions
[22], and can be used to efficiently train recurrent neural network structures
[23, 24]. Recently, it has become possible to train relatively large models
[25]. However, with few exceptions [26, 27], the neurons exhibit high firing
frequencies, and it is debatable whether the timing information is efficiently
utilized.

A timing-based learning method is a method that focuses on the displacement of
the spike time [28]. The coding most commonly used in this learning method is
time-to-first-spike (TTFS) coding, which has the property that each neuron
fires at most once [29, 30]. Because the information is contained in the
timing of a single spike and the gradient is computed directly using the spike
timing, this coding is expected to realize an ideal temporal coding. The high
learning performance of this method has been shown in various neuron models
[31, 32, 33, 34, 35, 36, 37]. Hardware implementation efforts are also
underway to achieve high power efficiency by taking advantage of its sparsity
characteristics [37, 38]. However, the constraint of firing at most once per
neuron in TTFS coding may not be sufficiently sparse in some situations. For
example, in the brain, there are many neurons that hardly fire at all [8].
Furthermore, in an extremely power-limited environment such as edge AI [39], a
sparse firing pattern that goes beyond the constraint of one firing per neuron
is desirable.

In this paper, we propose two methods to further improve the sparse firing
property of TTFS-coded SNNs. One method is derived from the loss in the value
of the membrane potential, and the other is derived from the firing
conditions. Both methods are characterized by the fact that they only require
information about the firing timing and the weights associated with it, as is
the case in timing-based learning. In the following, we describe the two
methods and show experimentally how they suppress firing effectively on the
MNIST, Fashion-MNIST, and CIFAR-10 datasets.

## Results

### Spike timing-based sparse-firing regularization methods

Figure 1: Derivation of the two SSR methods. (a) M-SSR is derived from the
membrane potential loss, where the loss occurs when the membrane potential is
larger than $\hat{v}$. Assuming that a neuron fires at $t_{i}^{(l)}$ and
$\hat{v}$ is sufficiently large, the loss occurs only in the time interval
$[\hat{t}_{i}^{(l)},t_{i}^{(l)}]$, where $\hat{t}_{i}^{(l)}$ is the time at
which the membrane potential equals $\hat{v}$. M-SSR is obtained by setting
$\hat{v}\rightarrow V_{\text{th}}$. (b) F-SSR is derived from the firing
condition. Assuming that no input spike is accepted after the neuron’s firing
time $t_{i}^{(l)}$, the membrane potential asymptotes to
$\tau_{v}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}$ at $t\rightarrow\infty$. The
F-SSR is derived by formulating the loss such that this asymptotic membrane
potential becomes small.

We first summarize the proposed spike timing-based sparse-firing
regularization (SSR) methods. SSR methods are characterized by the fact that
they only require information about the firing timing and the weights
associated with it, as is the case in ordinary timing-based learning. In this
study, we propose two SSR method variants: membrane-potential-aware SSR
(M-SSR) and firing-condition-aware SSR (F-SSR). In both cases, we add a new
regularization term to the cost function used in supervised learning to
suppress the firing. Moreover, firing events can be suppressed by adding a new
regularization term to the cost function used in supervised learning. M-SSR is
based on the idea of reducing the value of the membrane potential, which is
realized by adding the membrane potential loss $V$ as a regularization term to
the cost function. F-SSR is based on the idea of breaking the firing
conditions, which is realized by adding the firing condition loss $Q$ as a
regularization term to the cost function. Figure 1 shows the outline of each
method. For simplicity, this paper adopts commonly used leaky integrate-and-
fire (LIF) neuron models. The LIF neuron model has as parameters the time
constant of the membrane potential $\tau_{v}$ and the time constant of the
synaptic current $\tau_{I}$ (see Method). Extensions to other neuron models
are straightforward.

First, we explain M-SSR, which is based on the idea of reducing the membrane
potential value. The membrane potential loss $V$ is defined as

$\displaystyle V$ $\displaystyle=\sum_{l}\xi^{l}\sum_{i}V_{i}^{(l)},$ (1)
$\displaystyle V_{i}^{(l)}$
$\displaystyle=\frac{1}{V_{\text{th}}-\hat{v}}\int_{0}^{T}dt\left(v_{i}^{(l)}(t)-\hat{v}\right)\theta\left(v_{i}^{(l)}(t)-\hat{v}\right)\theta\left(t_{i}^{(l)}-t\right),$
(2)

where $V_{i}^{(l)}$ is the loss relating to the membrane potential trajectory
of the $l$th-layer neuron $i$ and $\xi(>0)$ is the hyperparameter for leveling
the sparsity in each layer. In addition, $V_{\text{th}}$ is the firing
threshold, and $T$ is a parameter specifying the time interval during which
firing is suppressed. Figure 1 (a) shows the loss associated with the membrane
potential trajectory of a neuron. Note that $\hat{v}$ is sufficiently large
and there is only one point $\hat{t}_{i}^{(l)}$ at which the membrane
potential equals $\hat{v}$. In this case, the loss is nonzero only during
$[\hat{t}_{i}^{(l)},t_{i}^{(l)}]$. To perform integration in Eq. (2), we need
information about the value of the membrane potential at each time step.
However, by setting $\hat{v}\rightarrow V_{\text{th}}$, we can obtain
$\hat{t}_{i}^{(l)}\rightarrow t_{i}^{(l)}$. This avoids the integral
calculation, and Eq. (2) can be solved analytically. When computing the
gradient of this integral, it is important to fix the integration range
$[\hat{t},t_{i}^{(l)}]$. If we do not treat it as a fixed value, the more
rapidly the membrane potential rises, the smaller the membrane potential loss
in Eq. (2) becomes, and thus the firing is not effectively suppressed.
Finally, we obtain the following M-SSR:

$\displaystyle
V_{i}^{(l)}=\begin{cases}\frac{1}{\color[rgb]{0,0,1}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}}\left[{\color[rgb]{0,0,1}t_{i}^{(l)}}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}t_{j}^{(l-1)}\right],\text{
for }\tau_{v}=\tau_{I}=\infty,\\\
\tau\frac{1}{{\color[rgb]{0,0,1}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}}}-V_{\text{th}}\tau^{-1}}\left[\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{\tau}\right)a_{i}^{(l)}\right],\text{
for }\tau_{v}=\infty,\tau_{I}=\tau,\\\
2\tau{\color[rgb]{0,0,1}{\alpha_{i}^{(l)}}}\left[\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}}}}{2\tau}\right)b_{i}^{(l)}-\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}}}}{\tau}\right)a_{i}^{(l)}\right],\text{
for }\tau_{v}=2\tau,\tau_{I}=\tau.\end{cases}$ (3)

Note that in the above equations, the gradients are not calculated for the
variables shown in blue (they are considered to be constants in the gradient
calculations). This corresponds to fixing the integration range
$[\hat{t},t_{i}^{(l)}]$. Constant terms not involved in the learning are
excluded. $\Gamma_{i}^{(l)}$ denotes the index set of spikes that have been
input to the $l$th-layer neuron $i$ up to firing time $t_{i}^{(l)}$. In Eq.
(3), the following variables are defined

$\displaystyle a_{i}^{(l)}$
$\displaystyle=\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right),~{}b_{i}^{(l)}=\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{2\tau}\right),$
(4) $\displaystyle\alpha_{i}^{(l)}$
$\displaystyle=\frac{2a_{i}^{(l)}}{\left(b_{i}^{(l)}+\sqrt{\left(b_{i}^{(l)}\right)^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}\right)\left(\sqrt{\left(b_{i}^{(l)}\right)^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}\right)}.$
(5)

Appendix A provides a detailed derivation. In addition, Appendix B discusses
the consistency of the M-SSR gradient (Eq. (3)) with that of the integral-form
loss (Eq. (2)) when $\hat{v}\rightarrow V_{\text{th}}$.

Next, we explain F-SSR, a method that suppresses firing based on the firing
conditions. From the non-leaky integrate-and-fire neuron model
$(\tau_{v}=\infty,~{}\tau_{I}=\tau)$, we obtain the following firing
conditions:

$\displaystyle\text{firing
condition}_{i}^{(l)}:=\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\geq
V_{\text{th}}\tau^{-1}.$ (6)

Because the firing will be suppressed if this firing condition is not
satisfied, we define the F-SSR term $Q$ as follows:

$\displaystyle Q$ $\displaystyle=\sum_{l}\xi^{l}\sum_{i}Q_{i}^{(l)}$ (7)
$\displaystyle Q_{i}^{(l)}$
$\displaystyle=\begin{cases}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)},~{}\text{if
}t_{i}^{(l)}<T\\\ 0,~{}\text{otherwise.}\\\ \end{cases}$ (8)

We note that $V_{i}^{(l)}=Q_{i}^{(l)}=0$ if the neuron does not fire.

### Numerical simulations

We trained several SNNs on the MNIST dataset [40], Fashion-MNIST dataset [41],
and CIFAR-10 dataset [42] to investigate the effect of SSR on suppressing
firing. In the experiment, in addition to the multilayer perceptron (MLP)
structure, a convolutional neural network (CNN) structure was used. The image
data in the dataset were converted to input spikes, where the intensity of
each pixel is converted to the input time of each spike (see Method). We
define sparsity as the average number of spikes per neuron per input data in a
time window $[0,t^{\text{ref}}]$, where $t^{\text{ref}}$ is the reference time
of the output layer firing time (see Method). In addition, we set
$T=t^{\text{ref}}$ in Eqs. (2) and (8), and set the firing threshold
$V_{\text{th}}$ to 1. When the integration form Eq. (2) was used, the
integration was approximated by dividing the integral by the time width
$\Delta t$.

Figure 2: Typical results with M-SSR. We trained SNNs with one hidden layer
(784-400-10) on the MNIST dataset with various values of M-SSR strengths
$\gamma_{2}$. The upper figures present raster plots for a given input data,
showing the results of the input layer (top), hidden layer (middle), and
output layer (bottom). The lower figures show the time evolution of the
membrane potentials when given the same input data, with the results for the
hidden layer (top) and the output layer (bottom). In the panels displaying the
membrane potentials in the hidden layer, only 50 neurons are shown. In all
cases, the following hyperparameters were used:
$\tau_{v}=\tau_{I}=\infty,~{}t^{\text{ref}}=8,~{}\gamma_{1}=10^{-4},\gamma_{3}=0,\eta=10^{-4}$,
and $\tau_{\text{soft}}=0.9$.

Figure 2 shows the learning results of SNNs with one hidden layer (784-400-10)
trained on the MNIST dataset with various M-SSR strengths $\gamma_{2}$ (see
Method). Note that all output layer neurons were required to fire because the
loss function was defined by the spike timing of the neurons in the output
layer. Therefore, sparse firing regularization was applied only to the hidden
layers. The upper figures show the raster plots of the firing distribution of
each layer for a given input data, and the lower figures show the time
evolution of the membrane potentials of each layer. As the strength of M-SSR
was increased, the number of neurons that fired tended to decrease, and it can
be seen that most of the hidden layer neurons stopped firing when
$\gamma_{2}=1.3\times 10^{-5}$. By contrast, the firing distribution of the
output layer did not change significantly with respect to M-SSR strength. The
neurons corresponding to the correct index fired the earliest $(t\sim 4)$, and
the other neurons fired later $(t\sim 8)$. The membrane potentials of the
hidden layer neurons were suppressed as the regularization strength increased,
whereas the output layer solved the task using fewer spikes from the hidden
layer. This indicates that M-SSR regularization can suppress the firing of the
hidden layer without significantly compromising recognition task performance.

Figure 3: Sparsity–accuracy tradeoff for different regularization forms. We
evaluated the integral-form regularization (Eq. (2)) and M-SSR (Eq. (3)) in a
784-400-10 SNN based on various neuron models. The tradeoff in
sparsity–accuracy is shown when the regularization strength $\gamma_{2}$ was
varied from $0$ to $10^{-4}$. The standard deviations were obtained over 10
trials. The hyperparameters were
$t^{\text{ref}}=8,\gamma_{1}=10^{-4},\gamma_{3}=0,\eta=10^{-4}$, and
$\tau_{\text{soft}}=0.9$. For the integral-form regularization, we used
various values of $\hat{v}$, and we set $\Delta t=t^{\text{ref}}/1000$.

Figure 3 shows the sparsity–accuracy tradeoff results obtained when using the
integral-form regularization (Eq. (2)) and M-SSR (Eq. (3)). We trained SNNs
with a single hidden layer (784-400-10) using various regularization
strengths. The standard deviations were obtained over 10 trials. The upper
figures show the results for the MNIST dataset, and the lower figures show the
results for the Fashion-MNIST dataset. Results are also shown for different
neuron models $(\tau_{v},\tau_{I})$. For the MNIST dataset, the tradeoff
curves show that a larger $\hat{v}$ led to a better tradeoff for all neuron
models, and the best tradeoff was obtained by M-SSR (corresponding to
$\hat{v}=1$). Similar results were obtained for Fashion-MNIST, although the
advantage was not as pronounced as it was for MNIST. This result demonstrates
that the integral-form regularization (Eq. (2)) smoothly transitioned to the
limit form (M-SSR, Eq. (3)). Furthermore, taking the limit of
$\hat{v}\rightarrow 1$, the tradeoff between sparsity and accuracy was
improved. In addition, good sparsity–accuracy tradeoff properties were
obtained for various neuron models ($\tau_{v}$ and $\tau_{I}$), suggesting
that M-SSR is preferable to the integral-form regularization for TTFS-coded
SNNs.

Figure 4: Comparison of SSR methods on the MNIST and Fashion-MNIST datasets.
The top figures show the sparsity–accuracy tradeoffs for the 2-layer SNNs
(784-400-10), whereas the bottom figures show the sparsity–accuracy tradeoffs
for the 4-layer SNNs (784-400-400-400-10). In the bottom figures, each panel
represents the accuracy–sparsity tradeoff for the first, second, and third
hidden layers, from the top. The bottom panel presents the sparsity–accuracy
tradeoffs for the sparsity averaged over the three hidden layers. The standard
deviations were obtained over 10 trials. We used the following
hyperparameters: $t^{\text{ref}}=8,\gamma_{1}=10^{-4},\eta=10^{-4}$, and
$\tau_{\text{soft}}=0.9$ for the 2-layer SNNs,
$t^{\text{ref}}=9,\gamma_{1}=10^{-4},\eta=10^{-4}$, and
$\tau_{\text{soft}}=0.9$ for the 4-layer SNNs. For the MNIST dataset, we set
$\xi$ to 6 for M-SSR, and 1 for F-SSR. For the Fashion-MNIST dataset, we set
$\xi$ to 6 for M-SSR, and 4 for F-SSR.

Figure 4 compares the results of the two proposed SSRs, M-SSR (Eq. (3)) and
F-SSR (Eq. (8)). The results for an SNN with one hidden layer (784-400-10) are
shown in the top figures. On the MNIST benchmark, F-SSR obtained a better
sparsity–accuracy tradeoff than M-SSR, whereas on the Fashion-MNIST benchmark,
both F-SSR and M-SSR yielded a similar sparsity–accuracy tradeoff. The results
for an SNN with three hidden layers are shown in the lower figures. For each
SSR method for each dataset, the value of $\xi$ (Eq. (1)) was set to obtain
the best tradeoff averaged over the whole layer (except for the output layer).
See Appendix C for the tradeoff properties for various values of $\xi$. For an
SNN with three hidden layers, M-SSR and F-SSR showed similar sparsity–accuracy
tradeoff characteristics, but the optimal value of $\xi$ differed
significantly in M-SSR and F-SSR. For the MNIST dataset, the optimal value was
$\xi=6$ for M-SSR and $\xi=1$ for F-SSR, whereas for the Fashion-MNIST
dataset, the optimal value was $\xi=6$ for M-SSR and $\xi=4$ for F-SSR. This
difference may be due to the characteristics of the regularization function.
In M-SSR (Eq. (3)), the error that occurred in the $l$-th layer propagates
back to the previous layer via spike timing $t_{j}^{(l-1)}$. In this case, if
the weights of neurons from $t_{j}^{(l-1)}$ to the $l$th layer were positive
overall, $t_{j}^{(l-1)}$ increased during training, and consequently the
$l-1$th layer also became more sparse. Similarly, the $l-2$th layer was
expected to become sparse. Therefore, to counteract this effect, a relatively
large value of $\xi$ was optimal. By contrast, in the case of F-SSR (Eq. (8)),
the losses that occurred in the $l$th layer do not propagate back to previous
layers. Therefore, a relatively small value of $\xi$ was optimal.

Table 1: Convolutional architectures used for each dataset. “Conv($a$,$b$)” represents a convolutional layer with a kernel size of $a\times a$, number of output channels $b$, and stride 1. “Pool” represents a pooling layer with kernel size of $2\times 2$, stride 2. The padding of the convolutional layer was set to 0 for MNIST and set to 1 for Fashion-MNIST and CIFAR-10.  Dataset | Network
---|---
MNIST | Conv(5, 6)-Pool-Conv(5, 16)-Pool-400-400-10
Fashion-MNIST | Conv(5, 6)-Pool-Conv(5, 16)-Pool-400-400-10
CIFAR-10 | Conv(3, 24)-Pool-Conv(3, 48)-Pool-Conv(3, 96)-Pool-600-10
Figure 5: Effects of SSR with a CNN-architecture on the MNIST dataset (a) and
the Fashion-MNIST dataset (b). From top to bottom, the panels represent the
sparsity–accuracy tradeoffs for the first convolutional layer, second
convolutional layer, first fully connected layer, second fully connected
layer, and whole network (not including pooling layers). The standard
deviations were obtained over 10 trials. We used the following
hyperparameters: $t^{\text{ref}}=16,\gamma_{1}=10^{-4},\eta=10^{-4},$ and
$\tau_{\text{soft}}=0.9$. For the MNIST dataset, we set $\xi$ to 6 for M-SSR,
and 2 for F-SSR. For the Fashion-MNIST dataset, we set $\xi$ to 6 for M-SSR,
and 4 for F-SSR.

Next, we applied the SSR methods to spiking CNNs. Table 1 shows the network
structure used for each dataset. Figure 5 shows the effect of SSR
regularization on MNIST and Fashion-MNIST. The overall sparsity–accuracy
tradeoff for the CNN structure was worse than that for the MLP structure. The
SNNs with MLP structures reduced the average number of firings per neuron to
about $0.1-0.2$ with almost no loss in accuracy, whereas the SNNs with
convolutional structures only reduced the number of firings to about $0.4$. It
can be seen that the first convolutional layers are not very sparse, with the
exception of the results of F-SSR on the MNIST dataset. This is considered to
be caused by the fact that it is difficult to force only a portion of the
neurons to fire in a convolutional layer because of the weight-sharing
property.

Figure 6: Sparsity–accuracy tradeoff on the CIFAR-10 task. The best case over
10 trials is plotted. From the top, the sparsity–accuracy tradeoffs for the
first, second, and third convolutional layers, and then fully connected layer
are represented, respectively. The sparsity—accuracy tradeoff for the whole
network (not including pooling layers) is presented in the bottom panel. We
inverted the effects of F-SSR, so the neurons tend to fire more frequently.

On the CIFAR-10 dataset, firing tended to be suppressed during training even
when SSR was not applied, and we confirmed that neurons relating to some
channels did not fire at all image locations and for all training data. To
avoid this problem, we trained the SNNs to attempt to satisfy the firing
conditions. This was achieved by making the regularization strength
$\gamma_{3}$ of the F-SSR term negative (see Method). The results of learning
by promoting firing are shown in Fig.6. Because the variance over the 10
trials was large, we plotted the results for the trial with the highest test
accuracy for better visualization. Appendix D presents the results for all
trials. By sacrificing sparsity, we were able to observe a noticeable
improvement in performance. We conducted an experiment incorporating M-SSR
regularization, but the combination of promoting firing and sparse firing did
not yield a favorable sparsity-accuracy tradeoff.

## Discussion

SNNs with TTFS coding can realize ideal temporal coding by constraining each
neuron to fire at most once. Due to this mechanism, the SNNs with TTFS coding
have high firing sparsity, and this approach has been applied in energy-
efficient hardware implementations [37, 38]. To further improve this sparse
firing characteristic, we developed the SSR methods. The two SSR methods were
derived from two different perspectives. The first one is M-SSR, which was
derived by assigning a penalty each time the membrane potential exceeds a
threshold $\hat{v}$ and taking the limit when the threshold equals the firing
threshold. The other is F-SSR, which was obtained from the firing conditions
of neurons. Both SSR methods are characterized by the fact that they do not
require information about the membrane potential itself, only the firing time
and associated weights. The sparsity–accuracy properties of these two methods
were investigated using the MNIST, Fashion-MNIST, and CIFAR-10 datasets.
Interestingly, although some differences were observed depending on the
datasets and the network structure, both F-SSR and M-SSR showed equally good
sparsity–accuracy properties, even though the regularization methods were
derived from different perspectives. In particular, for the fully connected
layer, it was found that the average number of firings for each neuron could
be lowered to 0.1 to 0.2. From the experiments conducted in this study, it is
difficult to determine which method is superior. We can at least conclude that
F-SSR has the advantage of a somewhat smaller computational load than M-SSR
due to its simpler formula. To understand the difference between F-SSR and
M-SSR, in addition to the sparsity–accuracy property, a detailed analysis of
the changes in the firing characteristics and in the information processing
mechanism associated with the sparse firing mechanism will be required in
future.

For the CNN structures, we found the SSR methods had more difficulty
suppressing the firing of neurons in the convolutional layer than in the fully
connected layer. In CIFAR-10 in particular, we observed that firing is
suppressed too much and learning becomes difficult even without SSR. This may
be because it is difficult to flexibly decide whether the outputs belonging to
a certain kernel should fire depending on the position because of the weight-
sharing property. To prevent this, we found that, in CIFAR-10, learning
performance can be improved by promoting firing. Similar firing promotion
terms have been introduced in previous studies [31, 43].In timing-based
learning of large-scale CNN structures, one way to obtain better
sparsity–accuracy properties is to utilize models that allow multiple neuron
firings [44] combined with the SSR methods.

Previous studies have developed methods that suppress the firing of SNNs in
the framework of the surrogate gradient method [26, 27]. They applied direct
regularization to the spike variable $s(t)\in\\{0,1\\}$ represented at each
time step in the model to the time-discretized SNN. The gradient calculation
is made possible using the surrogate function
$\frac{ds(t)}{dv(t)}=\sigma(v(t))$ [17]. This method is closely related to the
M-SSR proposed in this paper. In the surrogate gradient method, the spike
variables above are treated as a function of the membrane potential
$s(t)=\int_{-\infty}^{v(t)}\sigma(v^{\prime})dv^{\prime}$. In this sense, the
idea is similar to the loss in Eq. (2), which integrates the membrane
potential. By contrast, M-SSR, unlike the previous method [26, 27], can be
transformed from the time-integration form to the timing form by setting
$\hat{v}\rightarrow V_{\text{th}}$. This may correspond to the fact that the
learning method with the surrogate function can transition to a timing-like
learning method by taking a limit [19]. Interestingly, as shown in Fig.3, the
sparsity–accuracy property improves as $\hat{v}$ gets closer to
$V_{\text{th}}$. This suggests that timing-based sparse regularization is more
effective in timing-based learning. We note that F-SSR is a regularization
method using firing conditions, which is unique to timing-based learning.

SNNs can operate efficiently on neuromorphic hardware [26, 45]. Because the
energy consumed by the spike transmission increases as the firing frequency
increases, reducing the firing frequency is an important issue in real-world
applications. SNNs with TTFS coding are expected to provide significant power
advantages in hardware implementation due to their extremely sparse firing
[34, 35, 46]. Several research groups have reported hardware implementations
of such SNNs [38, 37, 47]. The SSR methods are expected to further improve the
energy efficiency of SNNs. Moreover, unlike the methods in [26, 27], the SSR
methods can calculate the gradient without observing the membrane potential,
which may simplify the learning system on hardware. Finally, in addition to
the reduction in the firing rate, the combination of binarized weights [48]
and pruned weights [49, 27, 50] is expected to make the SNN model more
suitable for hardware implementation.

## Method

### SNN models

In this study, we constructed a multilayer SNN using the following LIF neuron
model:

$\displaystyle\frac{d}{dt}v_{i}^{(l)}(t)$
$\displaystyle=\frac{1}{\tau_{v}}v_{i}^{(l)}(t)+I_{i}^{(l)},$ (9)
$\displaystyle\frac{d}{dt}I_{i}^{(l)}(t)$
$\displaystyle=\frac{1}{\tau_{I}}I_{i}^{(l)}(t)+\sum_{j=1}^{N^{(l-1)}}w_{ij}^{(l)}\delta(t-t_{j}^{(l-1)}),$
(10)

where $v_{i}^{(l)}$ is the membrane potential of neuron $i$ in the $l$th
layer, $w_{ij}^{(l)}$ is the coupling strength from neuron $j$ in the $l-1$th
layer to neuron $i$ in the $l$th layer, and $t_{j}^{(l-1)}$ is the firing time
of neuron $j$ in the $l-1$th layer. Furthermore, $\tau_{v}$ is the time
constant of the membrane potential and $\tau_{I}$ is the time constant of the
synaptic current. $N^{(l)}$ is the number of neurons that make up the $l$th
layer. Neurons fire when the membrane potential reaches the firing threshold
$V_{\text{th}}$ and generate spikes. After firing, the membrane potential is
fixed at $0$ and never fires again. The membrane potential of the model
described in Eq. (9) is analytically obtained as follows:

$\displaystyle v_{i}^{(l)}(t)$
$\displaystyle=\frac{\tau_{v}\tau_{I}}{\tau_{v}-\tau_{I}}\sum_{j=1}^{N^{(l-1)}}w_{ij}^{(l)}\kappa(t-t_{j}^{(l-1)}),$
(11) $\displaystyle\kappa(t)$
$\displaystyle=\theta(t)\left[\exp\left(-\frac{t}{\tau_{v}}\right)-\exp\left(-\frac{t}{\tau_{I}}\right)\right],$
(12) $\displaystyle\theta(t)$ $\displaystyle=\begin{cases}0,\text{ for
}t<0,\\\ 1,\text{ for }0\leq t.\end{cases}$ (13)

The experiments in this study consider the three cases
$(\tau_{v},\tau_{I})\in\\{(2\tau,\tau),(\infty,\tau)$, and
$(\infty,\infty)\\}$. Note that the learning characteristics of SNNs with TTFS
coding were investigated for the cases of $\tau_{v}=\infty,\tau_{I}=\infty$,
[34], $\tau_{v}=\infty$ [31], and $\tau_{v}\neq\infty,\tau_{I}\neq\infty$ [33,
37].

The firing time in each case can be calculated from the condition
$v_{i}^{(l)}(t_{i}^{(l)})=V_{\text{th}}$ as follows:

$\displaystyle
t_{i}^{(l)}=\begin{cases}\frac{V_{\text{th}}+\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}t_{j}^{(l)}}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}},\text{
for }(\tau_{v},\tau_{I})=(\infty,\infty),\\\
\tau\ln\left[\frac{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-V_{\text{th}}\tau^{-1}}\right],\text{
for }(\tau_{v},\tau_{I})=(\infty,\tau),\\\
2\tau\ln\left[\tau\frac{b_{i}^{(l)}-\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}}{V_{\text{th}}}\right],\text{
for }(\tau_{v},\tau_{I})=(2\tau,\tau),\end{cases}$ (14)

where $\Gamma_{i}^{(l)}$ denotes the index set of spikes input to the $l$th
layer neuron $i$ up to firing time $t_{i}^{(l)}$. We define the following
variables:

$\displaystyle
a_{i}^{(l)}=\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right),~{}b_{i}^{(l)}=\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{2\tau}\right).$
(15)

A detailed derivation of the firing time in the case of
$(\tau_{v},\tau_{I})=(\tau,2\tau)$ is given in Appendix A.3.

### Learning algorithms

Supervised learning of the SNN was performed using the following cost function

$\displaystyle C$
$\displaystyle=L(t^{(M)})+\gamma_{1}T(t^{(M)})+\gamma_{2}V+\gamma_{3}Q,$ (16)
$\displaystyle L$ $\displaystyle=\sum_{i=1}^{N^{(M)}}\kappa_{i}\ln S_{i},$
(17) $\displaystyle S_{i}$
$\displaystyle=\frac{\exp\left(\frac{t_{i}^{(M)}}{\tau_{\text{soft}}}\right)}{\sum_{j=1}^{N^{(M)}}\exp\left(\frac{t_{j}^{(M)}}{\tau_{\text{soft}}}\right)},$
(18) $\displaystyle T$
$\displaystyle=\sum_{i=1}^{N^{(M)}}\left(t_{i}^{(M)}-t^{\text{ref}}\right)^{2},$
(19)

where $M$ represents the output layer and
$t^{(M)}=\left(t_{1}^{(M)},t_{2}^{(M)},\dots,t_{N^{(M)}}^{(M)}\right)$. The
value of the teacher label $\kappa_{i}$ is equal to one when the $i$th label
is assigned and zero otherwise. Parameters $\gamma_{1}$, $\gamma_{2}$, and
$\gamma_{3}$ are real numbers, and they respectively control the signiﬁcance
of the temporal penalty term $T$, the membrane potential loss $V$ (Eq. (1)),
and the firing condition term $Q$ (Eq. (7)). Parameter $\tau_{\text{soft}}$ is
a positive real number, which adjusts the softmax scaling. Learning was
performed by minimizing this cost function using the gradient method with the
Adam optimizer [51] at a learning rate of $\eta$. On the CIFAR-10 task, the
coefficient $\gamma_{3}$ was set to a negative number to promote firing. In
this case, the firing condition term $Q$ was modified as follows:

$\displaystyle Q_{i}^{(l)}$
$\displaystyle=\begin{cases}\sum_{i=1}^{N^{(l-1)}}w_{ij}^{(l)},~{}\text{if not
fired},\\\ 0,~{}\text{otherwise.}\\\ \end{cases}$ (20)

We promoted the neurons to fire only when the corresponding neurons did not
fire.

### Dataset

The MNIST, Fashon-MNIST, and CIFAR-10 datasets include 2-dimensional image
data. In the MNIST and Fashion-MNIST datasets, each image has one channel,
whereas in the CIFAR-10 dataset, the images have three channels. To process
such image data, we first normalized the pixel intensity to $[0,1]$. Then, we
obtained an input spike as follows:

$\displaystyle t_{ijk}^{(0)}=\tau_{\text{in}}(1-x_{ijk}),$ (21)

where $x_{ijk}$ is the normalized pixel intensity, the first and second
indices represent the coordinates of the pixel, and the third index represents
the channel number. Here, $\tau_{\text{in}}$ is a positive constant. We set
$\tau_{\text{in}}=5$ in all experiments. When spikes are input to a fully
connected layer, the input tensors are reshaped into one-dimensional tensors.
For the CIFAR-10 dataset, to avoid the problem of the first hidden layer
firing too early and ignoring later inputs, the number of channels was doubled
as follows:

$\displaystyle x_{i,j,k}=1-x_{i,j,k-3}~{}(k=3,4,5).$ (22)

Furthermore, we used data augmentation (horizontal flipping, rotation, and
cropping) as in the previous study [43].

## Acknowledgements

This work was partially supported by JST PRESTO Grant Number JPMJPR22C5, NEC
Corporation, SECOM Science and Technology Foundation, JST Moonshot R&D Grant
Number JPMJMS2021, AMED under Grant Number JP23dm0307009, the International
Research Center for Neurointelligence (WPI-IRCN) at The University of Tokyo
Institutes for Advanced Study (UTIAS), JSPS KAKENHI Grant Number JP20H05921.
Computational resource of AI Bridging Cloud Infrastructure (ABCI) provided by
National Institute of Advanced Industrial Science and Technology (AIST) was
used.

## Appendices

##  Appendix A Derivation of M-SSR (Eq. (3))

In this section, we present a detailed derivation of M-SSR (Eq. (3)). If
$\hat{v}$ is sufficiently close to $V_{\text{th}}$, $\hat{t}_{i}^{(l)}$, the
time at which the membrane potential is $\hat{v}$ can be assumed to be a
single point (see Fig. 1 (a)). Moreover, we can assume that the number of
input spikes can be assumed to be constant $\left(\Gamma_{i}^{(l)}\right)$ in
the time interval $[\hat{t}_{i}^{(l)},t_{i}^{(l)}]$. From the above, the
integral (Eq. (2)) can be transformed as follows

$\displaystyle V_{i}^{(l)}$
$\displaystyle=\frac{1}{V_{\text{th}}-\hat{v}}\int_{0}^{T}dt\left(v_{i}^{(l)}(t)-\hat{v}\right)\theta\left(v_{i}^{(l)}(t)-\hat{v}\right)\theta\left(t_{i}-t\right),$
(23)
$\displaystyle=\frac{1}{V_{\text{th}}-\hat{v}}\int_{\hat{t}_{i}^{(l)}}^{t_{i}^{(l)}}dt\left(v_{i}^{(l)}(t)-\hat{v}\right).$
(24)

Note that $V_{i}^{(l)}=0$ if the neuron does not fire. Importantly, when
computing the gradient of this integral, the integration range
$[\hat{t},t_{i}^{(l)}]$ should be fixed. This is because if they are not
treated as fixed values, the SNNs learn to make the membrane potential rise
rapidly, thus the effect of firing suppression cannot be obtained. We can omit
the constant term $\hat{v}(t_{i}^{(l)}-\hat{t}_{i}^{(l)})$, which is not
involved in the learning process, and we only need to calculate the following:

$\displaystyle
V_{i}=\frac{1}{V_{\text{th}}-\hat{v}}\int_{\hat{t}_{i}^{(l)}}^{t_{i}^{(l)}}v_{i}^{(l)}(t)dt.$
(25)

In the following, the limit of the above integral $(\hat{v}\rightarrow
V_{\text{th}})$ is calculated for various neuron models
($\left(\tau_{v},\tau_{I}\right)=\left(\infty,\infty\right),\left(\infty,\tau\right),\text{and}\left(2\tau,\tau\right)$).

###  A.1 Neuron model with $\tau_{v}=\infty$ and $\tau_{I}=\infty$

When $\tau_{I}=\tau_{v}=\infty$, the membrane potential and firing time can be
calculated as follows:

$\displaystyle v_{i}^{(l)}(t)$
$\displaystyle=\sum_{j=1}^{N^{(l-1)}}w_{ij}^{(l)}(t-t_{j}^{(l-1)})\theta(t-t_{j}^{(l-1)}),$
(26) $\displaystyle t_{i}^{(l)}$
$\displaystyle=\frac{V_{\text{th}}+\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}t_{j}^{(l)}}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}},$
(27) $\displaystyle\hat{t}_{i}^{(l)}$
$\displaystyle=\frac{\hat{v}+\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}t_{j}^{(l)}}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}},$
(28)
$\displaystyle\frac{t_{i}^{(l)}-\hat{t}_{i}^{(l)}}{V_{\text{th}}-\hat{v}}$
$\displaystyle=\frac{1}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}}.$ (29)

Using these equations, the limit $(\hat{v}\rightarrow V_{\text{th}})$ can be
calculated as follows:

$\displaystyle V_{i}^{(l)}$
$\displaystyle=\frac{1}{V_{\text{th}}-\hat{v}}\int_{\hat{t}_{i}^{(l)}}^{t_{i}^{(l)}}dtv_{i}^{(l)}(t)$
(30)
$\displaystyle=\frac{1}{V_{\text{th}}-\hat{v}}\int_{\hat{t}_{i}^{(l)}}^{t_{i}^{(l)}}dt\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\left(t-t_{j}^{(l-1)}\right)$
(31)
$\displaystyle=\frac{1}{V_{\text{th}}-\hat{v}}\left[\frac{1}{2}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\left(t-t_{j}^{(l-1)}\right)^{2}\right]_{\hat{t}_{i}^{(l)}}^{t_{i}^{(l)}}$
(32)
$\displaystyle=\frac{1}{2\left(V_{\text{th}}-\hat{v}\right)}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\left[\left({\color[rgb]{0,0,1}t_{i}^{(l)}}-t_{j}^{(l-1)}\right)^{2}-\left({\color[rgb]{0,0,1}\hat{t}_{i}^{(l)}}-t_{j}^{(l-1)}\right)^{2}\right]$
(33)
$\displaystyle=\frac{1}{2}\frac{{\color[rgb]{0,0,1}t_{i}^{(l)}}-{\color[rgb]{0,0,1}\hat{t}_{i}^{(l)}}}{V_{\text{th}}-\hat{v}}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}\left({\color[rgb]{0,0,1}t_{i}^{(l)}}+{\color[rgb]{0,0,1}\hat{t}_{i}^{(l)}}-2t_{j}^{(l-1)}\right)$
(34)
$\displaystyle\overset{\text{Eq.}(\ref{eq:t_diff_M_SSR})}{=}\frac{1}{\color[rgb]{0,0,1}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}}\left[{\color[rgb]{0,0,1}t_{i}^{(l)}}\sum_{j}w_{ij}-\sum_{j}w_{ij}t_{j}^{(l-1)}\right].$
(35)

Note that the blue variables are related to the integral range, which is
treated as a fixed value when calculating the gradient, as mentioned above.

###  A.2 Neuron model with $\tau_{v}=\infty$ and $\tau_{I}=\tau$

When $\tau_{v}=\infty$ and $\tau_{I}=\tau$, the membrane potential and firing
time are given by [31]

$\displaystyle v_{i}^{(l)}(t)$
$\displaystyle=\tau\sum_{j=1}^{N^{(l-1)}}w_{ij}^{(l)}\theta(t-t_{j}^{(l-1)})\left[1-\exp\left(-\frac{t-t_{j}^{(l-1)}}{\tau}\right)\right],$
(36) $\displaystyle t_{i}^{(l)}$
$\displaystyle=\tau\ln\left[\frac{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-V_{\text{th}}\tau^{-1}}\right],$
(37) $\displaystyle\hat{t}_{i}^{(l)}$
$\displaystyle=\tau\ln\left[\frac{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-\hat{v}\tau^{-1}}\right].$
(38)

We also obtain the following in $\hat{v}\rightarrow V_{\text{th}}$:

$\displaystyle\frac{t_{i}^{(l)}-\hat{t}_{i}^{(l)}}{V_{\text{th}}-\hat{v}}$
$\displaystyle=\frac{\tau}{V_{\text{th}}-\hat{v}}\ln\frac{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l)}}{\tau}\right)}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-V_{\text{th}}\tau^{-1}}-\frac{\tau}{V_{\text{th}}-\hat{v}}\ln\frac{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l)}}{\tau}\right)}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-\hat{v}\tau^{-1}}$
(39)
$\displaystyle=\frac{\tau}{V_{\text{th}}-\hat{v}}\ln\frac{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-\hat{v}\tau^{-1}}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-V_{\text{th}}\tau^{-1}}$
(40)
$\displaystyle=\frac{\tau}{V_{\text{th}}-\hat{v}}\ln\left(1+\frac{\left(V_{\text{th}}-\hat{v}\right)\tau^{-1}}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-V_{\text{th}}\tau^{-1}}\right)$
(41) $\displaystyle\overset{\hat{v}\rightarrow
V_{\text{th}}}{=}\frac{1}{\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}-V_{\text{th}}\tau^{-1}}.$
(42)

Using these, the limit $(\hat{v}\rightarrow V_{\text{th}})$ can be calculated
as follows:

$\displaystyle\frac{\int
dtv_{i}^{(l)}(t)}{V_{\text{th}}-\hat{v}}=\frac{\tau}{V_{\text{th}}-\hat{v}}\left[\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\left\\{t+\tau\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)\exp\left(-\frac{t}{\tau}\right)\right\\}\right]^{t_{i}^{(l)}}_{\hat{t}_{i}^{(l)}}$
(43)
$\displaystyle=\frac{\tau}{V_{\text{th}}-\hat{v}}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\left[({\color[rgb]{0,0,1}{t_{i}^{(l)}}}-{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}})+\tau\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)\left\\{\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{\tau}\right)-\exp\left(-\frac{{\color[rgb]{0,0,1}{\hat{t}_{i}}^{(l)}}}{\tau}\right)\right\\}\right]$
(44)
$\displaystyle=\frac{\tau({\color[rgb]{0,0,1}{t_{i}^{(l)}}}-{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}})}{V_{\text{th}}-\hat{v}}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\left[1+\frac{\tau}{({\color[rgb]{0,0,1}{t_{i}^{(l)}}}-{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}})}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{\tau}\right)\left\\{1-\exp\left(\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}-{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}}}{\tau}\right)\right\\}\right]$
(45) $\displaystyle\overset{\hat{v}\rightarrow
V_{\text{th}}}{=}\tau{\frac{1}{{\color[rgb]{0,0,1}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}}-V_{\text{th}}\tau^{-1}}}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\left\\{1-\exp\left(\frac{t_{j}^{(l)}}{\tau}\right)\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{\tau}\right)\right\\}$
(46)
$\displaystyle=\tau{\frac{1}{{\color[rgb]{0,0,1}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}}-V_{\text{th}}\tau^{-1}}}\left[\left(\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\right)-\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{\tau}\right)\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)\right].$
(47)

###  A.3 Alpha-synaptic neuron model with $\tau_{v}=2\tau_{I}=2\tau$

When $\tau_{v}=2\tau_{I}=2\tau$, the membrane potential is given by

$\displaystyle v_{i}^{(l)}(t)$
$\displaystyle=2\tau\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\theta(t-t_{j}^{(l-1)})\left[\exp\left(-\frac{t-t_{j}^{(l-1)}}{2\tau}\right)-\exp\left(-\frac{t-t_{j}^{(l-1)}}{\tau}\right)\right].$
(48)

From the firing condition $v_{i}^{(l)}\left(t_{i}^{(l)}\right)=V_{\text{th}}$,
we obtain the following:

$\displaystyle\exp\left(-\frac{t_{i}^{(l)}}{2\tau}\right)\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{2\tau}\right)-\exp\left(-\frac{t_{i}^{(l)}}{\tau}\right)\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)$
$\displaystyle=\frac{V_{\text{th}}}{2\tau}$ (49)
$\displaystyle\left[\exp\left(-\frac{t_{i}^{(l)}}{2\tau}\right)\right]^{2}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)-\exp\left(-\frac{t_{i}^{(l)}}{2\tau}\right)\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{2\tau}\right)+\frac{V_{\text{th}}}{2\tau}$
$\displaystyle=0$ (50) $\displaystyle
a_{i}\left[\exp\left(-\frac{t_{i}^{(l)}}{2\tau}\right)\right]^{2}-b_{i}^{(l)}\exp\left(-\frac{t_{i}^{(l)}}{2\tau}\right)+\frac{V_{\text{th}}}{2\tau}$
$\displaystyle=0,$ (51)

where we defined the following variables:

$\displaystyle
a_{i}^{(l)}=\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right),~{}b_{i}^{(l)}=\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{2\tau}\right).$
(52)

From the formula to solve a quadratic function, the firing time $t_{i}^{(l)}$
can be calculated from Eq. (51) as follows [37]:

$\displaystyle\exp\left(-\frac{t_{i}^{(l)}}{2\tau}\right)$
$\displaystyle=\frac{b_{i}^{(l)}+\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}}{2a_{i}^{(l)}}$
(53) $\displaystyle t_{i}^{(l)}$
$\displaystyle=-2\tau\ln\left[\frac{b_{i}^{(l)}+\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}}{2a_{i}^{(l)}}\right],$
(54)

where the other solution of the quadratic function is ignored because it
indicates the time at which the membrane potential decreases from a value
greater than $V_{\text{th}}$ to a smaller value [37]. We also obtain the
following:

$\displaystyle\frac{t_{i}^{(l)}-\hat{t}_{i}^{(l)}}{V_{\text{th}}-\hat{v}}$
$\displaystyle=\frac{2\tau}{V_{\text{th}}-\hat{v}}\ln\frac{b_{i}^{(l)}+\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}\hat{v}}}{b_{i}^{(l)}+\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}}$
(56)
$\displaystyle=\frac{2\tau}{V_{\text{th}}-\hat{v}}\ln\frac{b_{i}^{(l)}+\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}+2a_{i}^{(l)}\tau^{-1}\left(V_{\text{th}}-\hat{v}\right)}}{b_{i}^{(l)}+\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}}$
(57) $\displaystyle\overset{\hat{v}\rightarrow
V_{\text{th}}}{=}2\tau\frac{a_{i}^{(l)}\tau^{-1}}{\left(b_{i}^{(l)}+\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}\right)\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}}$
(58) $\displaystyle=\alpha_{i}^{(l)}.$ (59)

We define the following variable:

$\displaystyle\alpha_{i}^{(l)}=\frac{2a_{i}}{\left(b_{i}^{(l)}+\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}\right)\left(\sqrt{(b_{i}^{(l)})^{2}-2a_{i}^{(l)}\tau^{-1}V_{\text{th}}}\right)}.$
(60)

Using these, the limit $(\hat{v}\rightarrow V_{\text{th}})$ can be calculated
as follows:

$\displaystyle\frac{\int
v_{i}^{(l)}(t)dt}{V_{\text{th}}-\hat{v}}=\frac{2\tau}{V_{\text{th}}-\hat{v}}\sum_{j\in\Gamma_{i}^{(l)}}\left[w_{ij}\left\\{-2\tau\exp\left(-\frac{t-t_{j}^{(l)}}{2\tau}\right)+\tau\exp\left(-\frac{t-t_{j}^{(l-1)}}{\tau}\right)\right\\}\right]^{t_{i}^{(l)}}_{{\hat{t}_{i}^{(l)}}}$
(61)
$\displaystyle=\frac{2\tau^{2}}{V_{\text{th}}-\hat{v}}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\left[-2\exp\left(\frac{t_{j}^{(l-1)}}{2\tau}\right)\exp\left(-\frac{t}{2\tau}\right)+\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)\exp\left(-\frac{t}{\tau}\right)\right]^{t_{i}^{(l)}}_{{\hat{t}_{i}^{(l)}}}$
(62)
$\displaystyle=\frac{2\tau^{2}}{V_{\text{th}}-\hat{v}}\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\Bigg{[}-2\exp\left(\frac{t_{j}^{(l-1)}}{2\tau}\right)\left\\{\exp\left(-\frac{\color[rgb]{0,0,1}{t_{i}^{(l)}}}{2\tau}\right)-\exp\left(-\frac{{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}}}{2\tau}\right)\right\\}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)\left\\{\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{\tau}\right)-\exp\left(-\frac{{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}}}{\tau}\right)\right\\}\Bigg{]}$
(63)
$\displaystyle=\frac{2\tau^{2}\left({\color[rgb]{0,0,1}{t_{i}^{(l)}}}-{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}}\right)}{V_{\text{th}}-\hat{v}}\sum_{j\in\Gamma_{i}^{(l)}}\frac{w_{ij}^{(l)}}{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}-{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}}}\Bigg{[}-2\exp\left(\frac{t_{j}^{(l-1)}}{2\tau}\right)\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{2\tau}\right)\left\\{1-\exp\left(\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}-{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}}}{2\tau}\right)\right\\}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{\tau}\right)\left\\{1-\exp\left(\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}-{\color[rgb]{0,0,1}{\hat{t}_{i}^{(l)}}}}{\tau}\right)\right\\}\Bigg{]}$
(64) $\displaystyle\overset{\hat{v}\rightarrow
V_{\text{th}}}{=}2\tau{\color[rgb]{0,0,1}{\alpha_{i}^{(l)}}}\left[\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{2\tau}\right)\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{2\tau}\right)-\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{\tau}\right)\sum_{j\in\Gamma_{i}^{(l)}}w_{ij}^{(l)}\exp\left(\frac{t_{j}^{(l-1)}}{\tau}\right)\right]$
(65)
$\displaystyle=2\tau{\color[rgb]{0,0,1}{\alpha_{i}^{(l)}}}\left[\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{2\tau}\right)b_{i}^{(l)}-\exp\left(-\frac{{\color[rgb]{0,0,1}{t_{i}^{(l)}}}}{\tau}\right)a_{i}^{(l)}\right].$
(66)

##  Appendix B Numerical convergence analysis

M-SSR was derived by taking the limit of $\hat{v}\rightarrow V_{\text{th}}$
for the integral-form regularization (Eq. 2). In this section, we confirm that
the gradient of the integral-form regularization (Eq. 2) converges to that of
M-SSR (Eq. 3) when $\hat{v}$ approaches $V_{\text{th}}$ through numerical
simulation. The Iris dataset [52] was used as input data. The Iris dataset
consists of 150 instances of three-class data with four features. To process
the data with an SNN, each feature was normalized to $[0,1]$ and each element
of the vector ($x_{i}$) was transformed to the time
$t_{i}^{(0)}=\tau_{\text{in}}x_{i}~{}(i=0,1,2,3)$ of the input layer spikes,
where $\tau_{\text{in}}=5$. We also introduced a bias spike $t_{4}^{(0)}=0$.
This input spike was input to an SNN with two hidden layers (5-10-10-3). The
cost function (Eq. (16)) was active only for $V$, with $\gamma_{2}=1$. The
error for the gradient of the weights for each layer was defined as follows:

$\displaystyle
Error^{(l)}(x)=\frac{1}{N^{(l)}N^{(l-1)}}\sum_{i=1}^{N^{(l)}}\sum_{j=1}^{N^{(l-1)}}\frac{\left|\frac{\partial
C(\hat{v}=x)}{\partial w_{ij}^{(l)}}-\frac{\partial C(\hat{v}\rightarrow
V_{\text{th}})}{\partial w_{ij}^{(l)}}\right|}{\left|\frac{\partial
C(\hat{v}=x)}{\partial w_{ij}^{(l)}}\right|}.$ (67)

The integral-form regularization (Eq. (2)) was computed by setting
$T=t^{\text{ref}}$ and dividing the integral by the time step width $\Delta
t=\frac{T}{N_{\text{steps}}}$.

Figure 7: Numerical conversion analysis. The gradient errors in the SNNs
(5-10-10-3) for various values of $\hat{v}$ and
$N_{\text{steps}}\in\\{10^{4},10^{5},10^{6},10^{7}\\}$ are plotted. The
results are shown for the following neuron models: (a)
$\tau_{v}=\tau_{I}=\infty$, (b) $\tau_{v}=\infty$, $\tau_{I}=5$, (c)
$\tau_{v}=2\tau_{I}=10$, and (d) $\tau_{v}=2\tau_{I}=20$. In each subfigure,
the upper panel shows the gradient error relating to the first hidden layer
and the lower panel shows the gradient error relating to the second hidden
layer. The gradient error was calculated by averaging over the 150 data. We
set $t^{\text{ref}}=10$.

Figure 7 shows the results of the gradient error for various neuron models.
The gradient error was calculated by averaging over the 150 data. For all
neuron models and $N_{\text{steps}}$, we can confirm that the errors decrease
as $\hat{v}$ is brought closer to $V_{\text{th}}$, and then the errors tend to
increase. This increase in error is due to the fact that if $\hat{v}$ is too
close to $V_{\text{th}}$, the integration range becomes extremely narrow,
making it difficult to evaluate the integral. As $N_{\text{steps}}$ is
increased, it is found that the integrals can be evaluated more precisely up
to the point where $\hat{v}$ is closer to $V_{\text{th}}$. As
$N_{\text{steps}}$ is increased, the minimum value of the gradient error
decreases uniformly, and it is found to decrease to about $0.1\%$ for
$N_{\text{steps}}=10^{7}$.

##  Appendix C Effects of $\xi$

Figure 8: Effect of $\xi$ on the sparsity–accuracy tradeoff. The
sparsity–accuracy tradeoffs are plotted for SNNs (784-400-400-400-10) with
M-SSR and F-SSR for various values of $\xi$. For each dataset, the first three
panels present the sparsity–accuracy tradeoff for the first, second, and third
hidden layers, from the top. The bottom panel presents the sparsity–accuracy
tradeoff for the sparsity averaged over the three hidden layers. We used the
following hyperparameters:
$t^{\text{ref}}=9,~{}\gamma_{1}=10^{-4},\eta=10^{-4}$, and
$\tau_{\text{soft}}=0.9$.

In the proposed SSR methods (Eqs. (1) and (7)), the regularization loss is
larger in the latter layers as the coefficient $\xi$ increases, resulting in
sparse firing characteristics in the latter layers. Figure 8 shows the effect
of the value of $\xi$ on the sparsity–accuracy tradeoff with M-SSR (Eq. (3))
and F-SSR (Eq. (8)). We trained the SNNs with three hidden layers
(784-400-400-400-10), and the results are shown for the MNIST and Fashion-
MNIST datasets. For each dataset, the top three panels show the
sparsity–accuracy tradeoffs for the first, second, and third hidden layers,
and the bottom panel shows the average results for all hidden layers. In the
case of M-SSR, when $p_{\text{layer}}=1$, the latter layer is not sparse, but
when $p_{\text{layer}}=6$, both layers are sparse. By contrast, in the case of
F-SSR, even when $p_{\text{layer}}=1$, the latter layers are relatively
sparse. This difference can be attributed to the fact that in M-SSR, the
regularization error propagates backward to the previous layer, whereas in
F-SSR, the regularization error is local and does not propagate backward.

##  Appendix D Raw plots of sparsity–accuracy tradeoff for CNN-SNNs

Figure 9: Raw plot of the sparsity–accuracy tradeoff on the CIFAR-10 dataset.

Figure 9 shows the sparsity–accuracy tradeoffs for SNNs trained on the
CIFAR-10 dataset. The regularization strength of F-SSR $\gamma_{2}$ ranged
from $0$ to $-10^{-6}$. The training was performed with 10 different initial
weights for each value of $\gamma_{2}$. All 190 data points obtained are shown
in the figure.

## References

  * [1] Kaushik Roy, Akhilesh Jaiswal, and Priyadarshini Panda.  Towards spike-based machine intelligence with neuromorphic computing.  Nature, 575:607–617, 2019.
  * [2] Wenzhe Guo, Mohammed E. Fouda, Ahmed M. Eltawil, and Khaled Nabil Salama.  Neural coding in spiking neural networks: A comparative study for robust neuromorphic systems.  Frontiers in Neuroscience, 15:638474, 2021.
  * [3] Daniel Auge, Julian Hille, Etienne Mueller, and Alois Knoll.  A survey of encoding techniques for signal processing in spiking neural networks.  Neural Processing Letters, 53:4693–4710, 2021.
  * [4] P. U. Diehl, D. Neil, J. Binas, M. Cook, S. Liu, and M. Pfeiffer.  Fast-classifying, high-accuracy spiking deep networks through weight and threshold balancing.  In 2015 International Joint Conference on Neural Networks (IJCNN), pages 1–8, July 2015.
  * [5] Bodo Rueckauer, Iulia-Alexandra Lungu, Yuhuang Hu, Michael Pfeiffer, and Shih-Chii Liu.  Conversion of continuous-valued deep networks to efficient event-driven networks for image classification.  Frontiers in Neuroscience, 11:682, 2017.
  * [6] Seijoon Kim, Seongsik Park, Byunggook Na, and Sungroh Yoon.  Spiking-YOLO: spiking neural network for energy-efficient object detection.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 11270–11277, 2020.
  * [7] M. Davies et al.  Loihi: A neuromorphic manycore processor with on-chip learning.  IEEE Micro, 38(1):82–99, January 2018.
  * [8] Alison L. Barth and James F.A. Poulet.  Experimental evidence for sparse firing in the neocortex.  Trends in Neurosciences, 35(6):345–355, 2012.
  * [9] Hiroshi Fujii, Hiroyuki Ito, Kazuyuki Aihara, Natsuhiro Ichinose, and Minoru Tsukada.  Dynamical cell assembly hypothesis ― Theoretical possibility of spatio-temporal coding in the cortex.  Neural Networks, 9(8):1303–1350, 1996.
  * [10] Tim Gollisch and Markus Meister.  Rapid neural coding in the retina with relative spike latencies.  Science, 319(5866):1108–1111, 2008.
  * [11] Geoffrey Portelli, John M Barrett, Gerrit Hilgen, Timothée Masquelier, Alessandro Maccione, Stefano Di Marco, Luca Berdondini, Pierre Kornprobst, and Evelyne Sernagor.  Rank order coding: a retinal information decoding strategy revealed by large-scale multielectrode array retinal recordings.  eNeuro, 3(3):e0134, 2016.
  * [12] Phase precession: a neural code underlying episodic memory?  Current Opinion in Neurobiology, 43:130–138, 2017.
  * [13] Amirhossein Tavanaei, Masoud Ghodrati, Saeed Reza Kheradpisheh, TimothÃ©e Masquelier, and Anthony Maida.  Deep learning in spiking neural networks.  Neural Networks, 111:47–63, 2019.
  * [14] Michael Pfeiffer and Thomas Pfeil.  Deep learning with spiking neurons: Opportunities and challenges.  Frontiers in Neuroscience, 12(774):1–18, 2018.
  * [15] Manon Dampfhoffer, Thomas Mesquida, Alexandre Valentian, and Lorena Anghel.  Backpropagation-based learning techniques for deep spiking neural networks: A survey.  IEEE Transactions on Neural Networks and Learning Systems, pages 1–16, 2023.
  * [16] Jason K. Eshraghian, Max Ward, Emre Neftci, Xinxin Wang, Gregor Lenz, Girish Dwivedi, Mohammed Bennamoun, Doo Seok Jeong, and Wei D. Lu.  Training spiking neural networks using lessons from deep learning.  arXiv:2109.12894, 2023.
  * [17] Emre O. Neftci, Hesham Mostafa, and Friedemann Zenke.  Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks.  IEEE Signal Processing Magazine, 36(6):51–63, 2019.
  * [18] Emre O. Neftci, Charles Augustine, Somnath Paul, and Georgios Detorakis.  Event-driven random back-propagation: Enabling neuromorphic deep learning machines.  Frontiers in Neuroscience, 11:324, 2017.
  * [19] Dongsung Huh and Terrence J. Sejnowski.  Gradient descent for spiking neural networks.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 1440â–1450, 2018.
  * [20] Friedemann Zenke and Surya Ganguli.  SuperSpike: Supervised Learning in Multilayer Spiking Neural Networks.  Neural Computation, 30(6):1514–1541, 2018.
  * [21] Yujie Wu, Lei Deng, Guoqi Li, Jun Zhu, and Luping Shi.  Spatio-temporal backpropagation for training high-performance spiking neural networks.  Frontiers in Neuroscience, 12:331, 2018.
  * [22] Friedemann Zenke and Tim P. Vogels.  The Remarkable Robustness of Surrogate Gradient Learning for Instilling Complex Function in Spiking Neural Networks.  Neural Computation, 33(4):899–925, 2021.
  * [23] Bojian Yin, Federico Corradi, and Sander M Bohté.  Accurate and efficient time-domain classification with adaptive spiking recurrent neural networks.  Nature Machine Intelligence, 3(10):905–913, 2021.
  * [24] Guillaume Bellec, Franz Scherr, Anand Subramoney, Elias Hajek, Darjan Salaj, Robert Legenstein, and Wolfgang Maass.  A solution to the learning dilemma for recurrent networks of spiking neurons.  Nature communications, 11(1):1–15, 2020.
  * [25] Youngeun Kim and Priyadarshini Panda.  Revisiting batch normalization for training low-latency deep spiking neural networks from scratch.  Frontiers in Neuroscience, 15:773954, 2021.
  * [26] Benjamin Cramer, Sebastian Billaudelle, Simeon Kanya, Aron Leibfried, Andreas GrÃ¼bl, Vitali Karasenko, Christian Pehle, Korbinian Schreiber, Yannik Stradmann, Johannes Weis, Johannes Schemmel, and Friedemann Zenke.  Surrogate gradients for analog neuromorphic computing.  Proceedings of the National Academy of Sciences, 119(4):e2109194119, 2022.
  * [27] Yulong Yan, Haoming Chu, Yi Jin, Yuxiang Huan, Zhuo Zou, and Lirong Zheng.  Backpropagation with sparsity regularization for spiking neural network learning.  Frontiers in Neuroscience, 16:760298, 2022.
  * [28] Sander M. Bohte, Joost N. Kok, and Han La PoutrÃ©.  Error-backpropagation in temporally encoded networks of spiking neurons.  Neurocomputing, 48(1):17–37, 2002.
  * [29] Simon Thorpe, Arnaud Delorme, and Rufin Van Rullen.  Spike-based strategies for rapid processing.  Neural Networks, 14(6):715–725, 2001.
  * [30] Lina Bonilla, Jacques Gautrais, Simon Thorpe, and TimothÃ©e Masquelier.  Analyzing time-to-first-spike coding schemes: A theoretical approach.  Frontiers in Neuroscience, 16:971937, 2022.
  * [31] H. Mostafa.  Supervised learning based on temporal coding in spiking neural networks.  IEEE Transactions on Neural Networks and Learning Systems, 29(7):3227–3235, 2018.
  * [32] S. R. Kheradpisheh and TimothÃ©e Masquelier.  Temporal backpropagation for spiking neural networks with one spike per neuron.  International Journal of Neural Systems, 30(6):2050027, 2020.
  * [33] Iulia-Maria ComÅa, Krzysztof Potempa, Luca Versari, Thomas Fischbacher, Andrea Gesmundo, and Jyrki Alakuijala.  Temporal coding in spiking neural networks with alpha synaptic function: Learning with backpropagation.  IEEE Transactions on Neural Networks and Learning Systems, 33(10):5939–5952, 2022.
  * [34] Yusuke Sakemi, Kai Morino, Takashi Morie, and Kazuyuki Aihara.  A supervised learning algorithm for multilayer spiking neural networks based on temporal coding toward energy-efficient vlsi processor design.  IEEE Transactions on Neural Networks and Learning Systems, 34(1):394–408, 2023.
  * [35] Yusuke Sakemi, Takashi Morie, Takeo Hosomi, and Kazuyuki Aihara.  Effects of VLSI circuit constraints on temporal-coding multilayer spiking neural networks.  arXiv:2106.10382, 2021.
  * [36] Malu Zhang, Jiadong Wang, Jibin Wu, Ammar Belatreche, Burin Amornpaisannon, Zhixuan Zhang, Venkata Pavan Kumar Miriyala, Hong Qu, Yansong Chua, Trevor E. Carlson, and Haizhou Li.  Rectified linear postsynaptic potential function for backpropagation in deep spiking neural networks.  IEEE Transactions on Neural Networks and Learning Systems, 33(5):1947–1958, 2022.
  * [37] J Göltz et al.  Fast and energy-efficient neuromorphic deep learning with first-spike times.  Nature Machine Intelligence, 3:823–835, 2021.
  * [38] Seongbin Oh, Dongseok Kwon, Gyuho Yeom, Won-Mook Kang, Soochang Lee, Sung Yun Woo, Jaehyeon Kim, and Jong-Ho Lee.  Neuron circuits for low-power spiking neural networks using time-to-first-spike encoding.  IEEE Access, 10:24444–24455, 2022.
  * [39] M. G. Sarwar Murshed, Christopher Murphy, Daqing Hou, Nazar Khan, Ganesh Ananthanarayanan, and Faraz Hussain.  Machine learning at the network edge: A survey.  ACM Comput. Surv., 54(8):1–37, oct 2021.
  * [40] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [41] H. Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-MNIST: a novel image dataset for benchmarking machine learning algorithms.  arXiv:1708.07747, 2017.
  * [42] Alex Krizhevsky.  Learning multiple layers of features from tiny images.
  * [43] Shibo Zhou, Xiaohua Li, Ying Chen, Sanjeev T. Chandrasekaran, and Arindam Sanyal.  Temporal-coded deep spiking neural network with easy training and robust performance.  Proceedings of the AAAI Conference on Artificial Intelligence, 35(12):11143–11151, 2021.
  * [44] Kakei Yamamoto, Yusuke Sakemi, and Kazuyuki Aihara.  Timing-based backpropagation in spiking neural networks without single-spike restrictions.  arXiv:2211.16113, 2022.
  * [45] S. Moradi, N. Qiao, F. Stefanini, and G. Indiveri.  A scalable multicore architecture with heterogeneous memory structures for dynamic neuromorphic asynchronous processors (DYNAPs).  IEEE Transactions on Biomedical Circuits and Systems, 12(1):106–122, 2018.
  * [46] Yusuke Sakemi, Kai Morino, Takashi Morie, Takeo Hosomi, and Kazuyuki Aihara.  A spiking neural network with resistively coupled synapses using time-to-first-spike coding towards efficient charge-domain computing.  In 2022 IEEE International Symposium on Circuits and Systems (ISCAS), pages 2152–2156, 2022.
  * [47] Seongbin Oh, Dongseok Kwon, Gyuho Yeom, Won-Mook Kang, Soochang Lee, Sung Yun Woo, Jang Saeng Kim, Min Kyu Park, and Jong-Ho Lee.  Hardware implementation of spiking neural networks using time-to-first-spike encoding.  arXiv:2006.05033, 2020.
  * [48] Saeed Reza Kheradpisheh, Maryam Mirsadeghi, and Timothée Masquelier.  BS4NN: Binarized spiking neural networks with temporal coding and learning.  Neural Processing Letters, 54(2):1255–1273, 2022.
  * [49] Faramarz Faghihi, Hany Alashwal, and Ahmed A Moustafa.  A synaptic pruning-based spiking neural network for hand-written digits classification.  Frontiers in Artificial Intelligence, 5:680165, 2022.
  * [50] Bing Han, Feifei Zhao, Yi Zeng, and Wenxuan Pan.  Adaptive sparse structure development with pruning and regeneration for spiking neural networks.  arXiv preprint arXiv:2211.12219, 2022.
  * [51] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv, 1412.6980, 2014.
  * [52] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.