2cm2cm2cm2cm

# Electrical activity of fungi: Spikes detection and complexity analysis

Mohammad Mahdi Dehshibi  Department of Computer Science, Universitat Oberta de
Catalunya, Barcelona, Spain  Andrew Adamatzky  Unconventional Computing
Laboratory, University of the West England, Bristol, UK

###### Abstract

Oyster fungi _Pleurotus djamor_ generate actin potential like spikes of
electrical potential. The trains of spikes might manifest propagation of
growing mycelium in a substrate, transportation of nutrients and metabolites
and communication processes in the mycelium network. The spiking activity of
the mycelium networks is highly variable compared to neural activity and
therefore can not be analysed by standard tools from neuroscience. We propose
original techniques for detecting and classifying the spiking activity of
fungi. Using these techniques, we analyse the information-theoretic complexity
of the fungal electrical activity. The results can pave ways for future
research on sensorial fusion and decision making of fungi.

_Keywords:_ Pleurotus djamor, electrical activity, spikes, complexity

##  1 Introduction

Extracellular (EC) recordings of action potentials have been widely used to
record and measure neural activity in a number of species. The broad
functionality of this method has been shown for studying neural activity in
several applications, ranging from single nerve fibres in invertebrate sensory
organs to cortical neurons involved in cognition, learning and memory [19, 40,
50].

Excitation is an essential property of all living creatures, from bacteria
[34], Protists [14, 23, 9], fungi [35] and plants [51, 18, 58] to vertebrates
[24, 7, 36, 11]. Waves of excitation could be also found in various physical
[28, 52, 48, 20], chemical [8, 56, 57] and social systems [16, 15]. When
recorded with differential electrodes, a propagating excitation wave is
manifested by spike.

In our recent studies [4, 5], we demonstrated that the oyster fungi _Pleurotus
djamor_ generate action potential like impulses of electrical potential. We
observed trains of the spontaneous spikes111Calling the spikes spontaneous
means that they are not invoked by an intentional external stimulation.
Otherwise, the spikes indeed reflect physiological and morphological processes
ongoing in mycelial networks. with two types of activity, i.e., high-frequency
(period 2.6 min) and low-frequency (period 14 min). Appropriate utilisation of
this information is, however, subject to the accurate extraction of the EC
spike waveform, separating it from the background activity of neighbouring
cells, and sorting the characteristics.

Lack of an algorithmic framework for exhaustive characterisation of the
electrical activity of a substrate colonised by mycelium of oyster fungi
_Pleurotus djamor_ motivated us to develop this framework to extract spike
patterns, quantify the diversity of spiking events, and measure the complexity
of fungal electrical communication. We evidenced the spiking activity of the
mycelium (see an example in Fig. 1), which will enable us to develop an
experimental prototype of fungi-based information processing devices.

Figure 1: The electrical activity of the mycelium of the grey oyster fungi.
(a) Example of a dynamics of electrical potentials recorded from eight
channels of the same cluster during 63 hours. (b) Three channels are zoomed in
the inserts to show the rich combination of slow (hours) drift of base
electrical potential combined with relatively fast (minutes) oscillations of
the potential. (c) All ‘classical’ parts of a spike, i.e., depolarisation,
repolarisation and refractory period, can be found in this exemplar spike.
This spike has a period of 220 s, from base-level potential to refractory-like
period, and refractory period of 840 s. The depolarisation and repolarisation
rates are 0.03 and 0.009 mV/s, respectively.

We evaluated the proposed framework in comparison to the existing, in
neuroscience, techniques of spike detection [37, 47], and observed
considerable improvement in extracting spike activity periods. Evaluation of
the proposed method for detecting spikes events compared to the determined
spikes’ arrival time by an expert shows true-positive and false-positive rates
of 76% and 16%, respectively. We found that the average dominant duration of
an action-potential-like spike is 402 sec. The spikes’ amplitude varies from
0.5 mV to 6 mV and depends on the location of the electrical activity source
(the position of electrodes). We observed that the Kolmogorov complexity of
fungal spiking varies from 11$\times 10^{-4}$ to 57$\times 10^{-4}$. This
might indicate mycelium sub-networks in different parts of the substrate have
been transmitting different information to other parts of the mycelium
network, i.e., more extended propagation of excitation wave corresponds to
higher values of complexity.

The rest of this paper is organised as follows. Sect. 2 presents the
experimental setup. The details of the proposed methods for spike detection
are explained in Sect. 3. Experimental results and complexity analysis are
presented in Sect. 4. Finally, the discussion is given in Sect. 5.

##  2 Experimental set-up

A wood shavings substrate was colonised by the mycelium of the grey oyster
fungi, _Pleurotus ostreatus_ (Ann Miller’s Speciality Mushrooms Ltd, UK). The
substrate was placed in a hydroponic growing tent with a silver Mylar
lightproof inner lining (Green Box Tents, UK). Figure 2 shows three examples
of the experimental set-up.

Figure 2: Three examples of the experimental set-up with (a) in lines
placement of electrodes (1 cm distance), (b) in lines placement of electrodes
(2 cm distance), (c) random electrode placement.

We inserted pairs of iridium-coated stainless steel sub-dermal needle
electrodes (Spes Medica SRL, Italy), with twisted cables into the colonised
substrate to obtain electrical activity. Using an ADC-24 (Pico Technology, UK)
high-resolution data logger with a 24-bit A/D converter, galvanic isolation
and software-selectable sample rates all contribute to a superior noise-free
resolution. We recorded electrical activity one sample per second, where the
minimum and maximum logging times were 60.04 and 93.45 hours, respectively.
During the recording, the logger makes as many measurements as possible
(typically up to 600 per second) and saves the average value. We set the
acquisition voltage range to 156 mV with an offset accuracy of 9 $\mu$V at 1
Hz to maintain a gain error of 0.1%. Each electrode pair was considered
independently with the noise-free resolution of 17 bits and conversion time of
60 ms. In our experiments, electrode pairs were arranged in one of two
configurations: random placement or in lines. Distance between electrodes was
1-2 cm. In each cluster, we recorded 5–16 electrode pairs (channels)
simultaneously.

##  3 Proposed method

A spike event can be formally defined as an extracellular signal that exceeds
a simple amplitude threshold and passes through a subsequent pair of user-
specified time-voltage boxes. The spike, which includes depolarisation,
repolarisation, and refractory periods, reflects physiological and
morphological processes ongoing in mycelial networks. To extract spike events,
we proposed an unsupervised method which consists of three major steps.

In the first step, we split the whole recording period, $F(t)$, into $k$
chunks, $f_{k}(t)$, with respect to the signal’s transitions. To determine the
transitions, we estimated the state levels of the signal by its histogram and
identified all regions that cross the upper-state boundary of the low state
and the lower-state boundary of the high state. Then, we calculated scale-to-
frequency conversions of the analytic signal in each chunk using Morse wavelet
basis [31]. To assess the presence of spike-like events, we scaled the wavelet
coefficients at each frequency and obtained the sum of the scales that were
less than the threshold defined in Algorithm 1. Finally, we selected regions
of interest (ROI) enclosed between a consecutive local minimum and maximum
whose lengths were greater than 30 sec.

