11institutetext: King’s College London, London, UK
11email<EMAIL_ADDRESS>11email:
<EMAIL_ADDRESS>11email<EMAIL_ADDRESS>22institutetext: Imperial College London, London, UK
22email<EMAIL_ADDRESS>
# Average Calibration Error: A Differentiable Loss for Improved Reliability in
Image Segmentation

Theodore Barfoot 11 Luis C. Garcia Peraza Herrera 11 Ben Glocker 22
Tom Vercauteren 11

###### Abstract

Deep neural networks for medical image segmentation often produce
overconfident results misaligned with empirical observations. Such
miscalibration, challenges their clinical translation. We propose to use
marginal L1 average calibration error (mL1-ACE) as a novel auxiliary loss
function to improve pixel-wise calibration without compromising segmentation
quality. We show that this loss, despite using hard binning, is directly
differentiable, bypassing the need for approximate but differentiable
surrogate or soft binning approaches. Our work also introduces the concept of
_dataset reliability histograms_ which generalises standard reliability
diagrams for refined visual assessment of calibration in semantic segmentation
aggregated at the dataset level. Using mL1-ACE, we reduce average and maximum
calibration error by 45% and 55% respectively, maintaining a Dice score of 87%
on the BraTS 2021 dataset. We share our code here:
https://github.com/cai4cai/ACE-DLIRIS.

##  1 Introduction

Deep neural networks (DNNs) have significantly advanced the field of semantic
segmentation. However, their large capacity, makes them susceptible to
overfitting, leading to overconfident predictions that do not accurately
reflect the underlying uncertainties inherent to the task [1]. In medical
image segmentation, where the confidence of predictions can be as crucial as
the predictions themselves, and overfitting is exacerbated by smaller
datasets, such overconfidence poses significant risks, particularly when this
miscalibration is not communicated to the end-user. Appropriate management of
confidence can also have benefits for patient outcomes, for example in
radiotherapy contouring [2].

The Dice Similarity Coefficient (DSC) loss, is a popular choice in medical
image segmentation due to its robustness against class imbalance [3], but it
is known to yield poorly calibrated, overly confident predictions [4].

Reliability diagrams are the reference means of assessing and visualising the
calibration performance of a predictor over a dataset and this has extensively
been used for classification tasks [5]. They are constructed by binning the
samples according to a discretisation of the predicted probabilities, also
referred to as confidence. Within each bin, the average confidence and average
accuracy is computed. The reliability diagram is obtained by plotting the
Difference between Confidence and Accuracy (DCA) as a function of confidence.
An example is shown in Figure 1. An aggregation of these DCA values can then
be used to estimate measures of calibration error. While reliability diagrams
and associated calibration error metrics are typically computed at the dataset
level, in semantic segmentation, these diagrams can also be constructed per
image, due to multiple pixels/voxels being present. Capturing both image-
specific and dataset-wide reliability is of particular relevance in medical
imaging, due to the safety critical nature of the task. Box-plot diagrams have
been used in combination with image-specific reliability diagrams to represent
calibration variability within an image segmentation dataset [6].

Predicted Foreground Probability / Confidence
---

Count
---

Empirical Foreground
---

Frequency / Accuracy
---

Figure 1: Image-specific reliability diagrams for whole tumour segmentation
(BraTS case 00054) under different training losses. Left: DSC; Middle: DSC +
temperature scaling (Ts); Right: DSC + mL1-ACE (proposed).

Traditional post-hoc calibration strategies, such as temperature scaling (Ts)
[1, 7], offer a partial remedy by adjusting confidence levels uniformly. While
more nuanced post-hoc calibration methods offer class-specific adjustments [8,
9, 10], they still fail to exploit the extensive parameter space of DNNs [11,
12].

Train-time calibration methods are more effective in improving calibration
than post-hoc methods [1]. They leverage the large model capacity in DNNs to
provide fine-grained calibration adjustment with no detriments for the primary
task. For image classification tasks, DCA-based auxiliary losses such as
expected calibration error (ECE) have been shown to improve calibration [11].
However, the accuracy term in the DCA-based losses is not differentiable due
to the argmax operator. Furthermore, the usage of hard binning in standard DCA
has led to previous work seeking for alternative differentiable calibration
losses [11]. Multi-class DCA (MDCA) avoids binning, and its associated
differentiability challenges, by calculating DCA averaged over each class
within a mini-batch rather than averaged within confidence bins [12]. A
differentiable surrogate for expected calibration error (DECE), has been
proposed [13]. However, the DECE approach relies on an approximation of
accuracy, as well as soft binning to avoid the non-differentiability of hard
binning. The DECE work also exploits non-trivial meta-learning rather than
direct optimisation of the loss thereby complicating its use and limiting its
wider applicability.

In a different stream of work, authors have introduced specialised loss
functions which do not rely on DCA methods to tackle overconfidence in
biomedical image segmentation. For example, DSC++ loss, which uses a modified
DSC loss to enhance calibration by adjusting the false positive (FP) and false
negative (FN) terms of DSC based on class frequency [14]. Neighbor-Aware
Calibration (NACL) is another example focusing on spatial consistency, via
penalty constraints on logits, to improve reliability and performance in
segmentation tasks [15, 16]. While these methods have demonstrated some
improvement in calibration, their indirect approach make them more difficult
to interpret and potentially lead to suboptimal results.

In this work, we present a novel auxiliary loss function, marginal L1 Average
Calibration Error (mL1-ACE), designed to enhance calibration accuracy in
medical image segmentation without compromising on segmentation quality. This
represents the first clear attempt to use established calibration metrics
directly as an auxiliary loss in the context of medical image segmentation.
Despite presumptions related in previous DCA-based methods, our approach is
inherently differentiable, even with hard-binning of probabilities,
eliminating the need for surrogates or soft-binning techniques. mL1-ACE
employs a class-wise calibration strategy, opting for ACE over ECE to better
address the class imbalances in medical datasets. When combined with standard
loss functions, in particular DSC, mL1-ACE demonstrates superior calibration
performance on the BraTS2021 dataset [17], significantly surpassing
conventional post-hoc calibration methods like Ts. Additionally, we introduce
dataset reliability histograms, offering a comprehensive visual tool for
evaluating and communicating model calibration across image segmentation
datasets, enhancing transparency in reporting model reliability beyond
predictive performance.

##  2 Method

For calculation of image-specific calibration metrics, and our subsequent
mL1-ACE loss, we discretise the continuous predicted probabilities of each of
the $C$ classes into $M$ bins. For any given image, and class $c$, the
expected foreground probability $e^{c}_{m}$ within bin $m$ (also referred to
as confidence in DCA) is compared to the observed frequency $o^{c}_{m}$ within
that same bin (also referred to as accuracy in DCA). This provides us with an
analogous representation of DCA but extended to multi-class setting. An
ideally calibrated model is one where the gaps, $|o^{c}_{m}-e^{c}_{m}|$ equal
zero over all $M$ bins and all $C$ classes. That is, a perfectly calibrated
model is one where the predicted foreground probabilities align precisely with
the observed frequency of that class across numerous samples. For example, if
a well-calibrated brain tumour segmentation model segments $100$ voxels, each
with a probability/confidence of $0.7$, then we expect $70$ of those voxels to
be tumour.

Expected Calibration Error (ECE) [18] distils the reliability diagram
described above into a single scalar value that can be used to assess the
calibration of the model, by considering the weighted sum of the absolute
difference of the gaps over all the bins. We extend this approach to a multi-
class and image-specific variant, which we refer to as marginal L1 expected
calibration error:

$\textrm{mL1-ECE}=\frac{1}{C}\sum_{c}^{C}\sum_{m}^{M}\frac{n^{c}_{m}}{N}\lvert
o^{c}_{m}-e^{c}_{m}\rvert$ (1)

where $n^{c}_{m}$ is the number of in bin $B^{c}_{m}$, the $m^{\textrm{th}}$
bin for class $c$, and $N$ is the total number of voxels for an image. More
formally, we define $o^{c}_{m}$ and $e^{c}_{m}$ as:

$o^{c}_{m}=\operatorname{\mathbb{P}}[\bm{Y}^{c}_{i}\,\rvert\,i\in
B^{c}_{m}],\quad
e^{c}_{m}=\operatorname{\mathbb{E}}[f_{\bm{\theta}}(\bm{X})^{c}_{i}\,\rvert\,i\in
B^{c}_{m}]$ (2)

