Electrical activity of fungi: Spikes detection and complexity analysis
EElectrical activity of fungi: Spikes detection and complexity analysis
Mohammad Mahdi Dehshibi and Andrew Adamatzky Department of Computer Science, Universitat Oberta de Catalunya, Barcelona, Spain Unconventional Computing Laboratory, University of the West England, Bristol, UKAugust 25, 2020
Abstract
Oyster fungi
Pleurotus djamor generate actin potential like spikes of electrical potential. Thetrains of spikes might manifest propagation of growing mycelium in a substrate, transportationof nutrients and metabolites and communication processes in the mycelium network. The spikingactivity of the mycelium networks is highly variable compared to neural activity and therefore cannot be analysed by standard tools from neuroscience. We propose original techniques for detectingand 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 researchon sensorial fusion and decision making of fungi.
Keywords:
Pleurotus djamor, electrical activity, spikes, complexity
Extracellular (EC) recordings of action potentials have been widely used to record and measure neuralactivity in a number of species. The broad functionality of this method has been shown for studyingneural activity in several applications, ranging from single nerve fibres in invertebrate sensory organsto 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 foundin various physical [28, 52, 48, 20], chemical [8, 56, 57] and social systems [16, 15]. When recordedwith differential electrodes, a propagating excitation wave is manifested by spike.In our recent studies [4, 5], we demonstrated that the oyster fungi
Pleurotus djamor generateaction potential like impulses of electrical potential. We observed trains of the spontaneous spikes 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 ECspike waveform, separating it from the background activity of neighbouring cells, and sorting thecharacteristics.Lack of an algorithmic framework for exhaustive characterisation of the electrical activity of asubstrate colonised by mycelium of oyster fungi Pleurotus djamor motivated us to develop this frame-work to extract spike patterns, quantify the diversity of spiking events, and measure the complexityof fungal electrical communication. We evidenced the spiking activity of the mycelium (see an exam-ple in Fig. 1), which will enable us to develop an experimental prototype of fungi-based informationprocessing devices. Calling 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. a r X i v : . [ q - b i o . N C ] A ug a)(b) (c) Figure 1: The electrical activity of the mycelium of the grey oyster fungi. (a) Example of a dynamicsof electrical potentials recorded from eight channels of the same cluster during 63 hours. (b) Threechannels are zoomed in the inserts to show the rich combination of slow (hours) drift of base electricalpotential combined with relatively fast (minutes) oscillations of the potential. (c) All ‘classical parts ofa 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 refractoryperiod 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 ofspike 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.2e 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 variesfrom 11 × − to 57 × − . This might indicate mycelium sub-networks in different parts of thesubstrate 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. Thedetails of the proposed methods for spike detection are explained in Sect. 3. Experimental results andcomplexity analysis are presented in Sect. 4. Finally, the discussion is given in Sect. 5. 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 tentwith a silver Mylar lightproof inner lining (Green Box Tents, UK). Figure 2 shows three examples ofthe experimental set-up. (a) (b)(c)
Figure 2: Three examples of the experimental set-up with (a) in lines placement of electrodes (1 cmdistance), (b) in lines placement of electrodes (2 cm distance), (c) random electrode placement.3e 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 isolationand software-selectable sample rates all contribute to a superior noise-free resolution. We recordedelectrical activity one sample per second, where the minimum and maximum logging times were 60.04and 93.45 hours, respectively. During the recording, the logger makes as many measurements aspossible (typically up to 600 per second) and saves the average value. We set the acquisition voltagerange to 156 mV with an offset accuracy of 9 µ V at 1 Hz to maintain a gain error of 0.1%. Eachelectrode pair was considered independently with the noise-free resolution of 17 bits and conversiontime of 60 ms. In our experiments, electrode pairs were arranged in one of two configurations: randomplacement or in lines. Distance between electrodes was 1-2 cm. In each cluster, we recorded 516electrode pairs (channels) simultaneously.
A spike event can be formally defined as an extracellular signal that exceeds a simple amplitudethreshold and passes through a subsequent pair of user-specified time-voltage boxes. The spike, whichincludes depolarisation, repolarisation, and refractory periods, reflects physiological and morphologicalprocesses ongoing in mycelial networks. To extract spike events, we proposed an unsupervised methodwhich consists of three major steps.In the first step, we split the whole recording period, F ( t ), into k chunks, f k ( t ), with respectto the signal’s transitions. To determine the transitions, we estimated the state levels of the signalby its histogram and identified all regions that cross the upper-state boundary of the low state andthe lower-state boundary of the high state. Then, we calculated scale-to-frequency conversions ofthe analytic signal in each chunk using Morse wavelet basis [31]. To assess the presence of spike-likeevents, we scaled the wavelet coefficients at each frequency and obtained the sum of the scales thatwere 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 interpolationover local maximums. To determine the analytical signal, we first applied the discrete approxima-tion 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 thediscrete-time Fourier transform. From the average signal envelope, we extracted regions that fall in aconsecutive local minimum and maximum. These regions created constraints that observing them ledto the identification of spike events.In the third step, we preserved ROIs from the first step where satisfied constraints obtained inthe second step. The signal envelope could guide wavelet decomposition in an unsupervised way tocluster signal into the spike, pseudo-spike, and background activity of neighbouring cells. We detailedthe proposed method in the following sub-sections. To split the fungi electrical activity, F ( t ), with a length of t second into k chunks f k ( t ) , ≤ k ≤ t − 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 statesof the signal, we followed the following steps.1. Determining the minimum, maximum and range of amplitudes.4. Sorting amplitude values into the histogram bins and determining the bin width by dividing theamplitude 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 his-togram bins are hb low ≤ hb ≤ ( hb high − hb low ) and hb low + ( 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 pulseand the next positive-going transition. Figure 3 shows slicing results for two channels. (a) (b) Figure 3: Slicing electrical potential recordings for two channels.
The electrical activity of the mycelium exhibits modulated behaviour with variation in amplitude andfrequency over time. This feature hints that the signal can be analysed with analytic wavelets, whichare naturally grouped into even or cosine-like and odd or sine-like pairs, allowing them to capturephase variability. A wavelet ψ ( t ) is a finite energy function which projects the f ( t ) onto a family oftime-scale waveforms by translation and scaling. The Morse wavelet, ψ β,γ ( t ), is an analytic waveletwhose Fourier transforms is supported only on the positive real axis [31, 29]. This wavelet is definedin the frequency domain for β ≥ γ > ψ β,γ ( t ) = 12 π (cid:90) ∞−∞ Ψ β,γ ( ω ) e iωt d ω, Ψ β,γ ( ω ) ≡ a β,γ ω β e − ω γ × ω > ω = 00 ω < . (1)5here ω is the angular frequency and a β,γ ≡ (cid:16) eγβ (cid:17) γ is the amplitude coefficient used as a real-valuednormalised 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. φ β,γ ( τ, s ) ≡ (cid:90) ∞−∞ s ψ ∗ β,γ ( t − τs ) f ( t ) dt = 12 π (cid:90) ∞−∞ e i ωτ Ψ ∗ β,γ (s ω )F( ω ) d ω. (2)where F ( ω ) is the Fourier transform of f ( t ), and ∗ denotes the complex conjugate. When Ψ ∗ β,γ ( ω ) isreal-valued, the conjugation may be omitted. The scale variable s causes stretching or compressionof the wavelet in time. In order to reflect the energy of f ( t ) and normalise the time-domain waveletsto preserve constant energy, √ s is usually used. However, we used s instead, since we describe time-localised signals by the amplitude. To recover the time-domain representation, we can use the inverseFourier transform by f ( t ) = π (cid:82) ∞−∞ e iωt F ( ω ) d ω and ψ β,γ ( t ) = (cid:82) ∞−∞ e iωt dt = 2 πδ ( ω ), where δ ( ω ) isthe Dirac delta function.The representation of Morse wavelets can be more oscillatory when both β and γ increase, andmore localised with impulses when these parameters decrease. On the other hand, increasing β andkeeping γ 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 long-time decay rate. Following explanations given in [30], we set the symmetry parameter γ to 3 and thetime-bandwidth product P = βγ to 60. We also used L Slice and Slice ) with their Morse wavelet scalograms. (a) (b) Figure 4: Annotated spikes by the expert over with the Morse wavelet scalogram for (a)
Slice and(b) Slice . We added black arrows to point to the spike identified by the expert. The scalogram isplotted as a function of time and frequency in which the maximum absolute value at each frequencyis 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 coef-ficients can help in identifying events that may contain spikes. Hence, we proposed to use Eq. 3 fornormalisation and subsequently set all zero entries to 1.6 β,γ ( τ, s ) = | φ β,γ ( τ, s ) | (cid:124) ,g β,γ ( τ, s ) = (cid:18) η × κ β,γ ( τ, s ) − min s ( κ β,γ ( τ, s ))max s ( κ β,γ ( τ, s )) (cid:19) (cid:124) . (3)where | • | and ( • ) (cid:124) return the absolute value and the matrix transpose, respectively. Here, η is ascaling factor that we empirically set it to 240. We used g β,γ ( τ, s ) in Algorithm 1 to extract candidateROIs, which are shown in Fig. 5. Algorithm 1:
Detecting candidate regions for time-localised events.
Input : g β,γ ( τ, s ) – Scaled wavelets coefficients. Output: B – set of candidate regions. begin (cid:15) = 0 . × (max( g β,γ ( τ, s )) − min( g β,γ ( τ, s ))); max g ← set of all LocalMaximum( g β,γ ( τ, s ) , (cid:15) ) ; // LocalMaximum() returns τ ∗ if ∀ τ ∈ ( τ ∗ ± (cid:15) ) , g β,γ ( τ ∗ , s ) ≥ g β,γ ( τ, s ) . min g ← set of all LocalMinimum( g β,γ ( τ, s ) , (cid:15) ) ; U ← sort ( min g (cid:83) max g ); n = card ( U ); // card ( A ) returns number of entries in A . if n ≡ then slack ← mean (difference of two consecutive entries); Add min ( U n + slack , τ ) to U ; n = n + 1; end B ← ( U i , U i +1 ) , ∀ i ∈ { , , · · · , n − } end return B (a) (b) c) (d) Figure 5: (a,b) Identified local maxima and minima over g β,γ ( τ, s ) in (a) Slice and (b) Slice . Thesecond-row of the plot is the inverse of the first row; therefore, the marked maximums are identical tothe local minima. (c,d) Candidate regions of interest which are alternately coloured purple and greento ease visual tracking.As shown in Fig. 5(c,d), some of the detected regions are either too short or lack repolarisationand depolarisation periods that should be removed from B . We proposed Algorithm 2 to remove theseregions, which we called them pseudo-spike and inflection regions, respectively. Figure 6 shows theresults. (a) (b) Figure 6: Results of applying Algorithm 2 to (a)
Slice and (b) Slice . Two spike events are missedin Slice . Two pseudo-spike and two inflection regions still remain in Slice .Applying Algorithm 2 led to the loss of two spikes in Slice (see Fig. 6(a)) and failure to removetwo pseudo-spike and two inflection regions in Slice (see Fig. 6(b)). We found that assessing theanalytic signal by its envelope can increase the accuracy of spike detection. We observed in our previous studies [4, 5] that minimum spike length was 5 mins. lgorithm 2: Excluding pseudo-spike and inflation regions form candidate ROI.
Input : B — set of ROI, i.e. , Algorithm 1 output, f — Electrical potential. Output: C — set of wavelet-based ROIs, D — set of pseudo-spike and inflection regions. begin for i = 1 to card ( B ) do lb ← B ( i, ub ← B ( i, if ( ub − lb ) > then chunk = f [ lb · · · ub ]; minima = min ( isLocalMinimum( chunk ) ); // isLocalMinimum() and isLocalMaximum() use spline interpolation inlocating local extreme [22]. maxima = max ( isLocalMaximum( chunk ) ); if f ( minima ) < min ( f ( lb ) , f ( ub )) or f ( maxima ) > max ( f ( lb ) , f ( ub )) then C ← [ lb, ub ]; else D ← [ lb, ub ]; end end end end return C , D To obtain the signal envelope, ξ , we calculated the magnitude of its analytic signal. The analyticsignal is found using the discrete Fourier transform as implemented in Hilbert transform. To intensifyeffective peaks in the signal and, specifically, inflection regions ineffective, we calculated the secondnumerical derivation of the signal as L = ∂ f ∂t .A frequency-domain approach to approximately generate a discrete-time analytic signal is proposedin [33]. In this approach, the negative frequency half of each spectral period is set to zero, resultingin 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 asfollows:1. Compute the N -point discrete-time Fourier transform using F ( ω ) = T (cid:80) N − n =0 L [ n ] e − i πωT n ,where | ω | ≤ / T Hz and L [ n ] for 0 ≤ n ≤ N − L ( nT ) = L [ n ] at periodic time intervals of T seconds toprevent aliasing.2. Form the N -point one-sided discrete-time analytic signal transform: Z [ m ] = F [0] , for m = 02 F [ m ] , for 1 ≤ m ≤ N − F [ N ] , for m = N , for N + 1 ≤ m ≤ N − . (4)9. 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 ] = 1 N T N − (cid:88) m =0 Z [ m ] e i πmnN (5)Obtaining analytic signal in this way can satisfy two properties: (1) The real part is identical to theoriginal discrete-time sequence; (2) the real and imaginary components are orthogonal. Calculatingthe magnitude of this analytic signal yields signal envelope, ξ [ n ], containing the upper, ξ H [ n ], andlower, ξ L [ n ], envelopes of L [ n ] (Eq. 6). ξ [ n ] = | z [ n ] | (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 spikesof electrical potential shorter than 60 seconds . We proposed Algorithm 3 to locate candidate regionsusing signal envelope. Algorithm 3:
Detecting candidate spike region from signal envelope.
Input : ξ [ n ] — Envelope of signal L [ t ], n p = 60 — Minimum distance between two consecutive local extreme. Output: R — set of envelope-based ROIs. begin ξ M [ n ] = ( ξ H [ n ] + ξ L [ n ]) / [ val min , ind min ] = isLocalMinimum( ξ M [ n ] , n p ) ; [ val max , ind max ] = isLocalMaximum( ξ M [ n ] , n p ) ; // isLocalMinimum() and isLocalMaximum() locate local minimum and maximum,respectively. j ← index of the first local maximum whose value is greater than the value of the first localminimum; for i = 1 to card ( ind min ) do if j ≤ card ( ind max ) then ∆ ← val max ( j ) − val min ( i ); Add ( ind min ( i ) , ind max ( j ) , ∆) to R ; j ← j + 1; end end // R has j rows and columns, as R , R , and R . ρ = mean ( R ) − std ( R ); // mean () and std () calculate the mean and standard deviation, respectively. Remove the k th entry from R where R ( k ) < ρ – see Fig. 7(b); end return R This threshold can be changed with respect to the context of experiments. a) (b) (c)(d) (e) (f) Figure 7: Results of applying Algorithm 3 to
Slice (first row) and Slice (second row). (a,d) Can-didate regions by finding local minima and maxima of the analytic signal envelope. The regionswith arrows are also highlighted in red on the bar chart. (b,e) The absolute prominence differenceof consecutive local minimum and maximum. Regions that do not satisfy R ( k ) < ρ are coloured inred. (c,f) Regions of Interest in R . Gray rectangle with dash edge shows the correct spike, includingrepolarisation, depolarisation, and refractory periods. The purple dashed rectangle shows the regionwhose 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 correctlyidentified spikes are consonant with our findings in [4, 3]. To eliminate non-spike regions, which arehighlighted 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 theirrefractory 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 ofROIs from the intersection point into two slices. Then, the slice with the minimum length is addedto F p . Finally, regions with a length of fewer than 60 seconds are removed from F s and F p . Figure 8shows results of applying Algorithm 4. 11 lgorithm 4: Extracting fungi spike and pseudo-spike events.
Input : C , D , R — Regions of interest. Output: F s , F p — Fungi spike and pseudo-spike events, respectively. begin foreach r e ∈ R do chunk e ← [ r e · · · r e ]; foreach r w ∈ ( C ∪ D ) do chunk w ← [ r w · · · r w ]; switch chunk w , chunk e do case chunk e ⊂ chunk w do chunk w ( end ) = chunk e ( end ); F s ← chunk w ; end case chunk w ⊂ chunk e do F p ← chunk w ; end case intersect( chunk w , chunk e ) do // intersect() checks if two chunks have an intersection point. Split the concatenation of chunk w and chunk e from intersection point intotwo sub-Chunks ; F p ← sub-Chunks ; end end end end foreach r ∈ ( F s ∪ F p ) do Remove r if | r | < end end return F s , F p This section comprises of objective and complexity analyses. In the objective analysis, we showedthe 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 thecomplexity analysis, we selected complexity measures used in previous studies [49, 6, 10, 46] toquantify activity patterns that are spatio-temporally integrated and differentiated.
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 informationlike the construction of templates and the supervised setting of thresholds to detect and sort spikeevents [37, 47]. Nenadic and Burdick [37] developed an unsupervised method to detect and localisespikes in noisy neural recordings. This method benefits from continues wavelet transform. They12 a) (b)
Figure 8: Results of applying Algorithm 4 to (a)
Slice and (b) Slice . We alternatively usedred/purple for colouring spike events and blue/cyan for colouring pseudo-spike events.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 performedBayesian hypothesis testing. Finally, they combined decisions at different scales to estimate thearrival times of individual spikes.Shimazaki and Shinomoto [47] proposed an optimisation technique for selecting the bin width ofthe time-histogram. This optimisation minimised the mean integrated square in the kernel densityestimation. This method benefited from variable kernel width, which allowed grasping non-stationaryphenomena, and stiffness constant to avoid possible overfitting due to excessive freedom in the band-width 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. (a) (b) Figure 9: Results of applying proposed algorithms in [37, 47] to (a)
Slice and (b) Slice . Note thatthe wavelet-based method can only locate spike arrival time. The kernel bandwidth optimisation can,however, extract the spike region. 13oth 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). Theadaptive bandwidth kernel-based method could detect one spike in Fig. 9(b) and one pseudo-spike inFig. 9(a) and (b). While our proposed method wrongly introduced one spike event in Fig. 8(a) andhad 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, theproposed method could correctly locate 21 spikes, introduce four pseudo-spike events, overestimate tworefractory 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 proposedspike detection method. (a)(b) Figure 10: (a) Spike arrival time located by the expert. Here we used augmented pink arrow to pointto these spikes. (b) Spike regions extracted by the proposed method. Spike regions are alternativelycoloured in orange and violet. The green areas point to pseudo-spike regions that are mistaken forspikes. Blue rectangles with dash edge show overestimated refractory periods. We used black arrowsto point to the missed spikes.We applied the proposed method to six experiments where the statistical results are shown inFigs. 11 – 12 and summarised in Table 1. It should be noted that the placement of the electrodes14n two experiments was in lines with a distance of 1 cm, in two experiments it was in lines with adistance of 2 cm, and in two experiments it was random with a distance of approximately 2 cm. Theimplementation in MATLAB R2020a and details of experiments are available at [12]. (a) (b) (c)(d) (e) (f)
Figure 11: Distribution of spike event lengths with superimposed Gaussian and Adaptive bandwidthkernels [47]. (a,b) In-line electrode arrangements with a distance of 1 cm. (c,d) In-line electrodearrangements with a distance of 2 cm. (e,f) Random electrode arrangements with an approximatedistance of 2 cm.Table 1: The dominant value and bandwidth for the spike’s length and amplitude in each experimentacross all recording channels. The duration and amplitude of spikes are estimated via probabilitydensity 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 absolutevalue since we have bidirectional changes in potential. -0.00117 -0.01536 -0.01462 a) (b) (c)(d) (e) (f) Figure 12: Distribution of spike maximum amplitudes with superimposed Gaussian and Adaptivebandwidth kernels [47] for (a,b) in lines electrode placement with a distance of 1 cm. (c,d) in lineselectrode placement with a distance of 2 cm. (e,f) random electrode placement with an approximatedistance of 2 cm.These findings are aligned with the previously reported results on electrical activity of
Physarumpolycephalum [2, 3] in which we reported that Physarum spike lengths are in the range of 60-120seconds. In terms of growth, Physarum is faster than fungi. Therefore, we can now hypothesise thatfungal spikes can not be less than 60-120 seconds, with more observations.
To quantify the complexity of the electrical signalling recorded, we used the following measurements:1. The Shannon entropy, H , is calculated as H = − (cid:80) w ∈ W ( ν ( w ) /η · ln ( ν ( w ) /η )), where ν ( w ) is anumber of times the neighbourhood configuration w is found in configuration W , and η is thetotal number of spike events found in all channels of an experiment.2. Simpson’s diversity, S , is calculated as S = (cid:80) w ∈ W ( ν ( w ) /η ) . It linearly correlates with Shannonentropy for H < H . The valueof 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 , theexpressiveness reflects the ‘economy of diversity’.16 a) (b)(c) (d)(e) (f) Figure 13: Barcode-like representation of spike events in different channels for (a,b) in-line electrodearrangements 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.17. Lempel–Ziv complexity (compressibility), LZ , is evaluated by the size of binary string, n , andused to assess temporal signal diversity. Here, we represented the spiking behaviour of myceliumwith a binary string where ‘1s’ indicates the presence of a spike and ‘0s’ otherwise (see Fig. 13).6. Perturbation complexity index P CI = LZ/H .To calculate Lempel–Ziv complexity, we saved each signal as a PNG image (see two examples inFig. 14), where the ‘deflation’ algorithm used in PNG lossless compression [13, 25, 42] is a variation ofthe classical LZ77 algorithm [59]. We employed this approach as the recorded signal is a non-binarystring. We take the largest PNG file size to normalise this measurement. (a) (b)
Figure 14: Two samples from input channels, which are saved in black and white PNG format withoutaxes and annotations.To assess the signal diversity across all channels and observations, we represented each experimentas a matrix with binary entries with a row for each channel and a column for each observation. Thisbinary matrix is then concatenated observation-by-observation to form one binary string. We appliedKolmogorov 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 acrosschannels 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 inputsequence randomly shuffled. Since the value of LZ for a binary sequence of fixed length is maximal ifthe sequence is entirely random, the normalised values indicate the level of signal diversity on a scalefrom 0 to 1. Results of calculating these complexity measurements for all six setups are illustrated inFig. 15, and summarised in Tab. 2.Table 2: The mean of complexity measurements for six experiments. × − × − × × − × − × × − × − × − × − × × − × − × × − × − × a) (b)(c) (d)(e) (f) Figure 15: (a) Shannon entropy, (b) Simpson’s diversity, (c) Space filling, (d) Expressiveness, (e)LempelZiv complexity, and (f) Perturbation complexity index. All measurements are scaled to therange of [0 , ,(ii) a random sequence of alphanumeric , and (iii) a periodic sequence of alphanumeric converted tobinary strings by applying Huffman coding [26] (see barcodes in Fig. 16). Results of comparing thecomplexity of fungi electrical activity with these three forms are reported in Tab. 3. (a) (b) (c) 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, aperiodic sequence of alphanumeric along with three chunks randomly selected from our experiments.
Length LempelZivcomplexity Shannonentropy Simpson’sdiversity Spacefilling Kolmogorov PCI ExpressivenessNews 36187 0.127919 4.421728 0.999941 0.465996 0.765382 0.173096 9.49Random sequence 36002 0.125465 5.770331 0.999941 0.469835 1.001850 0.173621 12.28Periodic sequence 36006 0.127090 3.882058 0.999937 0.442426 0.076508 0.019708 8.77Chunk 1 36000 0.067611 16.194914 0.947368 0.000556 0.006307 0.000389 29150.84Chunk 2 36000 0.007250 15.478087 0.944444 0.000528 0.006727 0.000435 29326.90Chunk 3 36000 0.068417 31.680374 0.976190 0.001194 0.012613 0.000398 26523.10
We developed algorithmic framework for exhaustive characterisation of electrical activity of a sub-strate colonised by mycelium of oyster fungi
Pleurotus djamor . We evidenced spiking activity of themycelium. We found that average dominant duration of an action-potential like spike is 402 sec. Thespikes amplitudes’ depends on the location of the source of electrical activity related to the position ofelectrodes, thus the amplitudes provide less useful information. The amplitudes vary from 0.5 mV to6 mV. This is indeed low compared to 50-60 mV of intracellular recording, nevertheless understand-able due to the fact the electrodes are inserted not even in mycelium strands but in the substratecolonised 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 × − to 57 × − . This might indicated mycelium sub-networks in different parts of the substrate have been Z c o m p l e x it y Figure 17: Lempel-Ziv complexity of European languages (data from [45]) with average complexity offungal (‘fu’) electrical activity language added.transmitting different information to other parts of the mycelium network. This is somehow echoesexperimental results on communication between ants analysed with Kolmogorov complexity: longerpaths communicated ants corresponds to higher values of complexity [43].LZ complexity of fungal language (Tab. 2) is much higher than of news, random or periodicsequences (Tab. 3). The same can be observed for Shannon entropy. Kolmogorov complexity of thefungal language is much lower than that of news sampler or random or periodic sequences. Complexityof European languages based on their compressibility [45] is shown in Fig. 17, French having lowestLZ complexity 0.66 and Finnish highest LZ complexity 0.79. Fungal language of electrical activity hasminimum LZ complexity 0.61 and maximum 0.91 (media 0.85, average 0.83). Thus, we can speculatethat a complexity of fungal language is higher than that of human languages (at least for Europeanlanguages).
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 innovationprogramme 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.216] Andrew Adamatzky and Mohammad Mahdi Dehshibi. Exploring tehran with excitable medium.In Andrew Adamatzky, Selim Akl, and Georgios Ch. Sirakoulis, editors,
From Parallel to Emer-gent Computing , chapter 22, pages 475–488. CRC Press, 2019.[7] David J Aidley and DJ Ashley.
The physiology of excitable cells , volume 4. Cambridge UniversityPress Cambridge, 1998.[8] Boris P Belousov. A periodic reaction and its mechanism.
Compilation of Abstracts on RadiationMedicine , 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´elanie Boly, Simone Sarasso, Karina RCasali, Silvia Casarotto, Marie-Aur´elie Bruno, Steven Laureys, Giulio Tononi, et al. A theo-retically based index of consciousness independent of sensory processing and behavior.
Sciencetranslational medicine , 5(198):198ra105–198ra105, 2013.[11] Jorge M Davidenko, Arcady V Pertsov, Remy Salomonsz, William Baxter, and Jos´e Jalife. Sta-tionary 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 ac-tivity 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. Technicalreport, RFC 1950, May, 1996.[14] Roger Eckert and Paul Brehm. Ionic mechanisms of excitation in paramecium.
Annual review ofbiophysics and bioengineering , 8(1):353–383, 1979.[15] I Farkas, Dirk Helbing, and T Vicsek. Human waves in stadiums.
Physica A: statistical mechanicsand its applications , 330(1-2):18–24, 2003.[16] Ill´es Farkas, Dirk Helbing, and Tam´as Vicsek. Social behaviour: Mexican waves in an excitablemedium.
Nature , 419(6903):131, 2002.[17] Felix Franke, Michal Natora, Clemens Boucsein, Matthias HJ Munk, and Klaus Obermayer. Anonline spike detection and spike classification algorithm capable of instantaneous resolution ofoverlapping spikes.
Journal of computational neuroscience , 29(1-2):127–148, 2010.[18] J¨org 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. Spatialrepresentation 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.2222] Charles A Hall and W Weston Meyer. Optimal error bounds for cubic spline interpolation.
Journalof Approximation Theory , 16(2):105–122, 1976.[23] Helen G Hansma. Sodium uptake and membrane excitation in paramecium.
The Journal of cellbiology , 81(2):374–381, 1979.[24] Alan L Hodgkin and Andrew F Huxley. A quantitative description of membrane current and itsapplication 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 ofthe IRE , 40(9):1098–1101, 1952.[27] F Kaspar and HG Schuster. Easily calculable measure for the complexity of spatiotemporalpatterns.
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 eventsin noisy time series.
Proceedings of the Royal Society A: Mathematical, Physical and EngineeringSciences , 473(2200):20160776, 2017.[30] Jonathan M Lilly and Sofia C Olhede. Higher-order properties of analytic wavelets.
IEEETransactions on Signal Processing , 57(1):146–160, 2008.[31] Jonathan M Lilly and Sofia C Olhede. Generalized morse wavelets as a superfamily of analyticwavelets.
IEEE Transactions on Signal Processing , 60(11):6036–6041, 2012.[32] Zuozhi Liu, Xiaotian Wang, and Quan Yuan. Robust detection of neural spikes using sparsecoding based features.
Mathematical Biosciences and Engineering , 17(4):4257, 2020.[33] Lawrence Marple. Computing the discrete-time” analytic” signal via fft.
IEEE Transactions onsignal processing , 47(9):2600–2603, 1999.[34] Elisa Masi, Marzena Ciszak, Luisa Santopolo, Arcangela Frascella, Luciana Giovannetti, Em-manuela Marchi, Carlo Viti, and Stefano Mancuso. Electrical spiking in bacterial biofilms.
Journalof 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-machineinterface application.
IEEE Transactions on Biomedical Engineering , 51(6):905–911, 2004.2339] R Quian Quiroga, Zoltan Nadasdy, and Yoram Ben-Shaul. Unsupervised spike detection andsorting with wavelets and superparamagnetic clustering.
Neural computation , 16(8):1661–1687,2004.[40] Rodrigo Quian Quiroga, Alexander Kraskov, Christof Koch, and Itzhak Fried. Explicit encodingof multimodal percepts by single neurons in the human brain.
Current Biology , 19(15):1308–1313,2009.[41] Melinda R´acz, Csaba Liber, Erik N´emeth, Rich´ard Fi´ath, J´anos Rokai, Istv´an Harmati, Istv´anUlbert, and Gergely M´arton. Spike detection and sorting with deep learning.
Journal of NeuralEngineering , 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 tostudy 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. Inter-ictal spike detection in eeg using time series classification. In , pages 644–647. IEEE, 2020.[45] Markus Sadeniemi, Kimmo Kettunen, Tiina Lindh-Knuutila, and Timo Honkela. Complexityof 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 DMuthukumaraswamy. Increased spontaneous meg signal diversity for psychoactive doses of ke-tamine, lsd and psilocybin.
Scientific reports , 7:46421, 2017.[47] Hideaki Shimazaki and Shigeru Shinomoto. Kernel bandwidth optimization in spike rate estima-tion.
Journal of computational neuroscience , 29(1-2):171–182, 2010.[48] JC Slonczewski. Excitation of spin waves by an electric current.
Journal of Magnetism andMagnetic Materials , 195(2):L261–L268, 1999.[49] Nassim Taghipour, Hamid Haj Seyyed Javadi, Mohammad Mahdi Dehshibi, and AndrewAdamatzky. 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 spikewaveform 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-distancecommunication 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 amagnetic 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 spikedetection algorithm based on multi-channel of bect eeg signals.
IEEE Transactions on Circuitsand 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.2455] 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¨ofer. 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.