In the second step, we calculated the envelopes of the analytic signal using
spline interpolation over local maximums. To determine the analytical signal,
we first applied the discrete approximation of Laplace’s differential operator
to $f_{k}(t)$ to obtain a finite sequence of equally-spaced samples. Then, we
converted this finite sequence into a same-length sequence of equally-spaced
samples of the discrete-time Fourier transform. From the average signal
envelope, we extracted regions that fall in a consecutive local minimum and
maximum. These regions created constraints that observing them led to the
identification of spike events.

In the third step, we preserved ROIs from the first step where satisfied
constraints obtained in the second step. The signal envelope could guide
wavelet decomposition in an unsupervised way to cluster signal into the spike,
pseudo-spike, and background activity of neighbouring cells. We detailed the
proposed method in the following sub-sections.

###  3.1 Slicing fungi electrical activity

To split the fungi electrical activity, $F(t)$, with a length of $t$ second
into $k$ chunks $f_{k}(t),1\leq k\leq t-1$, we used signal transitions that
compose each pulse. To determine the transitions, we estimated the state
levels of $F(t)$ by a histogram method [1]. Then, we identified all regions
that cross the upper-state boundary of the low state and the lower-state
boundary of the high state. To estimate the states of the signal, we followed
the following steps.

  1. 1.

Determining the minimum, maximum and range of amplitudes.

  2. 2.

Sorting amplitude values into the histogram bins and determining the bin width
by dividing the amplitude range to the number of bins.

  3. 3.

Identifying the lowest- and highest-indexed histogram bins, $hb_{low}$,
$hb_{high}$, with non-zero counts.

  4. 4.

Dividing the histogram into two sub-histograms, where the indices of the lower
and upper histogram bins are $hb_{low}\leq
hb\leq\frac{1}{2}(hb_{high}-hb_{low})$ and
$hb_{low}+\frac{1}{2}(hb_{high}-hb_{low})\leq hb\leq hb_{high}$, respectively.

  5. 5.

Calculating the mean of the lower and upper histogram to compute the state
levels.

Each chunk is then enclosed between the last negative-going transitions of
every positive-polarity pulse and the next positive-going transition. Figure 3
shows slicing results for two channels.

Figure 3: Slicing electrical potential recordings for two channels.

###  3.2 Detecting time-localised events by Morse-based wavelets

The electrical activity of the mycelium exhibits modulated behaviour with
variation in amplitude and frequency over time. This feature hints that the
signal can be analysed with analytic wavelets, which are naturally grouped
into even or cosine-like and odd or sine-like pairs, allowing them to capture
phase variability. A wavelet $\psi(t)$ is a finite energy function which
projects the $f(t)$ onto a family of time-scale waveforms by translation and
scaling. The Morse wavelet, $\psi_{\beta,\gamma}(t)$, is an analytic wavelet
whose Fourier transforms is supported only on the positive real axis [31, 29].
This wavelet is defined in the frequency domain for $\beta\geq 0$ and
$\gamma>0$ using Eq. 1

$\displaystyle\psi_{\beta,\gamma}(t)=\frac{1}{2\pi}\,\int_{-\infty}^{\infty}\,\Psi_{\beta,\gamma}(\omega)\,e^{i\omega
t}\,\rm{d}\omega,$ $\displaystyle\Psi_{\beta,\gamma}(\omega)\equiv
a_{\beta,\gamma}\,\omega^{\beta}\,e^{-\omega^{\gamma}}\,\times\,\begin{cases}1&\omega>0\\\
\frac{1}{2}&\omega=0\\\ 0&\omega<0\end{cases}.$ (1)

where $\omega$ is the angular frequency and $a_{\beta,\gamma}\equiv
2\left(\frac{e\gamma}{\beta}\right)^{\frac{1}{\gamma}}$ is the amplitude
coefficient used as a real-valued normalised constant. Here, $e$ is Euler’s
number, $\beta$ characterises the low-frequency behaviour, and $\gamma$
defines the high-frequency decay. We can rewrite Eq. 1 in the Fourier domain,
parameterised by $\beta$ and $\gamma$ as Eq. 2.

$\displaystyle\phi_{\beta,\gamma}(\tau,s)\equiv\int_{-\infty}^{\infty}\frac{1}{s}\,\psi^{*}_{\beta,\gamma}(\frac{t-\tau}{s})f(t)\,\rm{d}t=\frac{1}{2\pi}\,\int_{-\infty}^{\infty}\,e^{i\omega\tau}\,\Psi^{*}_{\beta,\gamma}(s\omega)F(\omega)\,\rm{d}\omega.$
(2)

where $F(\omega)$ is the Fourier transform of $f(t)$, and $*$ denotes the
complex conjugate. When $\Psi^{*}_{\beta,\gamma}(\omega)$ is real-valued, the
conjugation may be omitted. The scale variable $s$ causes stretching or
compression of the wavelet in time. In order to reflect the energy of $f(t)$
and normalise the time-domain wavelets to preserve constant energy,
$\frac{1}{\sqrt{s}}$ is usually used. However, we used $\frac{1}{s}$ instead,
since we describe time-localised signals by the amplitude. To recover the
time-domain representation, we can use the inverse Fourier transform by
$f(t)=\frac{1}{2\pi}\,\int_{-\infty}^{\infty}\,e^{i\omega
t}\,F(\omega)\,\rm{d}\omega$ and
$\psi_{\beta,\gamma}(t)=\int_{-\infty}^{\infty}\,e^{i\omega
t}\,\rm{d}t=2\pi\delta(\omega)$, where $\delta(\omega)$ is the Dirac delta
function.

The representation of Morse wavelets can be more oscillatory when both $\beta$
and $\gamma$ increase, and more localised with impulses when these parameters
decrease. On the other hand, increasing $\beta$ and keeping $\gamma$ fixed
broaden the central portion of the wavelet and increase the long-time decay
rate. Whereas, increasing $\gamma$ by keeping $\beta$ constant expands the
wavelet envelope without affecting the long-time decay rate. Following
explanations given in [30], we set the symmetry parameter $\gamma$ to 3 and
the time-bandwidth product $P^{2}=\beta\gamma$ to 60. We also used $L1$
normalisation to have an equal magnitude in the wavelets when we have equal
amplitude oscillatory components at different scales. Figure 4 shows two
randomly selected 3000-second chunks of the fungi electrical activity (namely,
$Slice_{1}$ and $Slice_{2}$) with their Morse wavelet scalograms.

Figure 4: Annotated spikes by the expert over with the Morse wavelet scalogram
for (a) $Slice_{1}$ and (b) $Slice_{2}$. We added black arrows to point to the
spike identified by the expert. The scalogram is plotted as a function of time
and frequency in which the maximum absolute value at each frequency is used to
normalise coefficient. Frequency axis is displayed on a linear scale.