where $\bm{X}$ is the input image, $\bm{Y}$ the associated ground truth label-
map, and $f_{\bm{\theta}}$ is our predictor with trainable weights $\theta$.

Average Calibration Error (ACE) [19], performs this distillation in a similar
way, but considers an unweighted sum over all bins. We extend this approach to
a multi-class and image-specific variant, which we refer to as marginal L1
average calibration error:

$\textrm{mL1-ACE}=\frac{1}{CM}\sum_{c}^{C}\sum_{m}^{M}\lvert
o^{c}_{m}-e^{c}_{m}\rvert$ (3)

Finally, we report Maximum Calibration Error (MCE) [18] which is only
concerned with the maximum gap but extend it to a multi-class and image-
specific variant while referring to it as marginal L1 maximum calibration
error:

$\textrm{mL1-MCE}=\frac{1}{C}\sum_{c}^{C}\max_{m}\lvert
o^{c}_{m}-e^{c}_{m}\rvert$ (4)

We anticipate that all three measures, as shown in equations (1), (3), and
(4), should be differentiable almost everywhere. Combined with their image-
specific nature, this characteristic makes them trivially usable as auxiliary
losses that can be paired with any image segmentation loss. Specifically,
useful differentiability arises due to the large number of voxels that have
membership to each bin with averaging being performed within each bin rather
than across the bins themselves. However, given that mL1-MCE only considers
the maximum gap, for a single bin, it is less likely to provide a sufficient
supervisory signal compared to mL1-ECE and mL1-ACE.

Since mL1-ECE is weighted by the relative proportion of voxels in each bin, if
a model is performing sufficiently well, mL1-ECE is heavily skewed towards the
first and last bin as is evident in Figure 1. Therefore, we argue that mL1-ACE
is the most appropriate measure for calibration, and therefore an excellent
candidate for an auxiliary loss targeted at model calibration, in the context
of medical image segmentation. Indeed, mL1-ACE will be much more sensitive to
uncertain predictions such as those found on the periphery of segmented
regions.

Predicted Foreground Probability / Confidence
---

Empirical Foreground Frequency / Accuracy
---

CE
---
DSC
---
DSC + CE
---
\+ mL1-ACE
---
\+ mL1-ACE
---
\+ mL1-ACE
---
Figure 2: Dataset reliability histograms for BraTS 2021 dataset, with
different loss functions, without (top row) and with (bottom row) mL1-ACE.

To capture reliability information at the dataset-level rather than image-
level, we suggest to aggregate the image-specific metrics rather than
computing a dataset-level DCA. This not only provides an opportunity to gain
insights in the variability within the dataset but also offers a pragmatic
solution to deal with large datasets. To support visual assessment of dataset-
level reliability, we propose an extension of reliability diagrams illustrated
in Figure 2 which we refer to as dataset reliability histograms. By binning
the empirical foreground frequencies within each predicted foreground
probability bin across all cases in a dataset, we create a joint histogram
representation of reliability. This provides an informative visualisation of
calibration across an entire dataset.

##  3 Experiments

Training, validation and testing was performed using the Multimodal Brain
Tumour Segmentation (BraTS) Benchmark training dataset from 2021 [17].
Segmentations consist of three tumour sub-regions. These sub-regions then form
clinical tumour sub-regions, summarised as follows: Enhancing tumour (ET)
captures the gadolinium enhancement. Tumour core (TC) is comprised of the
necrotic tumour core (NCR) and ET. Whole tumour (WT) is comprised of the TC
and the peritumoral edematous and infiltrated tissue (ED).

The training dataset of $1251$ cases was randomly split into training,
validation and testing with a ratio of $1000$:$51$:$200$ cases. A basic UNet
architecture [20] implemented in MONAI [21] was used with $8$, $16$, $32$, and
$64$ channels and strides of 2. We use the same data augmentation pipeline,
for training, as in [22], closely matching the one used for the nnUNet [23].

A training patch size of $224\times 224\times 144$ was used, with a batch size
of $8$. Training was ran for $1000$ epochs, with $16$ iterations per epoch.
The final model for evaluation was the one with the best validation DSC.
Different training strategies were compared, using combinations of popular
segmentation losses with and without the mL1-ACE loss. We considered cross-
entropy (CE), DSC, and CE + DSC loss functions. Equal weighting for each loss
function component was used. For all calibration metrics and mL1-ACE loss,
$20$ bins were used to discretise the predicted probability space. For all
runs, an Adam optimiser with a learning rate of $0.001$ was used. Post-hoc Ts
was performed for each model, using CE loss function, with the temperature
parameter optimised using the validation dataset. All experiments were
performed using an NVidia DGX cluster with NVidia V100 (32GB) and Nvidia A100
(40GB) GPUs.

##  4 Results

We first consider the impact of our auxiliary loss in terms of segmentation
performance and then focus on calibration. The DSC scores for the different
loss functions are summarised in Table 1, showing that despite using a basic
network with no hyper-parameter tuning, we obtain competitive DSC values on
BraTS 2021. The addition of mL1-ACE as an auxiliary loss with CE or DSC-CE
losses, results in a small reduction in DSC scores (p < 0.01) However, for DSC
loss, the addition of mL1-ACE maintains DSC performance, showing no
statistically significant difference (p = 0.61). The DSC values of models with
post-hoc Ts are not included as Ts has no practically significant effect on
DSC.

1 - ACE
---
1 - MCE
---
1 - ECE
---

DSC
---

Figure 3: Comparative analysis of DSC against calibration metrics for average tumour component of BraTS 2021 dataset. Statistically significant difference (p < 0.01). Table 1: Dice Score across Loss Functions, higher is better Loss Function |  ET |  TC |  WT |  Avg
---|---|---|---|---
CE |  0.82 ± 0.19 |  0.87 ± 0.18 |  0.90 ± 0.11 |  0.86 ± 0.14
CE + L1-ACE |  0.80 ± 0.20 |  0.86 ± 0.19 |  0.87 ± 0.12 |  0.84 ± 0.15
DSC |  0.84 ± 0.17 |  0.87 ± 0.20 |  0.89 ± 0.12 |  0.87 ± 0.14
DSC + L1-ACE |  0.83 ± 0.18 |  0.88 ± 0.18 |  0.89 ± 0.11 |  0.87 ± 0.13
DSC-CE |  0.85 ± 0.17 |  0.89 ± 0.17 |  0.90 ± 0.11 |  0.88 ± 0.13
DSC-CE + L1-ACE |  0.82 ± 0.20 |  0.86 ± 0.21 |  0.88 ± 0.13 |  0.85 ± 0.16

We show the effect of Ts and mL1-ACE on the calibration of the models in Table
2. Our auxiliary loss has little effect on the calibration of the model
trained with CE. Our loss shows the biggest improvement in calibration when
paired with the DSC loss function, with a 45% reduction in ACE (p < 0.01).
Overall the use of Ts shows marginal improvement in calibration, with the
exception of CE loss where it has a small detrimental effect. Similar
conclusion can be drawn for MCE and ECE metrics, shown in Table A1 and Table
A2 in the appendix. We also note that the number of bins used to calculate our
calibration metrics, has little effect, with ACE increasing by less than 3%
when increasing the bins from 10 to 100 (p < 0.01), as shown in Figure A1 in
the appendix.

