Electrical activity of fungi: Spikes detection and complexity analysis

Oyster fungi \emph{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.


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].
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 spikes 1 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. 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×10 −4 to 57×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.

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. 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 µ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 516 electrode pairs (channels) simultaneously.

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 Laplaces 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.

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 ≤ k ≤ 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 upperstate 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. Determining the minimum, maximum and range of amplitudes.

2.
Sorting amplitude values into the histogram bins and determining the bin width by dividing the amplitude range to the number of bins.
3. Identifying the lowest-and highest-indexed histogram bins, hb low , hb high , with non-zero counts.
4. Dividing the histogram into two sub-histograms, where the indices of the lower and upper histogram bins are hb low ≤ hb ≤ 1 2 (hb high − hb low ) and hb low + 1 2 (hb high − hb low ) ≤ hb ≤ hb high , respectively. 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.

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 ψ(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, ψ β,γ (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 β ≥ 0 and γ > 0 using Eq. 1 where ω is the angular frequency and a β,γ ≡ 2 eγ β 1 γ is the amplitude coefficient used as a real-valued normalised constant. Here, e is Eulers number, β characterises the low-frequency behaviour, and γ defines the high-frequency decay. We can rewrite Eq. 1 in the Fourier domain, parameterised by β and γ as Eq. 2.
where F (ω) is the Fourier transform of f (t), and * denotes the complex conjugate. When Ψ * β,γ (ω) 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, 1 √ s is usually used. However, we used 1 s instead, since we describe timelocalised signals by the amplitude. To recover the time-domain representation, we can use the inverse The representation of Morse wavelets can be more oscillatory when both β and γ increase, and more localised with impulses when these parameters decrease. On the other hand, increasing β and keeping γ fixed broaden the central portion of the wavelet and increase the long-time decay rate. Whereas, increasing γ by keeping β constant expands the wavelet envelope without affecting the longtime decay rate. Following explanations given in [30], we set the symmetry parameter γ to 3 and the time-bandwidth product P 2 = βγ 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. 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.
where | • | and (•) return the absolute value and the matrix transpose, respectively. Here, η is a scaling factor that we empirically set it to 240. We used g β,γ (τ, s) in Algorithm 1 to extract candidate ROIs, which are shown in Fig. 5.
Output: B -set of candidate regions.
Add min(U n + slack, τ ) to U; 10 n = n + 1; As shown in Fig. 5(c,d), some of the detected regions are either too short 2 or lack repolarisation and depolarisation periods that should be removed from B. We proposed Algorithm 2 to remove these regions, which we called them pseudo-spike and inflection regions, respectively. Figure 6 shows the results. 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.

Analytical signal envelope for locating spike pattern
To obtain the signal envelope, ξ, 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 = ∂ 2 f 4∂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 (ω) from a real-valued N -point discrete time signal L[n] are as follows: 3. Compute the N -point inverse discrete-time Fourier transform to obtain the complex discretetime analytic signal of same sample rate as the original L[n] 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, ξ[n], containing the upper, ξ H [n], and lower, ξ L [n], envelopes of L[n] (Eq. 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 seconds 3 . We proposed Algorithm 3 to locate candidate regions using signal envelope.  14 Remove the k th entry from R where R 3 (k) < ρ -see Fig. 7 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 R before applying Step 13. At this stage, although 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 (C ∪D) are used in updating R. If any ROI in R is a subset of (C ∪ D), it is added to the spike event set, F s , with an updated length. If any ROI in (C ∪ D) is a subset of R, it is added to the pseudo-spike set, 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 F p . Finally, regions with a length of fewer than 60 seconds are removed from F s and F p . Figure 8 shows results of applying Algorithm 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.

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. Waveletbased 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,000second 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.  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].  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.  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.

Complexity Analysis
To quantify the complexity of the electrical signalling recorded, we used the following measurements: 1. The Shannon entropy, H, is calculated as H = − w∈W (ν(w)/η · ln(ν(w)/η)), where ν(w) is a number of times the neighbourhood configuration w is found in configuration W , and η is the total number of spike events found in all channels of an experiment.
2. Simpson's diversity, S, is calculated as S = w∈W (ν(w)/η) 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. Space filling, D, is the ratio of non-zero entries in W to the total length of string.
4. Expressiveness, E, is calculated as the Shannon entropy H divided by space-filling ratio D, the expressiveness reflects the 'economy of diversity'. 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).
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. 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.  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 news 4 , (ii) a random sequence of alphanumeric 5 , 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.

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×10 −4 to 57×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]. 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).