We observed that using the maximum absolute value at each frequency (level) to
normalise coefficients can help in identifying events that may contain spikes.
Hence, we proposed to use Eq. 3 for normalisation and subsequently set all
zero entries to 1.

$\displaystyle\kappa_{\beta,\gamma}(\tau,s)=|\phi_{\beta,\gamma}(\tau,s)|^{\intercal},$
$\displaystyle
g_{\beta,\gamma}(\tau,s)=\left(\eta\times\frac{\kappa_{\beta,\gamma}(\tau,s)-\min_{s}(\kappa_{\beta,\gamma}(\tau,s))}{\max_{s}(\kappa_{\beta,\gamma}(\tau,s))}\right)^{\intercal}.$
(3)

where $|\bullet|$ and $(\bullet)^{\intercal}$ return the absolute value and
the matrix transpose, respectively. Here, $\eta$ is a scaling factor that we
empirically set it to 240. We used $g_{\beta,\gamma}(\tau,s)$ in Algorithm 1
to extract candidate ROIs, which are shown in Fig. 5.

Input : $g_{\beta,\gamma}(\tau,s)$ – Scaled wavelets coefficients.

Output : $\mathcal{B}$ – set of candidate regions.

1 begin

2
$\epsilon=0.05\times(\max(g_{\beta,\gamma}(\tau,s))-\min(g_{\beta,\gamma}(\tau,s)))$;

3 $max_{g}\leftarrow$ set of all
LocalMaximum(_$g_{\beta,\gamma}(\tau,s),\epsilon$_);

// LocalMaximum(__) returns $\tau^{*}$ if
$\forall\tau\in(\tau^{*}\pm\epsilon),~{}g_{\beta,\gamma}(\tau^{*},s)\geq
g_{\beta,\gamma}(\tau,s)$.

4 $min_{g}\leftarrow$ set of all
LocalMinimum(_$g_{\beta,\gamma}(\tau,s),\epsilon$_);

5 $\mathcal{U}\leftarrow\mathbf{sort}(min_{g}\bigcup max_{g})$;

6 $n=\mathbf{card}(\mathcal{U})$;

// $\mathbf{card}(A)$ returns number of entries in $A$.

7 if  _$n\equiv 1\;(\bmod\;2)$_ then

8 slack $\leftarrow\mathbf{mean}($difference of two consecutive entries$)$;

9 Add $\mathbf{min}(\mathcal{U}_{n}+\mathrm{slack},\tau)$ to $\mathcal{U}$;

10 $n=n+1$;

11

12 end if

13 $\mathcal{B}\leftarrow(\mathcal{U}_{i},\mathcal{U}_{i+1}),~{}\forall
i\in\\{1,3,\cdots,n-1\\}$

14 end

15

return _$\mathcal{B}$_

Algorithm 1 Detecting candidate regions for time-localised events.

Figure 5: (a,b) Identified local maxima and minima over
$g_{\beta,\gamma}(\tau,s)$ in (a) $Slice_{1}$ and (b) $Slice_{2}$. The second-
row of the plot is the inverse of the first row; therefore, the marked
maximums are identical to the local minima. (c,d) Candidate regions of
interest which are alternately coloured purple and green to ease visual
tracking.

As shown in Fig. 5(c,d), some of the detected regions are either too
short222We observed in our previous studies [4, 5] that minimum spike length
was 5 mins. or lack repolarisation and depolarisation periods that should be
removed from $\mathcal{B}$. We proposed Algorithm 2 to remove these regions,
which we called them _pseudo-spike_ and _inflection_ regions, respectively.
Figure 6 shows the results.

Input : $\mathcal{B}$ — set of ROI, i.e., Algorithm 1 output,

$f$ — Electrical potential.

Output : $\mathcal{C}$ — set of wavelet-based ROIs,

$\mathcal{D}$ — set of _pseudo-spike_ and _inflection_ regions.

1 begin

2 for  _$i=1$ to $\mathbf{card}(\mathcal{B})$_ do

3 $lb\leftarrow\mathcal{B}(i,1)$;

4 $ub\leftarrow\mathcal{B}(i,2)$;

5 if  _$(ub-lb) >30$_ then

6 $chunk=f[lb\cdots ub]$;

7
$minima=\mathbf{min}(\textnormal{{isLocalMinimum(}}\textnormal{\emph{chunk}}\textnormal{{)}})$;

// isLocalMinimum(__) and isLocalMaximum(__) use spline interpolation in
locating local extreme[22].

8
$maxima=\mathbf{max}(\textnormal{{isLocalMaximum(}}\textnormal{\emph{chunk}}\textnormal{{)}})$;

9 if  _$f(minima)
<\mathbf{min}(f(lb),f(ub))~{}\mathbf{or}~{}f(maxima)>\mathbf{max}(f(lb),f(ub))$_
then

10 $\mathcal{C}\leftarrow[lb,~{}ub]$;

11

12 else

13 $\mathcal{D}\leftarrow[lb,~{}ub]$;

14

15 end if

16

17 end if

18

19 end for

20

21 end

22

return _$\mathcal{C}$ , $\mathcal{D}$_

Algorithm 2 Excluding pseudo-spike and inflation regions form candidate ROI.

Figure 6: Results of applying Algorithm 2 to (a) $Slice_{1}$ and (b)
$Slice_{2}$. Two spike events are missed in $Slice_{1}$. Two pseudo-spike and
two inflection regions still remain in $Slice_{2}$.

Applying Algorithm 2 led to the loss of two spikes in $Slice_{1}$ (see Fig.
6(a)) and failure to remove two _pseudo-spike_ and two _inflection_ regions in
$Slice_{2}$ (see Fig. 6(b)). We found that assessing the analytic signal by
its envelope can increase the accuracy of spike detection.

###  3.3 Analytical signal envelope for locating spike pattern

To obtain the signal envelope, $\xi$, we calculated the magnitude of its
analytic signal. The analytic signal is found using the discrete Fourier
transform as implemented in Hilbert transform. To intensify effective peaks in
the signal and, specifically, _inflection_ regions ineffective, we calculated
the second numerical derivation of the signal as
$L=\frac{\partial^{2}f}{4\partial t^{2}}$.

A frequency-domain approach to approximately generate a discrete-time analytic
signal is proposed in [33]. In this approach, the negative frequency half of
each spectral period is set to zero, resulting in a periodic one-sided
spectrum. The specific procedures for creating a complex-valued $N$-point ($N$
is even) discrete-time analytic signal $F(\omega)$ from a real-valued
$N$-point discrete time signal $L[n]$ are as follows:

  1. 1.

Compute the $N$-point discrete-time Fourier transform using
$F(\omega)=T\sum_{n=0}^{N-1}L[n]e^{-i2\pi\omega Tn}$, where $|\omega|\leq
1/2T$ Hz and $L[n]$ for $0\leq n\leq N-1$ is obtained by sampling a band-
limited real-valued continuous-time signal $L(nT)=L[n]$ at periodic time
intervals of $T$ seconds to prevent aliasing.

  2. 2.