Table 2: Average Calibration Error across Loss Functions, lower is better Loss Function |  ET |  TC |  WT |  Avg
---|---|---|---|---
CE |  0.13 ± 0.09 |  0.14 ± 0.10 |  0.12 ± 0.07 |  0.13 ± 0.06
CE + Ts |  0.15 ± 0.07 |  0.13 ± 0.08 |  0.13 ± 0.06 |  0.14 ± 0.05
CE + L1-ACE |  0.12 ± 0.10 |  0.14 ± 0.10 |  0.11 ± 0.09 |  0.12 ± 0.07
CE + L1-ACE + Ts |  0.12 ± 0.07 |  0.13 ± 0.08 |  0.11 ± 0.07 |  0.12 ± 0.05
DSC |  0.25 ± 0.07 |  0.26 ± 0.07 |  0.23 ± 0.06 |  0.25 ± 0.05
DSC + Ts |  0.24 ± 0.07 |  0.24 ± 0.07 |  0.22 ± 0.06 |  0.24 ± 0.05
DSC + L1-ACE |  0.14 ± 0.10 |  0.15 ± 0.09 |  0.12 ± 0.08 |  0.14 ± 0.07
DSC + L1-ACE + Ts |  0.12 ± 0.09 |  0.13 ± 0.09 |  0.11 ± 0.07 |  0.12 ± 0.06
DSC-CE |  0.20 ± 0.08 |  0.21 ± 0.08 |  0.18 ± 0.09 |  0.20 ± 0.06
DSC-CE + Ts |  0.18 ± 0.09 |  0.18 ± 0.09 |  0.15 ± 0.09 |  0.17 ± 0.07
DSC-CE + L1-ACE |  0.15 ± 0.11 |  0.15 ± 0.11 |  0.11 ± 0.08 |  0.14 ± 0.08
DSC-CE + L1-ACE + Ts |  0.13 ± 0.10 |  0.14 ± 0.10 |  0.11 ± 0.06 |  0.13 ± 0.07

The relationship between DSC and calibration metrics is shown in Figure 3, for
the DSC loss model, showing that the addition of mL1-ACE as an auxiliary loss
maintains DSC performance with a substantial improvement in MCE and ACE (p <
0.01).

Figure 1 shows a reliability diagram for the DSC loss model, with and without
mL1-ACE and Ts, for the WT segmentation component. As can be seen, Ts has a
limited effect on the calibration for an individual case, whereas the addition
of mL1-ACE has a significant effect, clearly improving calibration.

Similar, but more pronounced, results are shown in Figure 2, where the
addition of mL1-ACE has a significant qualitative improvement on calibration
of the model across the entire dataset, for all loss functions, showing a
marked reduction in under- and over-confident cases, by the narrowing of the
distribution along the diagonal. We supplement this visual qualitative
analysis, by showing an individual case in Figure 4, where the addition of
mL1-ACE has a significant effect on the calibration of the model, with a
marked reduction in the confidence of false positive and false negative
predictions.

Figure 4: Segmentation of BraTS 2021 case 00014, with DSC loss function. Left:
Ground truth WT on FLAIR image, Middle: Predicted WT foreground probabilities,
Right: WT foreground probabilities with mL1-ACE. True positive, false
positive, and false negative contours are shown in blue, yellow, and red,
respectively.

##  5 Conclusion and discussion

Prior work assumed that hard binning leads to non-differentiable DCA-based
metrics. We demonstrate that it is sufficiently differentiable and provide a
multi-class and image-specific variant for use in challenging image
segmentation tasks. By performing averaging inside of each bin across a large
enough number of voxels our approach indeed provides a sufficiently
informative gradient updates. We further show that the introduction of the
proposed mL1-ACE as an equally weighted auxiliary loss function to DSC loss,
can improve the calibration of the model, with a $45$% reduction in ACE and
$55$% reduction in MCE metrics, without any significant change in the Dice
coefficient metric, remaining at competitive $87$%, when evaluated on the
BraTS 2021 dataset.

We also introduced the concept of dataset reliability histograms to visualise
the reliability of a semantic segmentation models across an entire dataset,
providing valuable insights into modes of failure.

Additionally, we show that models trained with cross-entropy loss tend to have
much better calibration than those trained with a DSC-based loss, as has been
reported in prior literature [14]. Our method has limited benefit for CE-
trained models. However, with DSC-based losses outperforming CE over a large
set of tasks [3], the need for well calibrated DSC-based losses remains. We
also show that the post-hoc calibration method of temperature scaling has
limited effect on the calibration of the model.

While our experiments were focused on BraTS with a baseline U-Net
architecture, we believe that our mL1-ACE auxiliary loss has the potential to
offer an easy-to-use and transferable solution to improve model calibration in
any medical image segmentation approach. Exploring the use with other
datasets, base losses, and architectures, including transformers, is left for
future work.

## Acknowledgments

The authors would like to acknowledge support for TB from EPSRC CDT
[EP/S022104/1]. and Intel, UK. TV is supported by a Medtronic / RAEng Research
Chair [RCSRF1819$\backslash$7$\backslash$34]. This work was supported by core
funding from the Wellcome/EPSRC [WT203148/Z/16/Z; NS/A000049/1]. For the
purpose of open access, the authors have applied a CC BY public copyright
licence to any Author Accepted Manuscript version arising from this
submission. TV is co-founder and shareholder of Hypervision Surgical. This
article only uses publicly available datasets. Their re-use did not require
any ethical approval.

## References

  * [1] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q. Weinberger.  On calibration of modern neural networks.  In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 1321–1330. PMLR, August 2017.
  * [2] David Bernstein, Alexandra Taylor, Simeon Nill, and Uwe Oelfke.  New target volume delineation and PTV strategies to further personalise radiotherapy.  Physics in Medicine & Biology, 66(5):055024, February 2021.
  * [3] Jun Ma, Jianan Chen, Matthew Ng, Rui Huang, Yu Li, Chen Li, Xiaoping Yang, and Anne L. Martel.  Loss odyssey in medical image segmentation.  Medical Image Analysis, 71:102035, 2021.
  * [4] Alireza Mehrtash, William M. Wells, Clare M. Tempany, Purang Abolmaesumi, and Tina Kapur.  Confidence calibration and predictive uncertainty estimation for deep medical image segmentation.  IEEE Transactions on Medical Imaging, 39(12):3868–3878, 2020.
  * [5] Morris H. DeGroot and Stephen E. Fienberg.  The comparison and evaluation of forecasters.  Journal of the Royal Statistical Society. Series D (The Statistician), 32(1/2):12–22, 1983.
  * [6] Farina Kock, Felix Thielke, Grzegorz Chlebus, and Hans Meine.  Confidence histograms for model reliability analysis and temperature calibration.  In Ender Konukoglu, Bjoern Menze, Archana Venkataraman, Christian Baumgartner, Qi Dou, and Shadi Albarqouni, editors, Proceedings of The 5th International Conference on Medical Imaging with Deep Learning, volume 172 of Proceedings of Machine Learning Research, pages 741–759. PMLR, July 2022.
  * [7] John Platt.  Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods.  In Advances in Large Margin Classifiers, 1999.
  * [8] Meelis Kull, Miquel Perello Nieto, Markus Kängsepp, Telmo Silva Filho, Hao Song, and Peter Flach.  Beyond temperature scaling: Obtaining well-calibrated multi-class probabilities with Dirichlet calibration.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [9] Zhipeng Ding, Xu Han, Peirong Liu, and Marc Niethammer.  Local temperature scaling for probability calibration.  In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pages 6889–6899, October 2021.
  * [10] Mobarakol Islam, Lalithkumar Seenivasan, Hongliang Ren, and Ben Glocker.  Class-distribution-aware calibration for long-tailed visual recognition.  In International Conference on Machine Learning, 2021.
  * [11] Gongbo Liang, Yu Zhang, Xiaoqin Wang, and Nathan Jacobs.  Improved trainable calibration method for neural networks on medical imaging classification.  In British Machine Vision Conference (BMVC), 2020.
  * [12] Ramya Hebbalaguppe, Surya Prakash, et al.  A stitch in time saves nine: A train-time regularizing loss for improved neural network calibration.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2022.
  * [13] Ondrej Bohdal, Yongxin Yang, and Timothy Hospedales.  Meta-calibration: Learning of model calibration using differentiable expected calibration error.  Transactions on Machine Learning Research, 2023.
  * [14] Michael Yeung, L. Rundo, Yang Nan, E. Sala, C. Schönlieb, and Guang Yang.  Calibrating the Dice loss to handle neural network overconfidence for biomedical image segmentation.  Journal of Digital Imaging, 2021.
  * [15] Balamurali Murugesan, Sukesh Adiga Vasudeva, Bingyuan Liu, Herve Lombaert, Ismail Ben Ayed, and Jose Dolz.  Trust your neighbours: Penalty-based constraints for model calibration.  In Hayit Greenspan, Anant Madabhushi, Parvin Mousavi, Septimiu Salcudean, James Duncan, Tanveer Syeda-Mahmood, and Russell Taylor, editors, Medical Image Computing and Computer Assisted Intervention – MICCAI 2023, pages 572–581, Cham, 2023. Springer Nature Switzerland.
  * [16] Balamurali Murugesan, Sukesh Adiga Vasudeva, Bingyuan Liu, H. Lombaert, Ismail Ben Ayed, and J. Dolz.  Neighbor-aware calibration of segmentation networks with penalty-based constraints.  arXiv preprint arXiv:2401.14487, 2024.
  * [17] Ujjwal Baid, Satyam Ghodasara, and Suyash Mohan et al.  The RSNA-ASNR-MICCAI BraTS 2021 benchmark on brain tumor segmentation and radiogenomic classification.  arXiv, July 2021.
  * [18] Mahdi Pakdaman Naeini, Gregory F. Cooper, and Milos Hauskrecht.  Obtaining well calibrated probabilities using bayesian binning.  In Blai Bonet and Sven Koenig, editors, Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence, January 25-30, 2015, Austin, Texas, USA, pages 2901–2907. AAAI Press, 2015.
  * [19] Lukas Neumann, Andrew Zisserman, and Andrea Vedaldi.  Relaxed softmax: efficient confidence auto-calibration for safe pedestrian detection.  In 2018 NIPS MLITS Workshop: Machine Learning for Intelligent Transportation System. OpenReview, 2018.
  * [20] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, pages 234–241. Springer, 2015.
  * [21] M. Jorge Cardoso, Wenqi Li, Richard Brown, et al.  MONAI: An open-source framework for deep learning in healthcare.  arXiv preprint arXiv:2211.02701, 2022.
  * [22] Lucas Fidon, Suprosanna Shit, Ivan Ezhov, Johannes C. Paetzold, Sébastien Ourselin, and Tom Vercauteren.  Generalized Wasserstein Dice loss, test-time augmentation, and transformers for the BraTS 2021 challenge.  In Alessandro Crimi and Spyridon Bakas, editors, Brainlesion: Glioma, Multiple Sclerosis, Stroke and Traumatic Brain Injuries, pages 187–196, Cham, 2022. Springer International Publishing.
  * [23] Fabian Isensee, Paul F. Jaeger, Simon A. A. Kohl, Jens Petersen, and Klaus H. Maier-Hein.  nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation.  Nat Methods, 18(2):203–211, February 2021.