Form the $N$-point one-sided discrete-time analytic signal transform:

$Z[m]=\begin{cases}F[0],&\text{ for }m=0\\\ 2F[m],&\text{ for }1\leq
m\leq\frac{N}{2}-1\\\ F[\frac{N}{2}],&\text{ for }m=\frac{N}{2}\\\ 0,&\text{
for }\frac{N}{2}+1\leq m\leq N-1.\end{cases}$ (4)

  3. 3.

Compute the $N$-point inverse discrete-time Fourier transform to obtain the
complex discrete-time analytic signal of same sample rate as the original
$L[n]$

$z[n]=\frac{1}{NT}\sum_{m=0}^{N-1}Z[m]e^{\frac{i2\pi mn}{N}}$ (5)

Obtaining analytic signal in this way can satisfy two properties: (1) The real
part is identical to the original discrete-time sequence; (2) the real and
imaginary components are orthogonal. Calculating the magnitude of this
analytic signal yields signal envelope, $\xi[n]$, containing the upper,
$\xi_{H}[n]$, and lower, $\xi_{L}[n]$, envelopes of $L[n]$ (Eq. 6).

$\xi[n]=\left|z[n]\right|$ (6)

Envelopes are determined using spline interpolation over local maxima
separated by at least $n_{p}=60$ samples. We considered $n_{p}=60$ since we
did not witness in our previous studies [4, 5] fungal spikes of electrical
potential shorter than 60 seconds333This threshold can be changed with respect
to the context of experiments.. We proposed Algorithm 3 to locate candidate
regions using signal envelope.

Input : $\xi[n]$ — Envelope of signal $L[t]$,

$n_{p}=60$ — Minimum distance between two consecutive local extreme.

Output : $\mathcal{R}$ — set of envelope-based ROIs.

1 begin

2 $\xi_{M}[n]=\left(\xi_{H}[n]+\xi_{L}[n]\right)/2$;

3 $[val_{min},ind_{min}]=$ isLocalMinimum(_$\xi_{M}[n],n_{p}$_);

4 $[val_{max},ind_{max}]=$ isLocalMaximum(_$\xi_{M}[n],n_{p}$_);

// isLocalMinimum(__) and isLocalMaximum(__) locate local minimum and maximum,
respectively.

5 $j\leftarrow$ index of the first local maximum whose value is greater than
the value of the first local minimum;

6 for  _$i=1$ to $\mathbf{card}(ind_{min})$_ do

7 if  _$j\leq\mathbf{card}(ind_{max})$_ then

8 $\Delta\leftarrow val_{max}(j)-val_{min}(i)$;

9 Add $\left(ind_{min}(i),ind_{max}(j),\Delta\right)$ to $\mathcal{R}$;

10 $j\leftarrow j+1$;

11

12 end if

13

14 end for

// $\mathcal{R}$ has $j$ rows and $3$ columns, as $\mathcal{R}_{1}$,
$\mathcal{R}_{2}$, and $\mathcal{R}_{3}$.

15 $\rho=\mathbf{mean}(\mathcal{R}_{3})-\mathbf{std}(\mathcal{R}_{3})$;

// $\mathbf{mean}()$ and $\mathbf{std}()$ calculate the mean and standard
deviation, respectively.

16 Remove the $k^{th}$ entry from $\mathcal{R}$ where
$\mathcal{R}_{3}(k)<\rho$ – see Fig. 7(b);

17

18 end

19

return _$\mathcal{R}$_

Algorithm 3 Detecting candidate spike region from signal envelope.

Figure 7: Results of applying Algorithm 3 to $Slice_{1}$ (first row) and
$Slice_{2}$ (second row). (a,d) Candidate regions by finding local minima and
maxima of the analytic signal envelope. The regions with arrows are also
highlighted in red on the bar chart. (b,e) The absolute prominence difference
of consecutive local minimum and maximum. Regions that do not satisfy
$\mathcal{R}_{3}(k)<\rho$ are coloured in red. (c,f) Regions of Interest in
$\mathcal{R}$. Gray rectangle with dash edge shows the correct spike,
including repolarisation, depolarisation, and refractory periods. The purple
dashed rectangle shows the region whose refractory period attached to a
_pseudo-spike_ event.

Figure 7(a,d) shows candidate regions in $\mathcal{R}$ before applying Step
13. At this stage, although $\mathcal{R}$ includes regions that do not observe
the spike definition (pointed by arrow in plot), the correctly identified
spikes are consonant with our findings in [4, 3]. To eliminate non-spike
regions, which are highlighted in red in Fig. 7(b,e), we applied Steps 13 and
14. Nevertheless, the output of Algorithm 3 (see Fig .7(c,f)) still contains
regions that either belong to _pseudo-spike_ /_inflection_ regions or their
refractory periods attached to a _pseudo-spike_ region.

To resolve issues in Algorithms 2 and 3, we proposed Algorithm 4 in which
regions belong to $(\mathcal{C}\cup\mathcal{D})$ are used in updating
$\mathcal{R}$. If any ROI in $\mathcal{R}$ is a subset of
$(\mathcal{C}\cup\mathcal{D})$, it is added to the spike event set,
$\mathcal{F}_{s}$, with an updated length. If any ROI in
$(\mathcal{C}\cup\mathcal{D})$ is a subset of $\mathcal{R}$, it is added to
the _pseudo-spike_ set, $\mathcal{F}_{p}$. In a case of having intersection
without observing subset condition, we split the concatenation of ROIs from
the intersection point into two slices. Then, the slice with the minimum
length is added to $\mathcal{F}_{p}$. Finally, regions with a length of fewer
than 60 seconds are removed from $\mathcal{F}_{s}$ and $\mathcal{F}_{p}$.
Figure 8 shows results of applying Algorithm 4.

Input : $\mathcal{C},\mathcal{D},\mathcal{R}$ — Regions of interest.

Output : $\mathcal{F}_{s},\mathcal{F}_{p}$ — Fungi _spike_ and _pseudo-spike_
events, respectively.

1 begin

2 foreach  _$r_{e}\in\mathcal{R}$_ do

3 $chunk_{e}\leftarrow[r_{e}^{1}\cdots r_{e}^{2}]$;

4 foreach  _$r_{w}\in(\mathcal{C}\cup\mathcal{D})$_ do

5 $chunk_{w}\leftarrow[r_{w}^{1}\cdots r_{w}^{2}]$;

6 switch  _$chunk_{w},chunk_{e}$_ do

7 case  _$chunk_{e}\subset chunk_{w}$_ do

8 $chunk_{w}(end)=chunk_{e}(end)$;

9 $\mathcal{F}_{s}\leftarrow chunk_{w}$;

10

11 end case

12 case  _$chunk_{w}\subset chunk_{e}$_ do

13 $\mathcal{F}_{p}\leftarrow chunk_{w}$;

14

15 end case

16 case  _intersect(_$chunk_{w},chunk_{e}$_)_ do

// intersect(__) checks if two chunks have an intersection point.

17 Split the concatenation of $chunk_{w}$ and $chunk_{e}$ from intersection
point into two _sub-Chunks_ ;

18 $\mathcal{F}_{p}\leftarrow$ _sub-Chunks_ ;

19

20 end case

21

22 end switch

23

24 end foreach

25

26 end foreach

27 foreach  _$r\in(\mathcal{F}_{s}\cup\mathcal{F}_{p})$_ do

28 Remove $r$ if $|r|<60$;

29

30 end foreach

31

32 end

33

return _$\mathcal{F}_{s},\mathcal{F}_{p}$_

Algorithm 4 Extracting fungi _spike_ and _pseudo-spike_ events.

Figure 8: Results of applying Algorithm 4 to (a) $Slice_{1}$ and (b)
$Slice_{2}$. We alternatively used red/purple for colouring spike events and
blue/cyan for colouring _pseudo-spike_ events.

##  4 Experimental results

This section comprises of objective and complexity analyses. In the objective
analysis, we showed the efficiency of the spike event detection method in
comparison with the existing, in neuroscience, techniques of spike detection
[37, 47] and the expert opinion in locating spikes’ arrival time. In the
complexity analysis, we selected complexity measures used in previous studies
[49, 6, 10, 46] to quantify activity patterns that are spatio-temporally
integrated and differentiated.

###  4.1 Objective analysis

Various methods have been proposed to detect and sort spike events in EC
recordings [39, 37, 38, 54, 21, 55, 17, 41, 53, 44, 32]. However, only a few
of these methods do not require auxiliary information like the construction of
templates and the supervised setting of thresholds to detect and sort spike
events [37, 47]. Nenadic and Burdick [37] developed an unsupervised method to
detect and localise spikes in noisy neural recordings. This method benefits
from continues wavelet transform. They applied multi-scale decomposition of
the signal using ‘bior1.3,’ ‘bior1.5,’ ‘Haar,’ or ‘db2’ wavelet basis. To
assess the presence of spikes, they separated the signal and noise at each
scale and performed Bayesian hypothesis testing. Finally, they combined
decisions at different scales to estimate the arrival times of individual
spikes.

Shimazaki and Shinomoto [47] proposed an optimisation technique for selecting
the bin width of the time-histogram. This optimisation minimised the mean
integrated square in the kernel density estimation. This method benefited from
variable kernel width, which allowed grasping non-stationary phenomena, and
stiffness constant to avoid possible overfitting due to excessive freedom in
the bandwidth variability. The estimated bandwidth was then used to filter
spike event regions from the signal. Figure 9 shows the results of applying
these methods to two chunks with a length of 3000 seconds.

Figure 9: Results of applying proposed algorithms in [37, 47] to (a)
$Slice_{1}$ and (b) $Slice_{2}$. Note that the wavelet-based method can only
locate spike arrival time. The kernel bandwidth optimisation can, however,
extract the spike region.

Both methods could not correctly detect all spike events that were located by
the expert. Wavelet-based method could locate three spikes in Fig. 9(b)
without detecting any spike in Fig. 9(a). The adaptive bandwidth kernel-based
method could detect one spike in Fig. 9(b) and one pseudo-spike in Fig. 9(a)
and (b). While our proposed method wrongly introduced one spike event in Fig.
8(a) and had a total-error (wrong- or non-detection) of three in Fig. 8(b).

We also compared the proposed method with the expert opinion on a randomly
selected 36,000-second chunk, i.e., 10 hours of electrical activity
recordings. In this quantitative comparison, the proposed method could
correctly locate 21 spikes, introduce four pseudo-spike events, overestimate
two refractory periods; resulting in the true-positive and false-positive
rates of 76% and 16%, respectively. Figure 10(a) shows located spikes by the
expert and Fig. 10(b) indicates the results of the proposed spike detection
method.

Figure 10: (a) Spike arrival time located by the expert. Here we used
augmented pink arrow to point to these spikes. (b) Spike regions extracted by
the proposed method. Spike regions are alternatively coloured in orange and
violet. The green areas point to _pseudo-spike_ regions that are mistaken for
spikes. Blue rectangles with dash edge show overestimated refractory periods.
We used black arrows to point to the missed spikes.

We applied the proposed method to six experiments where the statistical
results are shown in Figs. 11 – 12 and summarised in Table 1. It should be
noted that the placement of the electrodes in two experiments was in lines
with a distance of 1 cm, in two experiments it was in lines with a distance of
2 cm, and in two experiments it was random with a distance of approximately 2
cm. The implementation in MATLAB R2020a and details of experiments are
available at [12].

Figure 11: Distribution of spike event lengths with superimposed Gaussian and
Adaptive bandwidth kernels [47]. (a,b) In-line electrode arrangements with a
distance of 1 cm. (c,d) In-line electrode arrangements with a distance of 2
cm. (e,f) Random electrode arrangements with an approximate distance of 2 cm.
Table 1: The dominant value and bandwidth for the spike’s length and amplitude
in each experiment across all recording channels. The duration and amplitude
of spikes are estimated via probability density function (PDF) and adaptive
bandwidth kernel (ABK) [47]. The bold-face blue and red entries indicate the
absolute minimum and maximum values, respectively. We considered the absolute
value since we have bidirectional changes in potential.

| #Channels | #Spikes | Length (sec) | Amplitude (V)
---|---|---|---|---
PDF | ABK | PDF | ABK
Dominant | Bandwidth | Dominant | Bandwidth | Dominant | Bandwidth | Dominant | Bandwidth
#1 | 8 | 565 | 84.00 | 75.61 | 84.00 | 60.22 | 0.00003 | 0.00048 | -0.00117 | 0.00576
#2 | 5 | 447 | 366.80 | 154.31 | 625.60 | 126.47 | 0.00642 | 0.00544 | 0.00642 | 0.00667
#3 | 4 | 124 | 84.00 | 75.61 | 84.00 | 60.22 | 0.00003 | 0.00048 | -0.00117 | 0.00576
#4 | 5 | 951 | 534.12 | 80.09 | 534.12 | 84.80 | -0.00239 | 0.00301 | -0.00239 | 0.00508
#5 | 5 | 573 | 334.25 | 74.52 | 334.25 | 80.9 | -0.01536 | 0.00218 | -0.01462 | 0.00357
#6 | 15 | 862 | 1014.72 | 99.53 | 1014.72 | 92.67 | -0.00172 | 0.00381 | -0.01277 | 0.00591

Figure 12: Distribution of spike maximum amplitudes with superimposed Gaussian
and Adaptive bandwidth kernels [47] for (a,b) in lines electrode placement
with a distance of 1 cm. (c,d) in lines electrode placement with a distance of
2 cm. (e,f) random electrode placement with an approximate distance of 2 cm.

These findings are aligned with the previously reported results on electrical
activity of _Physarum polycephalum_ [2, 3] in which we reported that Physarum
spike lengths are in the range of 60-120 seconds. In terms of growth, Physarum
is faster than fungi. Therefore, we can now hypothesise that fungal spikes can
not be less than 60-120 seconds, with more observations.