## Appendix

Table A1: Maximum Calibration Error |  ET |  TC |  WT |  Avg
---|---|---|---|---
CE |  0.26 ± 0.19 |  0.28 ± 0.18 |  0.23 ± 0.12 |  0.25 ± 0.12
CE + Ts |  0.30 ± 0.16 |  0.29 ± 0.16 |  0.25 ± 0.10 |  0.28 ± 0.10
CE + L1-ACE |  0.24 ± 0.19 |  0.26 ± 0.18 |  0.20 ± 0.15 |  0.23 ± 0.14
CE + L1-ACE + Ts |  0.23 ± 0.15 |  0.25 ± 0.16 |  0.20 ± 0.12 |  0.23 ± 0.11
Dice |  0.59 ± 0.15 |  0.61 ± 0.15 |  0.55 ± 0.13 |  0.58 ± 0.11
Dice + Ts |  0.56 ± 0.15 |  0.58 ± 0.15 |  0.51 ± 0.14 |  0.55 ± 0.11
Dice + L1-ACE |  0.29 ± 0.18 |  0.30 ± 0.17 |  0.23 ± 0.14 |  0.27 ± 0.13
Dice + L1-ACE + Ts |  0.23 ± 0.17 |  0.25 ± 0.16 |  0.21 ± 0.11 |  0.23 ± 0.12
Dice-CE |  0.46 ± 0.17 |  0.46 ± 0.17 |  0.36 ± 0.17 |  0.43 ± 0.13
Dice-CE + Ts |  0.26 ± 0.19 |  0.32 ± 0.19 |  0.19 ± 0.12 |  0.26 ± 0.13
Dice-CE + L1-ACE |  0.29 ± 0.21 |  0.30 ± 0.20 |  0.21 ± 0.14 |  0.27 ± 0.15
Dice-CE + L1-ACE + Ts |  0.25 ± 0.19 |  0.27 ± 0.18 |  0.20 ± 0.11 |  0.24 ± 0.13
Table A2: Expected Calibration Error ($\times 10^{-4}$) |  ET |  TC |  WT |  Avg
---|---|---|---|---
CE |  4.21 ± 7.49 |  5.66 ± 10.5 |  9.71 ± 11.3 |  6.53 ± 7.53
CE + Ts |  7.77 ± 7.16 |  10.2 ± 10.0 |  19.1 ± 9.45 |  12.4 ± 6.93
CE + L1-ACE |  4.00 ± 4.75 |  6.16 ± 9.41 |  12.0 ± 12.9 |  7.38 ± 6.94
CE + L1-ACE + Ts |  21.1 ± 4.69 |  27.2 ± 6.65 |  62.0 ± 11.1 |  36.8 ± 5.69
Dice |  5.09 ± 7.67 |  7.48 ± 16.5 |  15.5 ± 17.8 |  9.34 ± 11.1
Dice + Ts |  4.95 ± 7.63 |  7.36 ± 16.5 |  15.0 ± 17.8 |  9.10 ± 11.1
Dice + L1-ACE |  3.47 ± 7.56 |  4.50 ± 8.45 |  10.1 ± 12.8 |  6.03 ± 6.84
Dice + L1-ACE + Ts |  3.27 ± 7.17 |  4.41 ± 8.00 |  12.0 ± 10.9 |  6.54 ± 5.84
Dice-CE |  4.40 ± 7.57 |  5.70 ± 12.2 |  11.6 ± 14.4 |  7.21 ± 8.62
Dice-CE + Ts |  3.84 ± 7.43 |  5.20 ± 11.8 |  10.2 ± 13.4 |  6.41 ± 8.11
Dice-CE + L1-ACE |  3.91 ± 8.09 |  6.04 ± 14.6 |  10.3 ± 14.8 |  6.75 ± 9.84
Dice-CE + L1-ACE + Ts |  4.23 ± 7.68 |  6.56 ± 14.0 |  13.9 ± 12.6 |  8.23 ± 8.83
Number of Bins
---

Metric Value
---

Figure A1: Comparison of calibration metrics using different number of bins
when evaluated on the DSC + mL1-ACE loss model. Both ACE and ECE are very
stable with respect to number of bins (5, 10, 20, 50, 100, 1000). MCE shows a
larger increase from 100 to 1000 bins, this is due to the presence of empty
bins, a likely occurrence when using such a high number of bins.


# Statistical complexity and connectivity relationship in cultured neural
networks

A. Tlaie<EMAIL_ADDRESS>L.M. Ballesteros-Esteban  I. Leyva  I. Sendiña-
Nadal  Complex Systems Group & GISC, Universidad Rey Juan Carlos, 28933
Móstoles, Madrid, Spain  Center for Biomedical Technology, Universidad
Politécnica de Madrid, 28223 Pozuelo de Alarcón, Madrid, Spain

###### Abstract

We explore the interplay between the topological relevance of a neuron and its
dynamical traces in experimental cultured neuronal networks. We monitor the
growth and development of these networks to characterise the evolution of
their connectivity. Then, we explore the structure-dynamics relationship by
simulating a biophysically plausible dynamical model on top of each networks’
nodes. In the weakly coupling regime, the statistical complexity of each
single node dynamics is found to be anti-correlated with their degree
centrality, with nodes of higher degree displaying lower complexity levels.
Our results imply that it is possible to infer the degree distribution of the
network connectivity only from individual dynamical measurements.

###### keywords:

Complex Networks, Neuron Models, Network Inference, Cultured Networks

††journal: Journal of LaTeX Templates

##  1 Introduction

One of the main research lines in the study of the dynamics of complex
networks has been the deep relationship between the connectivity and the
dynamics of the nodes, and how this interaction shapes the emergence of a
collective state such as synchronisation [1, 2, 3, 4]. An enormous effort has
been devoted to the understanding of this phenomenon, and the knowledge
gathered so far has driven the advance in crucial applications in brain
dynamics [5], power grids [6], and in many other fields where synchronisation
is essential [7, 8] for the system’s proper functioning.

Commonly, studies have focused in states of full synchronisation [4].
Nevertheless, there are very relevant cases in which only a partial or weak
synchronisation level is achieved [9, 10, 11], and often this state becomes
optimal to balance functional integration and segregation in the system [12,
13, 14] while a complete coordination is evidencing the existence of a
pathological condition.

Several investigations [15, 16, 17] have shown that nodes play different roles
in the ensemble dynamics depending on their topological position and intrinsic
dynamics [18]. One of the most explored situation is that of the hubs acting
as coordinators of the dynamics of the whole system [19, 20, 21, 22], being
the first nodes to synchronise among them [23] and to the mean field [24],
while the rest of the nodes progressively locks the hubs dynamics.

This effect of the topology on the dynamics in the weakly synchronised regime
opens the question of whether it is possible to infer the network architecture
from statistical correlations among the coupled units [15, 25]. Currently, a
great amount of the research is being conducted in this sense. In particular,
the computational neuroscience field roots in the hypothesis that dynamical
correlations (which can be recorded in non-invasive ways) are greatly
constrained and induced by the anatomical structure of the brain [5]. From
these site-to-site correlation maps, the functional brain networks, it is
often possible to obtain information about the underlying topological networks
[26].

However, it has been less explored the fact that this structural-dynamical
interaction also plays in the other way around: just as the dynamics of each
node influences the ensemble, the ensemble imprints its structural marks into
the dynamics of each individual node [23, 24]. We make the assumption that,
long before the coupling strength is high enough to induce synchronisation,
the dynamical changes at the node level are encoding the imprint of its
structural role. This relevant feature could be used to extract information
about the network without making any reference to pairwise correlations,
particularly in those cases where the structure is unknown or unreliable, as
we showed in a previous work [27].

Here, we extend our study of the influence that the ensemble has over the node
dynamics to an experimental case. We culture networks of neurons coming from
Schistocerca gregaria and study the potential relationship between a simulated
dynamical model (the Morris-Lecar neuron) and the anatomical network structure
of a neuronal culture. The main motivation for this study is that in cultured
neuronal networks the simultaneous obtention of structural and dynamical
information is not possible, either because one recording technique influences
the other measurement or mainly because the culture is not able to survive to
both measurements.

##  2 Experimental setup: culturing the network

To analyse the spatial structure of a real network, we focus on the study of
cultured neuronal networks (CNNs), considered as a simplified version of a
more complex network of the central nervous system [28, 29, 30]. For this
purpose , we analyse the network structure in this CNN model by means of
optical microscopy techniques, extracting the detailed connectivity and its
statistical topological properties.

Our CNNs were obtained from Schistocerca gregaria specimens, also known as
desert locusts. As they share basic neuronal features with vertebrates, this
invertebrate model has been recently used in neuroscience as an easier
approach for the understanding of more complex neural systems [31]. The large
size of its neurons makes it ideal for observing the structure of the network,
as an alternative to the mammals.

In our experiments we follow the protocol described in [32]. Each locust is
dissected to extract its frontal ganglion, formed by approximately 100 neurons
[33]. To obtain an intermediate neuron density that allows us to study a
complex network morphology we extracted $12$ ganglia per culture. After the
dissection, the frontal ganglia endure a chemical and mechanical procedure to
remove all the connections and dissociate the neurons. The neuronal somata are
cultured in a Petri dish, in an enriched environment to allow the neurites to
regrowth and form a new connectivity network. The cultures are monitored in
vitro from day 0 (DIV0, DIV=days in vitro) to day 14 (DIV14). The data used in
this work correspond to $6$ cultures grown in the same conditions. We
inspected the morphological features of the cultured networks using a phase
contrast inverted microscope (Eclipse Ti-S, Nikon) with a $10$x air objective
(Achromat, ADL, NA $0.25$) and an automated motorised $XYZ$ stage controller.
High resolution images were obtained in a daily basis.

In order to analyse the spatial network, we need to extract the corresponding
mathematical graph. To do so, we process the culture images by means of an
image segmentation algorithm [34, 35]. In Fig. 1 we portray the whole process,
starting from a typical microscope image of the culture (Fig. 1(a) shows just
a small area) and ending up with the output of the segmentation detecting
neurons (and aggregates of neurons) and the neuronal processess connecting
them Fig. 1(a). The algorithm is summarised in the following steps:

Figure 1: Image segmentation processing steps and extraction of the network
graph. (a) Red layer of a RGB cut of a culture $6$ DIV old. (b) Output of the
segmentation algorithm of the region of interest. Single neurons and
aggregates of neurons are highlighted in red while the neurites are marked in
green. (c) Mapping of the segmentation objects in (b) into a full graph where
circles (neurons and neuronal clusters) and diamonds (branching and end points
of neurites, with no neurons) are the nodes, these being connected with green
lines, representing the neurites. (d) Graph representing the projection of the
full graph into the cluster graph with only neuronal clusters (red) and
neurites (green), where the links represent the existence of a path between
cluster nodes through junction nodes in the full graph in (c).

  1. 1.

The red layer of a RGB high-resolution image of the recorded cultured network
is processed (Fig.1(a)).

  2. 2.

The image is segmented and thresholded to separate background from foreground
areas. Then neurons and aggregates of neurons (red areas in Fig. 1(b)) and
neurites (green paths in Fig. 1(b)) are identified separately .

  3. 3.

Both neurons and neurites are connected and coded in the adjacency matrix
where single and clustered neurons are the nodes and neurites are the links
between them. Branching and end points of the neurites are also registered as
junction nodes in the graph, even when there is no neuron on them. This
provides a complete version of the graph that we call the full graph (Fig.
1(c)) with two types of nodes, those corresponding to neurons and the ones
denoting a branching point in the neuronal process path connecting two
neurons.

  4. 4.

The previous data is used to build a reduced version of the graph, where only
neurons (or neuronal clusters) are the nodes, and the links represent the
existence of a path between neuronal clusters, eventually through branching
points. In this graph, junction nodes have been removed and we observe a more
direct path between neuronal clusters, obtaining a simple version of the
matrix, the cluster graph (Fig. 1(d)).

We analyse the morphological and topological properties of the cultured
network using both full and cluster graphs, where the links are unweighted.
With the purpose of characterising the segregation and integration of the
cultured neuronal network along the experiment, in Fig. 2 we measure the
longitudinal progression of the averaged clustering coefficient ($C$) and
shortest path length $L$, normalised by the size of the largest connected
component $S_{1}$ [3], resulting $L/S_{1}$. The relationship between $C$ and
$L$ is often used as an indicator of the balance between the local and long-
distance connectivity in the network. These aforementioned parameters were
measured in both the full graph (with neuronal clusters and junction nodes)
and in the cluster graph (with only neuronal clusters nodes).

In the full graph (Fig. 2 (a)), the normalised shortest path $L/S_{1}$ shows a
high mean value at the early days of the culturing, where the connectivity is
still not fully developed. Between DIV3 and DIV6 there is a significant
decrease and showed no significant change thereafter, meaning a high
integration degree in the mature culture network. The clustering coefficient
was characterised by a very low mean value, showing a slight increase between
DIV3 and DIV6, when the network development occurs. The low mean values of $C$
are due to the fact that both neurons and branching points are considered
nodes meaning that the probability of forming triangles is reduced (see Fig.
1(b)).

In the case of the cluster graphs (Fig. 2 (d)) we observe a similar trend in
$L/S_{1}$ (see Fig. 2(b)) as the one described for the full graph. On the
contrary, $C$ exhibits higher mean values, with a more acute increase between
DIV 3 and DIV 6, that coincides with the most intense developmental phase. As
neuronal aggregates are the only nodes in this cluster graph and junctions are
not represented as nodes, the mean values of $C$ are more accurate with the
connected structure. After that point, in the mature neuronal network these
two parameters keep constant values [34, 35].

Figure 2: Results of the cultured network analysis. Mean value of the
clustering coefficient $C$ and mean path length $L$ normalised by the size of
the largest connected component $S_{1}$ (a) in a full graph, where clusters
and junctions are the nodes, and (b) in a cluster graph, where the only nodes
are clusters of neurons. Each point is the average over 6 experiments.

The analysis of these parameters in both types of graphs concludes with the
emergence of a mature cultured neuronal network from an initial random stage.
This evolved structure is characterised by high clustering coefficients and
low mean path values, indicating the presence of a mature network with high
segregation (favored by high clustering values) and integration (facilitated
by the existence of small shortest paths) levels. These are the
characteristics of a small world structure, where the high tendency to form
clusters of nodes in highly interconnected subgroups and short distance
between them contribute to an optimal functionality in the network [34].