###  4.2 Complexity Analysis

To quantify the complexity of the electrical signalling recorded, we used the
following measurements:

  1. 1.

The Shannon entropy, $H$, is calculated as $H=-\sum_{w\in W}(\nu(w)/\eta\cdot
ln(\nu(w)/\eta))$, where $\nu(w)$ is a number of times the neighbourhood
configuration $w$ is found in configuration $W$, and $\eta$ is the total
number of spike events found in all channels of an experiment.

  2. 2.

Simpson’s diversity, $S$, is calculated as $S=\sum_{w\in W}(\nu(w)/\eta)^{2}$.
It linearly correlates with Shannon entropy for $H<3$ and the relationships
becomes logarithmic for higher values of $H$. The value of $S$ ranges between
0 and 1, where 1 represents infinite diversity and 0, no diversity.

  3. 3.

Space filling, $D$, is the ratio of non-zero entries in $W$ to the total
length of string.

  4. 4.

Expressiveness, $E$, is calculated as the Shannon entropy $H$ divided by
space-filling ratio $D$, the expressiveness reflects the ‘economy of
diversity’.

Figure 13: Barcode-like representation of spike events in different channels
for (a,b) in-line electrode arrangements with a distance of 1 cm, (c,d) in-
line electrode arrangements with a distance of 2 cm, and (e,f) random
electrode arrangements with an approximate distance of 2 cm.

  1. 5.

Lempel–Ziv complexity (compressibility), $LZ$, is evaluated by the size of
binary string, $n$, and used to assess temporal signal diversity. Here, we
represented the spiking behaviour of mycelium with a binary string where ‘1s’
indicates the presence of a spike and ‘0s’ otherwise (see Fig. 13).

  2. 6.

Perturbation complexity index $PCI=LZ/H$.

To calculate Lempel–Ziv complexity, we saved each signal as a PNG image (see
two examples in Fig. 14), where the ‘deflation’ algorithm used in PNG lossless
compression [13, 25, 42] is a variation of the classical LZ77 algorithm [59].
We employed this approach as the recorded signal is a non-binary string. We
take the largest PNG file size to normalise this measurement.

Figure 14: Two samples from input channels, which are saved in black and white
PNG format without axes and annotations.

To assess the signal diversity across all channels and observations, we
represented each experiment as a matrix with binary entries with a row for
each channel and a column for each observation. This binary matrix is then
concatenated observation-by-observation to form one binary string. We applied
Kolmogorov complexity algorithm [27] to calculate the across channels
Lempel–Ziv complexity, $LZc$. $LZc$ captures temporal signal diversity of
single channels as well as spatial signal diversity across channels as the
result of the observation-by-observation concatenation of the binarised data
matrix. We also normalise $LZc$ by dividing the raw value by the value
obtained for the same binary input sequence randomly shuffled. Since the value
of $LZ$ for a binary sequence of fixed length is maximal if the sequence is
entirely random, the normalised values indicate the level of signal diversity
on a scale from 0 to 1. Results of calculating these complexity measurements
for all six setups are illustrated in Fig. 15, and summarised in Tab. 2.

Figure 15: (a) Shannon entropy, (b) Simpson’s diversity, (c) Space filling,
(d) Expressiveness, (e) Lempel–Ziv complexity, and (f) Perturbation complexity
index. All measurements are scaled to the range of $[0,1]$. Table 2: The mean
of complexity measurements for six experiments.

| #Channel | #Spike |  | Lempel–Ziv
---
complexity
| Shannon
---
entropy
| Simpson’s
---
diversity
| Space
---
filling
Kolmogorov | PCI | Expressiveness
#1 | 8 | 565 | 0.79 | 45.81 | 0.76 | 30.68$\times 10^{-5}$ | 30.36$\times 10^{-4}$ | 0.365 | 20.8$\times 10^{4}$
#2 | 5 | 447 | 0.91 | 63.27 | 0.98 | 35.20$\times 10^{-5}$ | 35.78$\times 10^{-4}$ | 0.021 | 18.6$\times 10^{4}$
#3 | 4 | 124 | 0.75 | 22.57 | 0.61 | 48.10$\times 10^{-5}$ | 10.94$\times 10^{-4}$ | 0.333 | 29.71
#4 | 5 | 951 | 0.93 | 123.11 | 0.89 | 57.30$\times 10^{-5}$ | 56.05$\times 10^{-4}$ | 0.072 | 23.8$\times 10^{4}$
#5 | 5 | 573 | 0.88 | 75.75 | 0.79 | 53.02$\times 10^{-5}$ | 52.80$\times 10^{-4}$ | 0.077 | 16.4$\times 10^{4}$
#6 | 15 | 862 | 0.69 | 39.96 | 0.71 | 24.20$\times 10^{-5}$ | 25.06$\times 10^{-4}$ | 0.207 | 20.4$\times 10^{4}$

In order to clarify the communication complexity in the mycelium substrate, we
also calculated the mentioned complexity measurements for the communications
in the forms of (i) pieces of
news444https://www.sciencemag.org/news/2020/07/meet-lizard-man-reptile-loving-
biologist-tackling-some-biggest-questions-evolution, (ii) a random sequence of
alphanumeric555We used available service at https://www.random.org/ , and
(iii) a periodic sequence of alphanumeric converted to binary strings by
applying Huffman coding [26] (see barcodes in Fig. 16). Results of comparing
the complexity of fungi electrical activity with these three forms are
reported in Tab. 3.

Figure 16: Binary representation of (a) pieces of news,(b) random sequence of
alphanumeric, and (c) periodic sequence of alphanumeric after applying Huffman
coding. Table 3: The complexity measurements for pieces of news, a random
sequence of alphanumeric, a periodic sequence of alphanumeric along with three
chunks randomly selected from our experiments.

| Length |  | Lempel–Ziv
---
complexity
| Shannon
---
entropy
| Simpson’s
---
diversity
| Space
---
filling
Kolmogorov | PCI | Expressiveness
News | 36187 | 0.127919 | 4.421728 | 0.999941 | 0.465996 | 0.765382 | 0.173096 | 9.49
Random sequence | 36002 | 0.125465 | 5.770331 | 0.999941 | 0.469835 | 1.001850 | 0.173621 | 12.28
Periodic sequence | 36006 | 0.127090 | 3.882058 | 0.999937 | 0.442426 | 0.076508 | 0.019708 | 8.77
Chunk 1 | 36000 | 0.067611 | 16.194914 | 0.947368 | 0.000556 | 0.006307 | 0.000389 | 29150.84
Chunk 2 | 36000 | 0.007250 | 15.478087 | 0.944444 | 0.000528 | 0.006727 | 0.000435 | 29326.90
Chunk 3 | 36000 | 0.068417 | 31.680374 | 0.976190 | 0.001194 | 0.012613 | 0.000398 | 26523.10