We also analized the degree distribution $P(k)$ of these networks, being $k$
the number of links that each node has. In Fig. 3 we plot an example of the
cumulative degree distribution $P_{c}(k)$ of an in-vitro clustered network at
DIV7 (black squares), compared to equivalent simulated networks (same number
of nodes and links) obtained from usual generative models: random Erdös-Renyi
(ER, blue diamonds), scale-free obtained by Barabasi-Albert algorithm (SF, red
dots) and a spatial network with a distance-dependence linkage pattern
(spatial-ER green diamonds). As described in[34], the cultured networks belong
to the single-scale type as they show a well defined $k$. The study of
$P_{cum}(k)$ reveals a fast decay with a large number of nodes with similar
number of connections, and a few ones with a different and large node degrees.
As it can be seen, the experimental connectivity largely differs from the pure
random ER, showing instead shared features between SF and spatial networks.

Figure 3: Comparative plot of the cumulative degree distribution $P_{c}(k)$ of
an experimental clustered neuronal network at DIV7 (black squares) and
equivalent synthetic networks: random Erdös-Renyi (ER, blue diamonds), scale-
free obtained with the Barabási-Albert algorithm (SF, red dots) and a spatial
network where distance-dependence linkage pattern (spatial-ER green diamonds).
We see that the cultured network’s distribution shows a similar behavior, for
low degrees, as the SF one, and there is a region in which it is similar to
the spatial network (namely, the decay for large values of $k$ is almost
identical).

##  3 Dynamical model

Once we have extracted the connectivity of real neuronal cultures , we can
provide a dynamical behavior to their nodes, in order to enable the
exploration of the potential interplay between structure and dynamics. For
this study we implemented the bio-inspired Morris-Lecar (ML) model [36, 27],
whose equations describing the membrane potential behavior for each unit read
[37, 17]:

$\displaystyle C\dot{V_{i}}$ $\displaystyle=-{}\overbrace{g_{\rm
X}M_{\infty}(V_{i}-V_{\rm X})}^{\text{Ionic
channels}}+q\xi_{i}+\underbrace{\frac{\sigma}{K}\sum_{j}a_{ij}\overbrace{e^{-2(t-t_{j})}(V_{0}-V_{i})}^{\text{Synaptic
function}}}_{I_{i}}+I^{ext}_{i},$ (1) $\displaystyle\dot{W_{i}}$
$\displaystyle=\phi\,\tau_{W}(W_{\infty}-W_{i})$

where $V_{i}$ and $W_{i}$ are, respectively, the membrane potential and the
fraction of open $\rm K^{+}$ channels of the $i$th neuron and
$M_{\infty},W_{\infty}$, and $\tau_{W}$ are hyperbolic functions dependent on
$V_{i}$ and $\phi$ is a reference frequency. The parameters $g_{\rm X}$ and
$V_{\rm X}$ account for the electric conductance and equilibrium potentials of
the $\rm X=\\{K,Ca,\text{leaky}\\}$ channels. The external current
$I_{i}^{ext}=50.0$ mA is the same for all the neurons and is chosen such that
neurons are sub-threshold to neuronal firing which is induced by the white
Gaussian noise $q\xi_{i}$ of zero mean and intensity $q$. The coupling of the
neuron $i$th with the neuron ensemble is described by the injected synaptic
current $I_{i}$, given by the superposition of all the post-synaptic
potentials emitted by the neighbours of node $i$ in the past, being $t_{j}$
the time of the last spike of node $j$, and the corresponding element of the
adyacency matrix is $a_{ij}=1$ if there is a link between nodes $i,j$ and
$a_{ij}=0$ otherwise. The synaptic conductance $\sigma$, normalised by the
largest node degree present in the network $K$, plays the role of coupling
intensity.

Additionally, the channel voltage-dependent saturation values are given by the
following functions:

$\displaystyle M_{\infty}(V_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1+\tanh\left(\frac{V_{i}-V_{1}}{V_{2}}\right)\right],$
$\displaystyle W_{\infty}(V_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1+\tanh\left(\frac{V_{i}-V_{3}}{V_{4}}\right)\right],$
(2) $\displaystyle\tau_{W}(V_{i})$ $\displaystyle=$
$\displaystyle\cosh\left(\frac{V_{i}-V_{3}}{2V_{4}}\right).$

We chose the parameters such that the neurons in the simulations corresponded
to type II class excitability for the neuron dynamics, which means that a
discontinuous transition is found in the dependence of the spiking frequency
on the external current. The values for all the parameters can be found in
Refs. [27, 38].

We have to remark, however, that this model was originally conceived for a
single neuron and in this work we are dealing with aggregates of $~{}20$ of
them as our individual nodes. Although this quantity is not enough for
employing a neural mass model, it should be noted that, strictly, it is not a
single neuron either.

##  4 Statistical characterisation of the dynamics

With the purpose of providing a solid description of the system, we use two
different quantities: a global and a local one. The global measure is the
synchronisation level of the network, i.e. how similar are the dynamical
outputs of our units, while the local one is an individual measure of the
complexity of a node’s time series.

###  4.1 Synchronization measure

In order to quantify the level of synchronization we estimate how many neurons
fire within the same time window. The total simulation time $T$ is divided in
$n=1,\dots,N_{b}$ bins of a convenient size $\tau$, such that $T=N_{b}\tau$,
and the binary quantity $B_{i}(n)$ is defined such that $B_{i}(n)=1$ if the
$i$th neuron spiked within $n$th interval and $0$ otherwise. The
synchronisation between the spiking sequences of neurons $i$ and $j$ is
therefore characterised with the pairwise correlation matrix $s_{ij}\in[0,1]$

$s_{ij}=\frac{\sum_{n=1}^{N_{b}}B_{i}(n)B_{j}(n)}{\sum_{n=1}^{N_{b}}B_{i}(n)\sum_{n=1}^{N_{b}}B_{j}(n)},$
(3)

where the term in the denominator is a normalisation factor and $s_{ij}=1$
means full coincidence between the two spiking series. The ensemble average of
$s_{ij}$, $S=\langle s_{ij}\rangle=\frac{2}{N(N-1)}\sum_{i,j=1,i\neq
j}^{N}s_{ij}$ is a measure of the global synchronisation in the network. In
Fig. 5(a) we plot an example of the averaged value of $S$ as a function of
$\sigma$ for an experimental CNN clustered network with $N=246$ nodes,
obtained at DIV7 as explained in Sec. 2. A transition from an asynchronous to
an almost synchronous firing is observed as the synaptic conductance $\sigma$
is increased, which confirms that the structure is suitable for inducing
synchronisation in the system.

###  4.2 Statistical complexity

Once we know the effect that the relationship between network topology and
dynamics has on the global state, we explore the effect that the presence of
the ensemble has on the single node dynamics by measuring the statistical
complexity of the single nodes along the synchronisation process. As the
typical neuronal dynamics exhibited by Eq. (1) consists of a sequence of $L$
spikes whose amplitude variability is negligible, we focused on the complexity
$C_{i}$ of the sequence of inter-spike intervals $(t_{l}-t_{l-1})$ (ISI) of
each neuron.

The ordinal patterns formalism [39] associates a symbolic sequence to a series
[40], transforming the actual values of the data series into a set of natural
numbers. To do that, the ISI series of each neuron is divided in sequences of
length $D$. In each sequence, the data values are ordered in terms of their
relative magnitudes [13], which provides the corresponding symbolic sequence.
The information content of these sequences is then evaluated as a function of
the complexity measure. The complete process is illustrated in Fig. 4.

Figure 4: Illustration of the Ordinal Patterns formalism. (a) Detection of the
maxima of the simulated (Morris-Lecar neuron) signal. (b) Extraction of the
timestamps at which these maxima have been attained. Consecutive timestamps
are subtracted (c) and their differences (ISIs) are stored in an array (d).
The time series of ISIs is divided in sequences of length $D$ ($D=3$ in this
example), and consecutive points in each sequence are compared and classified
based on their relative values (e). Each sequence of length $D$ is given a
natural number and the probability of the realization of each possible
$D$-symbols sequence is computed (f). The probabilities are then used to
construct the complexity measure.