##  5 Discussion

We developed algorithmic framework for exhaustive characterisation of
electrical activity of a substrate colonised by mycelium of oyster fungi
_Pleurotus djamor_. We evidenced spiking activity of the mycelium. We found
that average dominant duration of an action-potential like spike is 402 sec.
The spikes amplitudes’ depends on the location of the source of electrical
activity related to the position of electrodes, thus the amplitudes provide
less useful information. The amplitudes vary from 0.5 mV to 6 mV. This is
indeed low compared to 50-60 mV of intracellular recording, nevertheless
understandable due to the fact the electrodes are inserted not even in
mycelium strands but in the substrate colonised by mycelium. The spiking
events have been characterised with several complexity measures. Most
measures, apart of Kolmogorov complexity shown a low degree of variability
between channels (different sites of the recordings). The Kolmogorov
complexity of fungal spiking varies from 11$\times 10^{-4}$ to 57$\times
10^{-4}$. This might indicated mycelium sub-networks in different parts of the
substrate have been transmitting different information to other parts of the
mycelium network. This is somehow echoes experimental results on communication
between ants analysed with Kolmogorov complexity: longer paths communicated
ants corresponds to higher values of complexity [43].

Figure 17: Lempel-Ziv complexity of European languages (data from [45]) with
average complexity of fungal (‘fu’) electrical activity language added.

LZ complexity of fungal language (Tab. 2) is much higher than of news, random
or periodic sequences (Tab. 3). The same can be observed for Shannon entropy.
Kolmogorov complexity of the fungal language is much lower than that of news
sampler or random or periodic sequences. Complexity of European languages
based on their compressibility [45] is shown in Fig. 17, French having lowest
LZ complexity 0.66 and Finnish highest LZ complexity 0.79. Fungal language of
electrical activity has minimum LZ complexity 0.61 and maximum 0.91 (media
0.85, average 0.83). Thus, we can speculate that a complexity of fungal
language is higher than that of human languages (at least for European
languages).

## Acknowledgement

We are grateful to Dr. Anna L. Powell for assisting in collecting experimental
data.

This project has received funding from the European Union’s Horizon 2020
research and innovation programme FET OPEN “Challenging current thinking”
under grant agreement No 858132.

## References

  * [1] IEEE standard on transitions, pulses, and related waveforms.  IEEE Std 181-2003, pages 1–60, 2003.
  * [2] Andrew Adamatzky.  Tactile bristle sensors made with slime mold.  IEEE Sensors journal, 14(2):324–332, 2013.
  * [3] Andrew Adamatzky.  On spiking behaviour of oyster fungi Pleurotus djamor.  Scientific reports, 8(1):1–7, 2018.
  * [4] Andrew Adamatzky.  Towards fungal computer.  Interface focus, 8(6):20180029, 2018.
  * [5] Andrew Adamatzky.  Plant leaf computing.  Biosystems, 182:59–64, 2019.
  * [6] Andrew Adamatzky and Mohammad Mahdi Dehshibi.  Exploring tehran with excitable medium.  In Andrew Adamatzky, Selim Akl, and Georgios Ch. Sirakoulis, editors, From Parallel to Emergent Computing, chapter 22, pages 475–488. CRC Press, 2019.
  * [7] David J Aidley and DJ Ashley.  The physiology of excitable cells, volume 4.  Cambridge University Press Cambridge, 1998.
  * [8] Boris P Belousov.  A periodic reaction and its mechanism.  Compilation of Abstracts on Radiation Medicine, 147(145):1, 1959\.
  * [9] MS Bingley.  Membrane potentials in amoeba proteus.  Journal of Experimental Biology, 45(2):251–267, 1966.
  * [10] Adenauer G Casali, Olivia Gosseries, Mario Rosanova, Mélanie Boly, Simone Sarasso, Karina R Casali, Silvia Casarotto, Marie-Aurélie Bruno, Steven Laureys, Giulio Tononi, et al.  A theoretically based index of consciousness independent of sensory processing and behavior.  Science translational medicine, 5(198):198ra105–198ra105, 2013\.
  * [11] Jorge M Davidenko, Arcady V Pertsov, Remy Salomonsz, William Baxter, and José Jalife.  Stationary and drifting spiral waves of excitation in isolated cardiac muscle.  Nature, 355(6358):349, 1992.
  * [12] Mohammad Mahdi Dehshibi and Andrew Adamatzky.  Supplementary material for “Electrical activity of fungi: Spikes detection and complexity analysis”.  https://doi.org/10.5281/zenodo.3997031, 08 2020.  (Accessed on 24/08/2020).
  * [13] Peter Deutsch and J Gailly.  Zlib compressed data format specification version 3.3.  Technical report, RFC 1950, May, 1996.
  * [14] Roger Eckert and Paul Brehm.  Ionic mechanisms of excitation in paramecium.  Annual review of biophysics and bioengineering, 8(1):353–383, 1979\.
  * [15] I Farkas, Dirk Helbing, and T Vicsek.  Human waves in stadiums.  Physica A: statistical mechanics and its applications, 330(1-2):18–24, 2003.
  * [16] Illés Farkas, Dirk Helbing, and Tamás Vicsek.  Social behaviour: Mexican waves in an excitable medium.  Nature, 419(6903):131, 2002.
  * [17] Felix Franke, Michal Natora, Clemens Boucsein, Matthias HJ Munk, and Klaus Obermayer.  An online spike detection and spike classification algorithm capable of instantaneous resolution of overlapping spikes.  Journal of computational neuroscience, 29(1-2):127–148, 2010.
  * [18] Jörg Fromm and Silke Lautner.  Electrical signals and their physiological significance in plants.  Plant, cell & environment, 30(3):249–257, 2007.
  * [19] Marianne Fyhn, Sturla Molden, Menno P Witter, Edvard I Moser, and May-Britt Moser.  Spatial representation in the entorhinal cortex.  Science, 305(5688):1258–1264, 2004.
  * [20] LM Gorbunov and VI Kirsanov.  Excitation of plasma waves by an electromagnetic wave packet.  Sov. Phys. JETP, 66(290-294):40, 1987.
  * [21] J Gotman and LY Wang.  State-dependent spike detection: concepts and preliminary results.  Electroencephalography and clinical Neurophysiology, 79(1):11–19, 1991.
  * [22] Charles A Hall and W Weston Meyer.  Optimal error bounds for cubic spline interpolation.  Journal of Approximation Theory, 16(2):105–122, 1976.
  * [23] Helen G Hansma.  Sodium uptake and membrane excitation in paramecium.  The Journal of cell biology, 81(2):374–381, 1979.
  * [24] Alan L Hodgkin and Andrew F Huxley.  A quantitative description of membrane current and its application to conduction and excitation in nerve.  The Journal of physiology, 117(4):500–544, 1952.
  * [25] Paul Glor Howard.  The Design and Analysis of Efficient Lossless Data Compression Systems.  PhD thesis, Citeseer, 1993.
  * [26] David A Huffman.  A method for the construction of minimum-redundancy codes.  Proceedings of the IRE, 40(9):1098–1101, 1952.
  * [27] F Kaspar and HG Schuster.  Easily calculable measure for the complexity of spatiotemporal patterns.  Physical Review A, 36(2):842, 1987.
  * [28] Ch Kittel.  Excitation of spin waves in a ferromagnet by a uniform rf field.  Physical Review, 110(6):1295, 1958.
  * [29] Jonathan M Lilly.  Element analysis: a wavelet-based method for analysing time-localized events in noisy time series.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2200):20160776, 2017.
  * [30] Jonathan M Lilly and Sofia C Olhede.  Higher-order properties of analytic wavelets.  IEEE Transactions on Signal Processing, 57(1):146–160, 2008.
  * [31] Jonathan M Lilly and Sofia C Olhede.  Generalized morse wavelets as a superfamily of analytic wavelets.  IEEE Transactions on Signal Processing, 60(11):6036–6041, 2012\.
  * [32] Zuozhi Liu, Xiaotian Wang, and Quan Yuan.  Robust detection of neural spikes using sparse coding based features.  Mathematical Biosciences and Engineering, 17(4):4257, 2020.
  * [33] Lawrence Marple.  Computing the discrete-time” analytic” signal via fft.  IEEE Transactions on signal processing, 47(9):2600–2603, 1999.
  * [34] Elisa Masi, Marzena Ciszak, Luisa Santopolo, Arcangela Frascella, Luciana Giovannetti, Emmanuela Marchi, Carlo Viti, and Stefano Mancuso.  Electrical spiking in bacterial biofilms.  Journal of The Royal Society Interface, 12(102):20141036, 2015.
  * [35] Ann M. McGillviray and Neil A.R. Gow.  The transhyphal electrical current of _Neuruspua crassa_ is carried principally by protons.  Microbiology, 133(10):2875–2881, 1987.
  * [36] Phillip G Nelson and Melvyn Lieberman.  Excitable cells in tissue culture.  Springer Science & Business Media, 2012.
  * [37] Zoran Nenadic and Joel W Burdick.  Spike detection using the continuous wavelet transform.  IEEE transactions on Biomedical Engineering, 52(1):74–87, 2004\.
  * [38] Iyad Obeid and Patrick D Wolf.  Evaluation of spike-detection algorithms fora brain-machine interface application.  IEEE Transactions on Biomedical Engineering, 51(6):905–911, 2004\.
  * [39] R Quian Quiroga, Zoltan Nadasdy, and Yoram Ben-Shaul.  Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering.  Neural computation, 16(8):1661–1687, 2004.
  * [40] Rodrigo Quian Quiroga, Alexander Kraskov, Christof Koch, and Itzhak Fried.  Explicit encoding of multimodal percepts by single neurons in the human brain.  Current Biology, 19(15):1308–1313, 2009.
  * [41] Melinda Rácz, Csaba Liber, Erik Németh, Richárd Fiáth, János Rokai, István Harmati, István Ulbert, and Gergely Márton.  Spike detection and sorting with deep learning.  Journal of Neural Engineering, 17(1):016038, 2020.
  * [42] Greg Roelofs and Richard Koman.  PNG: the definitive guide.  O’Reilly & Associates, Inc., 1999.
  * [43] Boris Ryabko and Zhanna Reznikova.  Using shannon entropy and kolmogorov complexity to study the communicative system and cognitive capacities in ants.  Complexity, 2(2):37–42, 1996.
  * [44] Shlok Sablok, Githali Gururaj, Naushaba Shaikh, I Shiksha, and Antara Roy Choudhary.  Interictal spike detection in eeg using time series classification.  In 2020 4th International Conference on Intelligent Computing and Control Systems (ICICCS), pages 644–647. IEEE, 2020.
  * [45] Markus Sadeniemi, Kimmo Kettunen, Tiina Lindh-Knuutila, and Timo Honkela.  Complexity of european union languages: A comparative approach.  Journal of Quantitative Linguistics, 15(2):185–211, 2008.
  * [46] Michael M Schartner, Robin L Carhart-Harris, Adam B Barrett, Anil K Seth, and Suresh D Muthukumaraswamy.  Increased spontaneous meg signal diversity for psychoactive doses of ketamine, lsd and psilocybin.  Scientific reports, 7:46421, 2017.
  * [47] Hideaki Shimazaki and Shigeru Shinomoto.  Kernel bandwidth optimization in spike rate estimation.  Journal of computational neuroscience, 29(1-2):171–182, 2010.
  * [48] JC Slonczewski.  Excitation of spin waves by an electric current.  Journal of Magnetism and Magnetic Materials, 195(2):L261–L268, 1999\.
  * [49] Nassim Taghipour, Hamid Haj Seyyed Javadi, Mohammad Mahdi Dehshibi, and Andrew Adamatzky.  On complexity of persian orthography: L-systems approach.  Complex Systems, 25(2):127–156, 2016.
  * [50] Caterina Trainito, Constantin von Nicolai, Earl K Miller, and Markus Siegel.  Extracellular spike waveform dissociates four functionally distinct cell classes in primate cortex.  Current Biology, 29(18):2973–2982, 2019.
  * [51] Kazimierz Trebacz, Halina Dziubinska, and Elzbieta Krol.  Electrical signals in long-distance communication in plants.  In Communication in plants, pages 277–290. Springer, 2006.
  * [52] M Tsoi, AGM Jansen, J Bass, W-C Chiang, M Seck, V Tsoi, and P Wyder.  Excitation of a magnetic multilayer by an electric current.  Physical Review Letters, 80(19):4281, 1998.
  * [53] Zimeng Wang, Duanpo Wu, Fang Dong, Jiuwen Cao, Tiejia Jiang, and Junbiao Liu.  A novel spike detection algorithm based on multi-channel of bect eeg signals.  IEEE Transactions on Circuits and Systems II: Express Briefs, 2020\.
  * [54] Scott B Wilson and Ronald Emerson.  Spike detection: a review and comparison of algorithms.  Clinical Neurophysiology, 113(12):1873–1881, 2002.
  * [55] Scott B Wilson, Christine A Turner, Ronald G Emerson, and Mark L Scheuer.  Spike detection ii: automatic, perception-based detection and clustering.  Clinical neurophysiology, 110(3):404–411, 1999.
  * [56] AM Zhabotinsky.  Periodic processes of malonic acid oxidation in a liquid phase.  Biofizika, 9(306-311):11, 1964.
  * [57] Anatol M Zhabotinsky.  Belousov-zhabotinsky reaction.  Scholarpedia, 2(9):1435, 2007.
  * [58] Matthias R Zimmermann and Axel Mithöfer.  Electrical long-distance signaling in plants.  In Long-Distance Systemic Signaling and Communication in Plants, pages 291–308. Springer, 2013.
  * [59] Jacob Ziv and Abraham Lempel.  A universal algorithm for sequential data compression.  IEEE Transactions on information theory, 23(3):337–343, 1977.