This is a broad-field, well-established and known method, statistically
reliable and robust to noise, extremely fast in computation and with a clear
definition and interpretation in physical terms. It is derived from two also
well-established measures (divergence and entropy), also easily interpretable
when analysing non linear dynamical systems. In addition, it only requires
soft criteria, namely that the time series must be weakly-stationary, i.e.,
for $k\leq D$, the probability for $\rm{ISI}_{t}<\rm{ISI}_{t+k}$ should not
depend on time [39], and that $M>>D!$ (where $M$ is the number of points of
the entire time series of ISIs), which are easily checkable. We proceed in the
following way, as shown in Fig. 4:

  * 1.

From each single node time series in the simulated (Morris-Lecar neuron)
signal, we detect the spikes (a) and extract the duration between two
consecutive spikes (b).

  * 2.

We compute series of the inter-spike time intervals ISI (c) and save them in
an array (d), which will be our object of study.

  * 3.

The ISI series is divided in sequences of length $D$ ($D=3$ in this
illustration, (e)). We compare consecutive points in each sequence and
associate a natural number to each of them (f), ranking them based on their
relative size.

  * 4.

We count how many times a certain symbolic sequence (or pattern) $\pi$ of
length $D$ appears ($N_{\pi}$).

  * 5.

Then, we define a probability of occurrence for each pattern:
$P_{\pi}=\frac{N_{\pi}}{N_{T}}$, where $N_{T}$ is the total number of
sequences of length $D$ in which the time series is divided, i.e.
$N_{T}=(L-1)/D$, being $L$ the total number of spikes.

  * 6.

We construct a  probability distribution, which we call $P$ from now on, from
all possible symbolic sequences of length D with probability $P_{\pi}$.

Once the probability distribution $P$ is obtained, the statistical complexity
is defined. It is a measure that should be minimal both for pure noise and
absolute regularity, and provide a bounded value for other regimes. Being this
so, we need to characterise the disorder and a correcting term (i.e., a way of
comparing known probability distributions with the actual one). The
statistical complexity ($C$), as defined in Ref. [41], is the product of the
Permutation Entropy ($H$) and the Disequilibrium ($Q$).

To define the permutation entropy $H$, the first step is the evaluation of the
Shannon entropy, that gives an idea of the predictability of the series:

$S[P]=-\sum_{j=1}^{D!}p_{j}\cdot\log(p_{j})$ (4)

The permutation entropy corresponds to the normalisation of $S$ with respect
to the entropy of the uniform probability distribution, $S_{max}$:

$\displaystyle H=\frac{S}{S_{max}},\quad S_{max}=S[P_{e}],$ (5) $\displaystyle
P_{e}\equiv\\{p_{i}=1/D!\\}_{i=1,...,D!}\implies 0\leq H\leq 1$

Regarding the disequilibrium $Q$, it is a way of measuring the distance of the
actual probability distribution $P$ with the equilibrium probability
distribution $P_{e}$. This notion of distance can be acquired by several
means; in this work, we adopt the statistical distance given by the Kullback-
Leibler [42] relative entropy ($K$):

$\displaystyle K[P|P_{e}]$
$\displaystyle=-\sum_{j=1}^{D!}p_{j}\cdot\log(p_{e})+\sum_{j=1}^{D!}p_{j}\cdot\log(p_{j})=$
$\displaystyle=S[P|P_{e}]-S[P]$ (6)

where $S[P|P_{e}]$ is the Shannon cross entropy. If we now symmetrise Eq.
(4.2), we get the Jensen-Shannon divergence ($J$):

$J[P|P_{e}]=(K[P|P_{e}]+K[P_{e}|P])/2\underbracket{\rightarrow}_{(*)}J[P|P_{e}]=S[(P+P_{e})/2]-S[P]/2-S[P_{e}]/2$
(7)

where $(*)$ is simply the rewritten version in terms of $S$. Finally, we can
write the disequilibrium $Q$ as the normalised version of $J$ as:

$Q=Q_{0}J[P|P_{e}]$ (8)

with $Q_{0}={\frac{N+1}{N}\log(N+1)-2\log(2N)+\log(N)}^{-1}$, implying again
$0\leq Q\leq 1$. We then just have to multiply $H$ and $Q$ to obtain the
Complexity measure:

$C=H\cdot Q$ (9)

##  5 Results

We summarise our results for the statistical complexity in a cultured neuronal
network in Fig. 5. As commented above, in panel (a) we show as a reference the
synchronisation level vs the synaptic conductance $\sigma$ for the dynamics
simulated on top of a DIV7 experimental network. In panel Fig. 5(b) we plot
the value of $\langle C\rangle_{k}$ as a function of the conductance $\sigma$
for two nodes with high ($k=30$) and poor ($k=3$) connectivity, being $\langle
C\rangle_{k}=\sum_{[i|k_{i}=k]}C_{i}/N_{k}$, with $N_{k}$ the number of nodes
with degree $k$. The results evidences that, on that same route to
synchronisation, there exists differences between how hubs and peripheral
nodes behave due to the presence of the ensemble, even when the global
synchronisation level is still very low.

The main detail that catches our attention is that the peripheral nodes show a
greater complexity than the hubs ($\sigma=950$). To further explore this
finding, in the third panel we depict the statistical complexity vs the
degree, for this value of $\sigma$. We can extract an interesting result here:
for $\sigma=950$, there exists an anti-correlation between $\langle
C\rangle_{k}$ and $k$.

Figure 5: Dependence of the statistical complexity at the node level and its
topological role in an experimental neuronal network 7 days in vitro old with
$N=246$ Morris-Lecar neurons. (a) Synchronisation curve for the explored
values of the synaptic conductance. (b) Complexity values $\langle
C\rangle_{k}$ vs. $d$ for low ($k=3$) and high ($k=20$) degree node values.
(b) $\langle C\rangle_{k}$ vs. $k$ for the conductance value ($\sigma=950$)
marked in (b) with a vertical dashed line. Each point is the average of $6$
network realisations.

This anti-correlation observed in cultured neuronal cultures is not as evident
as the one reported in Ref. [27] in SF networks, but taking into account that
the structure, which has grown in a limited spatial domain, does not belong to
the class of power-law networks (as it was already discussed) and that the
model (Morris-Lecar) was not originally designed for this kind of neuronal
aggregates, one can conclude that the anti-correlation between the statistical
complexity $C$ and the degree $k$ is a quite robust feature.

##  6 Conclusions

In Ref. [27] we investigated the relationship between the statistical
complexity and topology in synthetically generated networks. Here, we focused
on the study of real-world topologies, as the ones exhibited by self-organised
neuronal cultures. The longitudinal study of the morphology of these networks
shows an evolution in the topology from isolated neurons to a percolated
heterogeneous topology with small-work properties.

In order to study the structure-dynamics interaction in these networks, we
simulated a dynamical model (the Morris-Lecar neuron) on top of experimental
neuronal networks at the mature developmental stage. We evidenced that, in the
weakly coupled regime, it is possible to anti-correlate the individual node
statistical complexity of the series of the neuronal inter-spike intervals
with the degree of a node. Therefore, it would be possible to infer the degree
distribution of the network from node dynamical measurements, which confirms
the result obtained for synthetic networks [27]. This approach based on the
computation of complexity values retrieved from single node dynamics, provides
a different perspective than the usual methods of network inference, since it
does not imply node-to-node calculations. Additionally, our method does not
impose the need of measuring the dynamics of every node: it can be an
incomplete measure, an still it will provide their relative roles. We hope
this approach will be useful in applications where the knowledge of the degree
distribution, instead of the detailed connectome, provides a sufficient
insight over an unknown topology and about the functioning of the underlying
system.

##  7 Acknowledgements

Financial support from the Ministerio de Economía y Competitividad of Spain
under project FIS2017-84151-P and from the Group of Research Excelence URJC-
Banco de Santander is acknowledged. A.T. and L.B.-E. acknowledge support from
the European Youth Employment Initiative.

## References

  * [1] L. M. Pecora, Synchronization conditions and desynchronizing patterns in coupled limit-cycle and chaotic systems, Phys. Rev. E 58 (1) (1998) 347–360.
  * [2] M. Barahona, L. M. Pecora, Synchronization in small-world systems, Phys. Rev. Lett. 89 (2002) 054101.
  * [3] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, D.-U. Hwang, Complex networks: Structure and dynamics, Phys. Rep. 424 (4–5) (2006) 175–308.
  * [4] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, C. Zhou, Synchronization in complex networks, Phys. Rep. 469 (3) (2008) 93 – 153.
  * [5] E. Bullmore, O. Sporns, Complex brain networks: graph theoretical analysis of structural and functional systems, Nat. Rev. Neurosci. 10 (3) (2009) 186–198.
  * [6] M. Rohden, A. Sorge, M. Timme, D. Witthaut, Self-organized synchronization in decentralized power grids, Phys. Rev. Lett. 109 (2012) 064101.
  * [7] A. Pluchino, V. Latora, A. Rapisarda, Changing opinions in a changing world: a new perspective in sociophysics, Int. J. Mod. Phys. C 16 (04) (2005) 515–531.
  * [8] N. Fujiwara, J. Kurths, A. Díaz-Guilera, Synchronization in networks of mobile oscillators, Phys. Rev. E 83 (2011) 025101.
  * [9] E. Rodriguez, N. George, J. P. Lachaux, J. Martinerie, B. Renault, F. J. Varela, Perception’s shadow: Long-distance synchronization of human brain activity, Nature 397 (6718) (1999) 430–433.
  * [10] L. Jean-Philippe, R. Eugenio, M. Jacques, V. F. J., Measuring phase synchrony in brain signals, Human Brain Mapping 8 (4) (1999) 194–208.
  * [11] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, R. Roy, Cluster synchronization and isolated desynchronization in complex networks with symmetries, Nature Communications 5 (2014) 4079.
  * [12] G. Tononi, O. Sporns, G. M. Edelman, A measure for brain complexity: relating functional segregation and integration in the nervous system., Proc. Natl. Acad. Sci. 91 (11) (1994) 5033–5037.  doi:10.1073/pnas.91.11.5033.
  * [13] A. A. Rad, I. Sendiña Nadal, D. Papo, M. Zanin, J. M. Buldú, F. del Pozo, S. Boccaletti, Topological measure locating the effective crossover between segregation and integration in a modular network, Phys. Rev. Lett. 108 (2012) 228701\.
  * [14] O. Sporns, Network attributes for segregation and integration in the human brain, Curr. Opinion Neurobiol. 23 (2013) 162–171.  doi:10.1016/j.conb.2012.11.015.
  * [15] J. Gómez-Gardeñes, Y. Moreno, A. Arenas, Paths to synchronization on complex networks, Phys. Rev. Lett. 98 (3) (2007) 034101.  doi:10.1103/PhysRevLett.98.034101.
  * [16] D. Li, I. Leyva, J. A. Almendral, I. Sendiña Nadal, J. M. Buldú, S. Havlin, S. Boccaletti, Synchronization interfaces and overlapping communities in complex networks, Phys. Rev. Lett. 101 (2008) 168701.
  * [17] A. Navas, J. A. Villacorta-Atienza, I. Leyva, J. A. Almendral, I. Sendiña Nadal, S. Boccaletti, Effective centrality and explosive synchronization in complex networks, Phys. Rev. E 92 (2015) 062820.
  * [18] P. S. Skardal, D. Taylor, J. Sun, Optimal synchronization of complex networks, Phys. Rev. Lett. 113 (14) (2014) 144101.
  * [19] M. P. Van Den Heuvel, O. Sporns, Network hubs in the human brain, Trends Cogn. Sci. 17 (2013) 683–696.
  * [20] D. Papo, M. Zanin, J. A. Pineda-Pardo, S. Boccaletti, J. M. Buldu, J. M. Buldú, Functional brain networks: great expectations, hard times and the big leap forward, Phil. Trans. R. Soc. B 369 (1653) (2014) 20130525–20130525.
  * [21] G. Zamora-López, Y. Chen, G. Deco, M. L. Kringelbach, C. Zhou, Functional complexity emerging from anatomical constraints in the brain: the significance of network modularity and rich-clubs, Sci. Rep. 6 (1) (2016) 38424\.
  * [22] G. Deco, T. J. Van Hartevelt, H. M. Fernandes, A. Stevner, M. L. Kringelbach, The most relevant human brain regions for functional connectivity: Evidence for a dynamical workspace of binding nodes from whole-brain computational modelling, NeuroImage 146 (2017) 197–210.
  * [23] T. Pereira, Hub synchronization in scale-free networks, Phys. Rev. E 82 (2010) 036201\.
  * [24] C. Zhou, J. Kurths, Hierarchical synchronization in complex networks with heterogeneous degrees, Chaos 16 (1) (2006) 015104.
  * [25] A. Arenas, A. Díaz-Guilera, C. J. Pérez-Vicente, Synchronization reveals topological scales in complex networks, Phys. Rev. Lett. 96 (2006) 114102.
  * [26] C. J. Honey, O. Sporns, L. Cammoun, X. Gigandet, J. P. Thiran, R. Meuli, P. Hagmann, Predicting human resting-state functional connectivity from structural connectivity, Proceedings of the National Academy of Sciences 106 (6) (2009) 2035–2040.
  * [27] A. Tlaie, I. Leyva, R. Sevilla-Escoboza, V. Vera-Avila, I. Sendina-Nadal, Dynamical complexity as a proxy for the network degree distribution, arXiv preprint arXiv:1807.09629 (2018).
  * [28] J.-P. Eckmann, L. Feinerman, O.and Gruendlinger, E. Moses, J. Soriano, T. Tlusty, The physics of living neural networks, Physics Reports 449 (2007) 54–76.
  * [29] E. Fuchs, A. Ayali, A. Robinson, E. Hulata, E. Ben-Jacob, Coemergence of regularity and complexity during neural network development, Dev. Neurobiol. 67 (2007) 1802–1814.
  * [30] J. van Pelt, I. Vadja, P. S. Wolters, M. Corner, G. J. A. Ramakers, Dynamics and plasticity in developing neuronal networks in vitro, Prog. Brain Res. 147 (2005) 173–188.
  * [31] M. Burrows, The Neurobiology of an Insect Brain, Oxford University Press, 1996.
  * [32] S. Anava, I. Saad, A. Ayali, The role of gap junction proteins in the development of neural network functional topology, Molecular Insect Biology 22 (2013) 457–472.
  * [33] A. Ayali, Y. Zilberstein, C. N., The role of gap junction proteins in the development of neural network functional topology, J. Exp. Biol. 205 (2002) 2825–2832.
  * [34] D. de Santos-Sierra, I. Sendiña Nadal, I. Leyva, J. Almendral, S. Anava, A. Ayali, D. Papo, S. Boccaletti, Emergence of small-world anatomical networks in self-organizing clustered neuronal cultures, PLoS ONE 9 (2014) E85828.
  * [35] D. de Santos-Sierra, I. Sendiña Nadal, I. Leyva, J. Almendral, S. Ayali, A.Anava, C. Sánchez-Ávila, S. Boccaletti, Graph-based unsupervised segmentation algorithm for cultured neuronal networks’ structure characterization and modeling, Citometry A 87 (2015) 513–523.
  * [36] C. Morris, H. Lecar, Voltage oscillations in the barnacle giant muscle fiber, Biophys. J. 35 (1) (1981) 193–213.
  * [37] B. Sancristóbal, R. Vicente, J. M. Sancho, J. García-Ojalvo, Emergent bimodal firing patterns implement different encoding strategies during gamma-band oscillations., Front. Comput. Neurosci. 7 (2013) 18.
  * [38] I. Leyva, A. Navas, I. Sendiña-Nadal, J. M. Buldu, J. A. Almendral, S. Boccaletti, Synchronization waves in geometric networks, Physical Review E 84 (6) (2011) 065101.
  * [39] C. Bandt, B. Pompe, Permutation entropy: A natural complexity measure for time series, Phys. Rev. Lett. 88 (2002) 174102.
  * [40] C. Letellier, Symbolic sequence analysis using approximated partition, Chaos, Solitons & Fractals 36 (1) (2008) 32 – 41.  doi:https://doi.org/10.1016/j.chaos.2006.06.025.
  * [41] M. Martin, A. Plastino, O. Rosso, Statistical complexity and disequilibrium, Physics Letters A 311 (2-3) (2003) 126–132.
  * [42] S. Kullback, R. A. Leibler, On information and sufficiency, The annals of mathematical statistics 22 (1) (1951) 79–86.