Emergence of behaviour in a selforganized living matter network
Abstract
What is the origin of behaviour? Although typically associated with a nervous system, simple organisms also show complex behaviours. Among them, the slime mold Physarum polycephalum, a giant single cell, is ideally suited to study emergence of behaviour. Here, we show how locomotion and morphological adaptation behaviour emerge from selforganized patterns of rhythmic contractions of the actomyosin lining of the tubes making up the networkshaped organism. We quantify the spatiotemporal contraction dynamics by decomposing experimentally recorded contraction patterns into spatial contraction modes. Notably, we find a continuous spectrum of modes, as opposed to a few dominant modes. Our data suggests that the continuous spectrum of modes allows for dynamic transitions between a plethora of specific behaviours with transitions marked by highly irregular contraction states. By mapping specific behaviours to states of active contractions, we provide the basis to understand behaviour’s complexity as a function of biomechanical dynamics.
Editor's evaluation
We have judged that the response to the referee's residual comments are sufficient to allow this paper to proceed to publication. In particular, the detailed analysis of the mode spectrum and its relationship to behavior is novel and possibly of general use in this field. Also, the experimental data per se should be interesting to a wide spectrum of readers.
https://doi.org/10.7554/eLife.62863.sa0Introduction
Survival in changing environments requires from organisms the ability to transition between diverse behaviours (Angilletta and Sears, 2011; Wong and Candolin, 2014). In higher organisms, a plethora of neural dynamics enable this capacity, ranging from almost random to strongly correlated firing patterns of neurons (Mochizuki et al., 2016). Decoding the origin of behaviour from neuronal activity has been called the ‘holy grail of neuroscience’ (Bando et al., 2019), a task especially challenging given the vastly complex networks of neurons (Berman, 2018). Significant progress has been made by simultaneous tracking of neuronal activity and behaviour – defined as trajectories through spaces of postural dynamics – in the fruit fly Drosophila melanogaster (Honegger et al., 2020) and the nematode Caenorhabditis elegans (Nguyen et al., 2016). Behaviours of these systems have been identified as lowdimensional (Stephens et al., 2008) and hierarchical (Berman et al., 2016).
While these discoveries have advanced our understanding of the origin of behaviour, the complexity and size of biological neural networks make the acquisition and interpretation of experimental data especially challenging. Curiously, organisms without a nervous system may offer an ideal intermediate step towards understanding behaviour. Certain nonneural organisms readily transition between a multitude of behaviors similar in dynamic variability to that of organisms with a nervous system (Berg and Brown, 1972; Otto and Kessin, 2001; McMains et al., 2008; BenJacob et al., 1994; BenJacob et al., 2000; Wan and Goldstein, 2014; Wan, 2018) and thus provide the opportunity to study the link between the underlying biophysical process and behaviour.
A nonneural organism with an exceptionally versatile behavioural repertoire is the slime mould Physarum polycephalum  a unicellular, networkshaped organism (Sauer, 1982) of macroscopic dimensions, typically ranging from a millimeter to tens of centimeters. P. polycephalum’s complex behaviour is most impressively demonstrated by its ability to solve spatial optimisation and decisionmaking problems (Nakagaki et al., 2000; Tero et al., 2010; Nakagaki and Guy, 2007; Dussutour et al., 2010; Reid et al., 2016), exhibit habituation to temporal stimuli (Boisseau et al., 2016), and use exploration versus exploitation strategy (Aono et al., 2014). Recently, P. polycephalum was found capable of encoding memory about food source locations in the hierarchy of its body plan (Kramar and Alim, 2021) in a process much reminding of synaptic facilitation– the brain’s way of creating memories (Jackman and Regehr, 2017). The generation of such rich behaviour requires a mechanism allowing not only for longrange spatial coordination but also the flexibility to enable switching between different specific behavioural states.
The behaviour generating mechanism in P. polycephalum are the active, rhythmic, crosssectional contractions of the actomyosin cortex lining the tube walls (Yoshimoto and Kamiya, 1984; Ueda et al., 1986; Kamiya et al., 1988). The contractions drive cytoplasmic flows throughout the organism’s network (Iima and Nakagaki, 2012; Alim et al., 2013), transporting nutrients and signalling molecules (Alim et al., 2017). Cytoplasmic flow is responsible for mass transport across the organism and thereby contractions directly control locomotion behaviour (Rieu et al., 2015; Lewis et al., 2015; Zhang et al., 2017; Bäuerle et al., 2020; Rodiek et al., 2015).
So far, only one type of networkspanning peristaltic contraction pattern has been described experimentally (Alim et al., 2013; Oettmeier et al., 2017). However, for small P. polycephalum plasmodial fragments various other shortrange contraction patterns have been observed (Lewis et al., 2015; Zhang et al., 2017) and predicted by theory of active contractions (Bois et al., 2011; Radszuweit et al., 2013; Radszuweit et al., 2014; Julien and Alim, 2018; Kulawiak et al., 2019). Similarly, up to now unknown complex, largescale contraction patterns might play a role in generating the behaviour of large P. polycephalum networks. Furthermore, transitions between such largescale patterns are needed to allow for switching between specific behaviours, for example taking sharp turns during migration in the absence of stimuli (Rodiek and Hauser, 2015).
Here, we decompose experimentally recorded contractions of a large P. polycephalum network of stable morphology into a set of physically interpretable contraction modes using Principal Component Analysis. We find a continuous spectrum of modes and high variability in the activation of modes along this spectrum. By perturbing the network with an attractive stimulus, we show that the resulting locomotion response is coupled to a selective activation of regular contraction patterns. Guided by these observations, we design an experiment on a P. polycephalum specimen reduced in morphological complexity to a single tube. This allows us to quantify the causal relation between locomotion behaviour, cytoplasmic flow rate and varying types of contraction patterns, thus revealing the central role of dynamical variability in generating different behaviours.
Results
Continuous spectrum of contraction modes reveals large variability in organism’s contraction dynamics
To characterize the contraction dynamics of a P. polycephalum network, we record contractions using brightfield microscopy (Video 1) and decompose this data into a set of modes using Principal Component Analysis (PCA). At first, networks in brightfield images are skeletonized, with every single skeleton pixel representing the local tube intensity as a measure of the local contraction state (Bäuerle et al., 2017). Thus, any network state at a time t_{i} is represented by a list of pixels, ${\overrightarrow{I}}^{{t}_{i}}$, along the skeleton, see Figure 1A and ‘Data processing’ (Appendix 1). Performing PCA on this data results in a linear decomposition of the intensity vectors ${\overrightarrow{I}}^{{t}_{i}}$ into a basis of modes ${\overrightarrow{\varphi}}_{i}$:
See ‘Principal Component Analysis (PCA) (Appendix 2)’ for details. The modes, ${\overrightarrow{\varphi}}_{\mu}$, are orthonormal eigenvectors of the covariance matrix of the data and represent linearly uncorrelated contraction patterns of the network, and ${a}_{\mu}^{{t}_{i}}$ denotes the timedependent coefficients of the modes.
We rank modes according to the magnitude of their eigenvalues. Contrary to the small number of large eigenvalues found in a number of biological systems (Stephens et al., 2008; Jordan et al., 2013; Gilpin et al., 2016), here the spectrum of relative eigenvalues, see ‘Principal component analysis (PCA)’ (Appendix 2) for technical details, is continuous with no clear cutoff (Figure 1B) and as a result the contraction dynamics is highdimensional. Notably, this is even the case when we disregard eigenvalues which lie below the upper noise bound (black line), computed from randomised data. Therefore, PCA does not directly lead to a dimensionality reduction of the data. Instead, we here investigate the characteristics of mode dynamics that result from a continuous spectrum and how these shape the organism’s behaviour.
The highestranking modes shown in Figure 1C(i) have a smooth spatial structure that varies on the scale of network size. As we will discuss below, such largescale modes are associated with the long wavelength peristalsis observed in Iima and Nakagaki, 2012; Alim et al., 2013. Interestingly, we also find modes highlighting specific morphological characteristics of the network. For example, the structure of mode ${\overrightarrow{\varphi}}_{4}$, Figure 1C(i), corresponds to the thickest tubes of the network Figure 1A, which suggests a special role of these tubes in the functioning of the network. Finally, as we go to lower ranked modes, the spatial structure of the modes becomes increasingly finer. Yet, despite lacking an obvious interpretation for their structures, like for mode ${\overrightarrow{\varphi}}_{30}$, Figure 1—figure supplement 1, it is not possible to ignore their contribution relative to highranking modes.
Next, we turn to the timedependent coefficients of modes shown in Figure 1C(ii). In accordance with the known rhythmic contractions (Kamiya, 1960) the coefficient $a}_{1$ of the highest ranked mode ${\overrightarrow{\varphi}}_{1}$ oscillates with a typical period of $T\sim 100\phantom{\rule{thinmathspace}{0ex}}\mathrm{sec}$ . Most strikingly, amplitudes of mode coefficients vary significantly over time  even on orders of magnitude, as shown in Figure 1—figure supplement 2.
To map out the complexity of contractions over time, we define a set of significant modes for every time point. We quantify the activity of a mode by its relative amplitude
where $\stackrel{~}{{a}_{\mu}^{2}}$ denotes the amplitude of the square of the mode’s coefficient. By definition the sum over the relative amplitudes of all modes is normalized to one at any given time, ${\sum}_{\mu}{p}_{\mu}^{{t}_{i}}=1$. For any given time point, we order the modes by their relative amplitude from largest to smallest and take the cumulative sum of their values until a chosen cutoff percentage is reached, see Figure 2A. We find that the percentage of modes required to reach a specified cutoff value varies considerably over time. For a 90% cumulative amplitude cutoff, we find that on average 6.06% ($\approx 70$ modes) of the 1500 modes are significant. As discussed in more detail in ‘Choice of the cutoff of mode coefficient amplitudes’ (Appendix 6), defining a cutoff for the cumulative sum of mode amplitudes is related to the problem of defining a cutoff for a continuous spectrum of eigenvalues. One common method is to define the cutoff with respect to the largest eigenvalue of the spectrum computed from a randomised version of the original data (Berman et al., 2014). In ‘Choice of the cutoff of mode coefficient amplitudes’ (Appendix 6), we find that the 90% cumulative amplitude cutoff considered above is consistent with this definition of cutoff for eigenvalues. As an important feature, we observe that there is large variation in the number of significant modes over time, with a standard deviation of 36.96% from the mean value. This is an indicator for the complexity of the contractions in the network.
Apart from the number of significant modes, the dynamics of the network depend on the temporal correlation of modes. While the modes form a spatially uncorrelated basis, the temporal correlation of mode activation is nontrivial. In Figure 2B, we show the distribution of temporal correlations between mode coefficients as a function of the number of significant modes, see ‘Distribution of temporal correlations’ (Appendix 3) for technical details. For a small number of significant modes, the coefficients are strongly (anti)correlated in time, while for a large number of significant modes, correlations values between coefficients are more uniformly distributed. Here, correlated coefficients result in coordinated pumping behaviour/contractions, while least correlated coefficients coincide with irregular networkwide contractions. The above analysis shows that the dynamics of network contractions covers a wide range in complexity, from superposition of few largescale modes strongly correlated in time, to superpositions of many modes of varying spatial scale and temporal correlations. This gives rise to strong variability in the regularity of the contraction dynamics over time. Up to now, we investigated an ‘idle’ network not performing a specific task, so we next stimulate the network to provoke a specific behaviour and scrutinize how the continuous spectrum of modes contributes to it.
Stimulus response behaviour is paired with activation of regular, largescale contraction patterns interspersed by manymode states
To probe the connection between a specific behaviour and network contraction dynamics, we next apply a food stimulus to the same network, see Figure 3A. Food acts as an attractant and causes locomotion of the organism toward the stimulus in the long term. The stimulus immediately triggers the tubes in the network to grow in a concentric region around the stimulus site. Also, the thick transport tubes oriented toward stimulus location increase their volume, see Figure 3A. Altogether these morphological changes are typical for the specific behaviour induced here, namely the generation of a new locomotion front.
In Figure 3B, we quantify this stimulus response behaviour by tracking the growth of the most active regions of the network, defined by the boxes shown in the 81 min in Figure 3A. The tracked regions are located on opposing sides of the network. Starting approximately at 85 min, the part of the network next to the stimulus site grows rapidly (burgundy curve in Figure 3B), at the expense of the fanshaped locomotion front in the lower left corner of the network (green curve in Figure 3B). In Figure 3—figure supplement 2, we additionally show that prior to the stimulus, the network grows the fanlike shaped locomotion front in the lower left corner. Taken together, the application of the stimulus leads to a reversal of the network’s growth direction.
To identify potential changes in the contraction dynamics due to stimulus application, we perform PCA on a 700 frames long subset of the data subsequent to the ‘idle’ data of the previous section. First, we rediscover a continuous spectrum of modes, see Figure 3—figure supplement 1, resembling that of the ‘idle’ dynamic state. However, now the highestranked contraction modes, see Figure 3C, show spatial patterns which can be directly related to the network’s growth behaviour. This includes activation of the upper region of the network close to the stimulus, as well as activation of the thick tubes extending from top to bottom of the network. In fact, for more than 500 frames after the stimulus has been applied, the rhythmic contraction dynamics of the network are dominated by the three highestranked modes, see Figure 3D and Figure 3—figure supplement 3 for the oscillatory dynamics of mode coefficients. During this period, every time a single mode is the most active one for a duration of $\S gt;30$ frames, its amplitude exceeds that of any other mode by 20–30%.
Next we link the stimulusinduced reversal in growth direction to the changes in the contraction pattern. Specifically, we observe that the time interval of growth reversal Figure 3B coincides with the activation of the thirdranked mode ${\overrightarrow{\varphi}}_{3}$, (orange curve in Figure 3D), as indicated by the pink shaded box extending across Figure 3D and B. The structure of this mode shows a clear distinction of the growth area close to the stimulus and an activation of the two thick tubes stretching from bottom to the top of the network. This mode is followed by an activation of mode ${\overrightarrow{\varphi}}_{2}$ (blue curve), clearly marking the growth region within its spatial structure.
Finally, over time the growth of the stimulus response region tapers off and we find reactivation of mode ${\overrightarrow{\varphi}}_{1}$ (red curve) which was the dominant mode before stimulus application. We note that the spatial structure of mode ${\overrightarrow{\varphi}}_{1}$ is remarkably similar to mode ${\overrightarrow{\varphi}}_{1}$, the topranked mode that we find for PCA on the prestimulus ‘idle’ data Figure 3—figure supplement 2B. The reactivation of this mode indicates that this contraction pattern is intrinsic to the network and is not simply erased by the stimulus.
Strikingly, the regular contraction dynamics shown in Figure 3D are interspersed with manymode states where the number of significant modes increases considerably, see Figure 3E. The number of significant modes oscillates after the stimulus. The oscillation maxima coincide with times at which the organism switches from one dominant contraction pattern to another, as indicated by the blueshaded boxes extending across Figure 3D and E. Our results suggest that prolonged regular dynamics dominated by a few or even a single mode are associated with specific behaviour like locomotion and growth, while the manymode states seem to serve as transition states between them.
While the network morphology is characteristic for P. polycephalum, reducing network complexity may help to conclude on the role of regular dynamics in driving specific behaviours, and the role of manymode states and the therefrom arising continuous distribution of modes.
Number of significant modes determines maximum cytoplasmic flow rate in the minimal morphological representation of the network
We next perform exactly the same course of experiments as before but on a P. polycephalum specimen reduced in complexity to a single tube with a locomotion front at either end, see inset in Figure 4 and Video 2. Strikingly, when performing PCA on this specimen of simple morphology we again find a continuous spectrum of modes (Figure 5—figure supplement 1) and large variability, including spikes of manymode states, in the number of significant modes (Figure 5A). This observation finally underlines that the continuous spectrum of modes and its variability in activation is intrinsic to the organism’s behaviour, ruling out that the complexity of contraction modes only mirrors morphological complexity. Foremost, this minimal constituent of a network allows us now to directly map the effect of variations in the contraction dynamics onto behaviour.
From the experimentally quantified tube contractions, we calculate the maximal flow rate at any point along the tube (Li and Brasseur, 1993) and over time correlate the strength of the flow rates, driving locomotion behaviour at the tube ends, with the number of significant modes, see ‘Flow rate calculation in a P. polycephalum cell with singletube morphology’ (Appendix 4). For both the flow rate at the left and right end of the tube, shown in Figure 4, and Figure 4—figure supplement 1, respectively, we find that large flow rates are only achieved when the number of significant modes is small. We had previously found that few significant modes are highly (anti)correlated, whereas states with many significant modes are not, see Figure 2B. This observation now confirms our physical intuition that the irregularity of states consisting of many modes goes hand in hand with reduced pumping efficiency and thus unspecific behaviour. Since a small number of significant modes not necessarily always implies a large flow rate, we next turn to analyze their exact spatial structure and instantaneous temporal correlation to determine how cytoplasmic flow rates impact behaviour.
Instantaneous coupling and selective activation of modes determine locomotion behaviour
We now demonstrate the impact of changes in the dynamics of a small number of modes on the organism’s behaviour. For this, we quantify the locomotion behaviour of the single tube by tracking the area of the locomotion fronts protruding from each end of the tube over time, see Figure 5A. The growth curves of the tube ends are shown in Figure 5B. While initially the right end is protruding faster at the expense of the left end, a food stimulus applied to the left end of the tube reverses the direction of locomotion.
As for the network, we use PCA to analyse the contraction dynamics of the single tube and link it to behaviour. We apply PCA to contraction data along the tube which we parameterize by a longitudinal coordinate. The spatial shapes of the two topranked modes ${\overrightarrow{\varphi}}_{1}$ and ${\overrightarrow{\varphi}}_{2}$ approximate Fourier modes, see Figure 5C and Figure 5—figure supplement 2. Examining the activation of modes, we find that over long time intervals, and in particular after the stimulus, the two topranked modes dominate the tube’s contraction dynamics, see Figure 5D. To illustrate the connection between the nature of tube contraction dynamics and locomotion behaviour, we pick two representative time intervals after the stimulus where either only mode ${\overrightarrow{\varphi}}_{1}$, or modes ${\overrightarrow{\varphi}}_{1}$ and ${\overrightarrow{\varphi}}_{2}$ equally, dominate overall, see vertical pink bars in Figure 5D. During the first interval when mode ${\overrightarrow{\varphi}}_{1}$ alone is dominating, the tube is driven by a standing wave contraction pattern  yielding only a low cytoplasmic flow rate. Correspondingly, the size of the locomotion front at either end shows no significant change in area during this interval. In contrast, during the interval when both modes ${\overrightarrow{\varphi}}_{1}$ and ${\overrightarrow{\varphi}}_{2}$ are equally active, the resulting superposition is a lefttraveling wave producing a large cytoplasmic flow rate in that direction. The lefttraveling wave is in accordance with the growth of the left and retraction of the right locomotion front as quantified in Figure 5B. See ‘Mode superpositions in a P. polycephalum cell with singletube morphology’ (Appendix 4) for more details. In Figure 5E, we highlight the most pronounced manymode states during changes of dominant contraction dynamics.
These two examples solve the conundrum of Figure 4, which shows that a small number of significant modes does not necessarily lead to high cytoplasmic flow rates. Yet, the direct mapping of contraction dynamics onto ensuing cytoplasmic flows confirms that a small number of significant modes is associated with specific behaviour. High cytoplasmic flow rates at the tube ends drive locomotion, while lower flow rates likely lead to other behaviours such as mixing. Furthermore, manymode states seem necessary for transitions in a multibehavioural space.
Our explanation of behaviour – from contractions via flows to locomotion behaviour – in the single tube is a template for an analogous explanation in the network morphology. The analogy is justified by the strong resemblance of the continuous mode spectrum, dynamics of significant modes, activation of regular contraction patterns and the nature of growth behaviour in both the network and single tube. Therefore, while it is beyond the scope of this study, we expect a detailed analysis of the link between contractions and flows in the network morphology to yield qualitatively similar results to those of the single tube, thus completing the mechanism of behaviour generation.
Discussion
To uncover the origin of behaviour in P. polycephalum, we quantified the dynamics of this living matter network and linked it to its emerging behaviour. The simple build of this nonneural organism allows us to trace contractions of the actomyosinlined tubes, compute cytoplasmic flows from the contractions and finally link these dynamics to the emerging mass redistribution and wholeorganism locomotion behavior. Decomposing the contractions across the network into individual modes, we discover a large intrinsic variability in the number of significant modes over time along a continuous spectrum of modes. By triggering locomotion through application of a stimulus, we identify that states with few significant modes and regular contraction patterns correspond to specific behaviors, in this case locomotion. Yet, irregular contraction patterns consisting of a large number of significant modes are also present, particularly marking the transitions between different regular contraction states. The use of an organism with a singletube morphology allows us to obtain quantitative insights into the mechanism connecting contraction dynamics and locomotion behavior and in first approximation serves as an analogue system for the large P. polycephalum with network morphology. Our findings suggest that a continuous spectrum of contraction modes allows the living matter network P. polycephalum to quickly transition between a multitude of behaviours using the superposition of multiple contraction patterns.
Networks are ubiquitous in biology, including examples such as ecological networks (García Martín and Goldenfeld, 2006) and biomolecular interaction networks (Albert, 2005). Measurable quantities of these networks, for instance the degree distribution of the network, typically follow continuous distributions and are oftentimes powerlaws. The spectrum of eigenvalues Figure 1B that we find for the contraction dynamics in P. polycephalum may similarly suggest a powerlaw. However, the presence of a powerlaw is generally difficult to prove and interpret. Instead, our sole focus is on the continuous nature of the spectrum. It is important to emphasise that the continuity of the eigenvalue spectrum is not simply the result of the organism’s complex network morphology. This is demonstrated by the fact that we find a similar spectrum also for the singletube morphology Figure 5—figure supplement 1. Therefore, here the continuous spectrum of eigenvalues is distinctively a property of the dynamic state of the organism.
Our observation of interlaced regular and irregular contraction patterns in P. polycephalum reminds of the strongly correlated or random firing patterns of neurons in higher organisms (Mochizuki et al., 2016). In neural organisms, stereotyped behaviours are associated with controlled neural activity, as for example for locomotion in C. elegans (Liu et al., 2018) or the behavioural states of the fruit fly Drosophila melanogaster (Berman et al., 2014; Berman et al., 2016). Variability in the dynamics of behaviour is also widely observed in these neuronal organisms (Grobstein, 1994; Renart and Machens, 2014; Werkhoven et al., 2021; Honegger et al., 2020; Ahamed et al., 2020). It is thus likely that the transition role of irregular states consisting of many significant modes observed here for P. polycephalum parallels the mechanisms of generating behaviour in the more complex forms of life.
P. polycephalum is renowned for its ability to make informed decisions and navigate a complex environment (Nakagaki et al., 2000; Tero et al., 2010; Nakagaki and Guy, 2007; Dussutour et al., 2010; Reid et al., 2016; Boisseau et al., 2016; Aono et al., 2014; Ueda et al., 1976; Miyake et al., 1991). It would be fascinating to next follow the variability of contraction dynamics during more complex decisionmaking processes. Furthermore, it would be interesting to observe ‘idle’ networks during foraging over tens of hours. It is likely that the contraction states with many significant modes here act as noisy triggers that can spontaneously cause the organism to reorient its direction of locomotion.
In the context of P. polycephalum’s foraging behaviour, another exciting line of research opened by our results is the link between contraction modes and the organism’s metabolic changes. The foraging networks displays a plethora of morphological patterns which are linked to the underlying metabolic states (Takamatsu et al., 2017; Lee et al., 2018). It has recently been shown that in the neural organism Drosophila melanogaster, behaviour stemming from neural activity causes largescale changes in metabolic activity (Mann et al., 2021). Exploring the relationship between behaviour emergence and metabolism in P. polycephalum will bring key insight about the interplay between the mechanical and the biochemical machinery of the organism.
P. polycephalum’s bodyplan as a fluidfilled living network with emerging behaviour finds its theoretical counterpart in theories for active flow networks developed recently (Woodhouse et al., 2016; Forrow et al., 2017). Strikingly, these theories predict selective activation of thick tubes which we observe in the living network as well, prominently appearing among the top ranking modes, see ${\overrightarrow{\varphi}}_{4}$ in Figure 1C(i) or ${\overrightarrow{\varphi}}_{3}$ in Figure 3C. This is a first hint that dynamics states arising from first principles in active flow networks could map onto behavioural and transition states observed here.
Likely our most broadly relevant finding in this work is that irregular dynamics, here arising in states with many significant modes, play an important role in switching between behaviours. This should inspire theoretical investigations to embrace irregularities rather than focusing solely on regular dynamic states. The most powerful aspect of P. polycephalum as a model organism of behaviour lies in the direct link between actomyosin contractions, resulting in cytoplasmic flows and emerging behaviours. The broad understanding of the theory of active contractions (Bois et al., 2011; Radszuweit et al., 2013; Radszuweit et al., 2014; Julien and Alim, 2018; Kulawiak et al., 2019) might therefore well be the foundation to formulate the physics of behaviour not only in P. polycephalum but also in other simple organisms. This would not only open up an new perspective on life but also guide the design of bioinspired soft robots with a behavioural repertoire comparable to higher organisms.
Materials and methods
Experiments
View detailed protocolThe specimen was prepared from fused microplasmodia grown in a liquid culture (Daniel et al., 1962) and plated on 1.5%agar. The network was trimmed and imaged in the bright field setting in Zeiss ZEN two imaging software with a Zeiss Axio Zoom V.16 microscope equipped with a Hamamatsu ORCAFlash 4.0 digital camera and Zeiss PlanNeoFluar 1 x/0.25 objective. The acquisition frame rate was 3 sec. The stimulus was applied in a form of a heatkilled HB101 bacterial pellet in close network proximity.
Appendix 1
Data
The typical thickness of tubes in a P. polycephalum network is $\sim 50100\mu \mathrm{m}$ and the contraction amplitude about ~10% of the tube’s typical thickness (Alim et al., 2013). This change in tube thickness can be detected from a brightfield microscopy recording. We record one brightfield frame every three seconds. Since the periodic contractions of the tubes take place on the time scale of 100 sec, they are thus well resolved by the selected frame rate. Typically an idle network keeps a stable morphology and does not move significantly over a period of 1.5 h to 2.5 h which we use for recording its contraction dynamics. Since no two P. polycephalum specimens ever have the same network morphology, we are naturally constrained to one biological and one technical replicate in our experiments.
Data processing
Our data is a stack of brightfield images recorded from the P. polycephalum network with a rate of one frame every three seconds (Video 1; Video 2). Each brightfield frame has a time label t_{i} and the total number of frames is given by $T$. We process this data in the following steps. First, we mask the network in the brightfield images through thresholding. It is important to note that we use the same mask for all the images in the stack. This is possible since we consider a network that does not significantly move or change its morphology. This is true even when we apply a stimulus to the network, since we only consider the initial stages of stimulus response, before the network starts to display strong movement. From the masked regions of the brightfield frames, we extract pixel intensity values which we convert to 8bit format. Since we are here primarily interested in the contraction dynamics of the organism and not in the the actual base thickness of tubes or its longterm growth dynamics, we detrend the data using a movingaverage filter (rational transfer function) with a window size of two contraction periods (~ 200 sec) (Bäuerle et al., 2017). This leaves us only with the desired information about contractions taking place in the time scale of several minutes. We store the intensity values of each frame in a vector ${\overrightarrow{I}}^{{t}_{i}}$ of dimension $M$ equal to the number of pixels in the network, and $i$ indexes the frames in the range $i=1,\mathrm{\dots},T$. From the postprocessed data, we define the following data matrix
where $t$ denotes the matrix transpose.
Appendix 2
Principal component analysis (PCA)
The contraction modes are computed from the covariance matrix of the data. We compute the covariance matrix from the data matrix $\mathbf{X}$ after subtracting the mean from each column. The covariance matrix is given by
The sought after contraction modes ${\overrightarrow{\varphi}}_{\mu}$ are the eigenvectors of the covariance matrix
and ${\lambda}_{\mu}$ is the eigenvalue. The number of nonzero eigenvalues is equal to the rank of the covariance matrix. The eigenvalue captures the variance of the data along the direction of mode ${\overrightarrow{\varphi}}_{\mu}$. We also define the relative eigenvalue as
The mode coefficient ${a}_{\mu}$ is obtained by projecting the data onto mode ${\overrightarrow{\varphi}}_{\mu}$.
We note that we perform PCA on data segments with at least 700 frames (=35 min). Since it is well known from the literature (Kamiya, 1960) that the period of the contraction dynamics in P. polycephalum is on the order of 100 sec the analysed data contains on the order of 20 contraction periods at minimum. We can therefore be sure that we use enough data to resolve the characteristic dynamical features investigated here. As a further reassuring result, we recover the typical contraction period of 100 sec in our analysis, see Figure 1C(ii).
We add a brief comment on Fourier analysis, as alternative decomposition method to PCA. First, in one dimension, PCA is equivalent to Fourier decomposition. This is indeed apparent in our PCA analysis of the singletube data set where the principal components shown in Figure 5C correspond precisely to half a period of a sine and cosine Fourier mode. In two dimensions, the situation is more complicated. While we could in principle apply 2D Fourier decomposition we would need to apply Fourier analysis separately to every frame in our data set. However, this would mean that we have no information about the temporal evolution of mode activation. The Fourier modes would be different from one frame to the next and the activation of largescale patterns over time would be obscured.
Appendix 3
Distribution of temporal correlations
For a given time point t_{i}, the significant modes are determined based on the 70% criterion curve from Figure 2A. Next, the temporal correlations among the coefficients are computed in a time interval of ±15 frames around the time point t_{i}. The correlations are then counted in bins of the appropriate row of Figure 2A. Repeating this processing for all time points and normalising each row by the total number of correlations in that row, we obtain the final distribution shown.
Appendix 4
Flow rate calculation in a P. polycephalum cell with singletube morphology
To compute the flow rate of the cytoplasm in a P. polycephalum specimen with singletube morphology we use the theory developed in Shapiro et al., 1969; Li and Brasseur, 1993. In that work the flow of an incompressible Newtonian fluid inside an axisymmetric tube of fixed length is considered and the equations for the flow velocity field are written in the lubrication theory approximation. Furthermore, a timedependent thickness profile of longitudinal waves is imposed in the tube. Assuming noslip boundary conditions, the flow field can be fully determined at every point along the tube as a function of the timedependent tube profile. For the case when the tube profile is a periodic train of waves, we compute the volume flow rate averaged over an oscillation period by evaluating equation (13) of Li and Brasseur, 1993. We express the flow rate in units of volume of the entire tube divided by the oscillation period. This serves to characterize the performance in pumping of the significantly contracting P. polycephalum cell. We determine the time period over which to average the volume flow rate directly from the flowrate curve. Furthermore, the thickness profile of the tube is given by the measured pixel intensity profile.
Appendix 5
Mode superpositions in a P. polycephalum cell with singletube morphology
We are interested in how the contraction dynamics of the cell controls the cell’s locomotion behavior. In our analysis we therefore focus on the tube segment connecting the locomotion fronts at either end of the tube and perform Principal Component Analysis only on this part of the cell. Since the tube is effectively onedimensional, we find that the modes we obtain closely approximate Fourier modes. This means that superpositions of these modes afford a clear interpretation in terms of different contractionwave patterns. Such an interpretation is even further facilitated by the fact that we find that over large time intervals after the stimulus, the number of significant modes is very small. Indeed, over such time intervals it is sufficient to approximate the contraction dynamics with only one or two modes, as can be seen from Figure 5D. Hence we are essentially studying a superposition of modes ${\overrightarrow{\varphi}}_{1}$ and ${\overrightarrow{\varphi}}_{2}$ shown in Figure 5C and Figure 5—figure supplement 2 of the main text with their oscillating mode coefficients shown in Figure 5—figure supplement 3. To develop intuitive understanding of the nature of the superposition, we note that the modes ${\overrightarrow{\varphi}}_{1}$ and ${\overrightarrow{\varphi}}_{2}$ approximate sine and cosine functions over the length of the tube. Given a sine and cosine spatial contraction profile, different types of superpositions can be formed depending on the nature of their timedependent coefficients. To illustrate further, let us assume the idealised case where both coefficients are sine functions that can have different phases and amplitudes. Then, if the coefficient of one contraction profile is very small compared to the other, the resulting superposition is a standing wave. In case the coefficients have equal amplitudes, but are phase shifted by $\pi /2$, the superposition is a traveling waves. Finally, if the coefficient amplitudes are not equal and the phase shift lies somewhere between zero and $\pi /2$, the nature of the superposition is a mix of standing and traveling wave. Extrapolating this idealised picture allows us to infer the contraction dynamics resulting from our twomode approximation. We see that the coefficients of the two modes ${\overrightarrow{\varphi}}_{1}$ and ${\overrightarrow{\varphi}}_{2}$ shown in Figure 5—figure supplement 3 change in amplitude and phase relative to each other. It is easy to identify from this plot together with the plot of relative amplitudes in Figure 5D time intervals which approximate one of the contraction dynamics that we have described for the idealised system. Therefore we conclude that the superposition of the two top modes changes in its nature over time, ranging form a pure standing wave to a pure traveling wave.
Appendix 6
Choice of the cutoff of mode coefficient amplitudes
Our analysis of contraction dynamics requires us to place a cutoff on the amplitude of mode coefficients. Our chosen cutoff of 90% is supported in two ways:
First, the problem of choosing a cutoff for the coefficients is related to the problem of choosing a cutoff for the eigenvalue spectrum, since an eigenvalue is the variance of the mode coefficient. Given the continuous nature of the eigenvalue spectrum, there is no unique way to choose a cutoff. In (Berman et al., 2014), it is proposed to define the cutoff by the largest eigenvalue of the spectrum of the randomised data, see the black line in Figure 1B. We tested their criterion on our data and find a 93% cutoff, equivalent to roughly 70 modes. This is consistent with our choice of a 90% cutoff for the amplitudes.
Second, our main qualitative observation  considerable variation in the number of significant modes over time  is robust to different choices of cutoff values. In Figure 2A and similarly in Figure 3E, we show the number of significant modes for two different values of the cutoff, namely 70% and 90%.
Data availability
The two datasets from which Figures 1,2 and 3 and Figures 4 and 5 were generated are included as videos of raw brightfield time series in the article.
References

Scalefree networks in cell biologyJournal of Cell Science 118:4947–4957.https://doi.org/10.1242/jcs.02714

Coordinating theoretical and empirical efforts to understand the linkages between organisms and environmentsIntegrative and Comparative Biology 51:653–661.https://doi.org/10.1093/icb/icr091

Cooperative selforganization of microorganismsAdvances in Physics 49:395–554.https://doi.org/10.1080/000187300405228

Mapping the stereotyped behaviour of freely moving fruit fliesJournal of the Royal Society, Interface 11:20140672.https://doi.org/10.1098/rsif.2014.0672

Pattern formation in active fluidsPhysical Review Letters 106:028103.https://doi.org/10.1103/PhysRevLett.106.028103

Habituation in nonneural organisms: evidence from slime mouldsProceedings. Biological Sciences 283:1829.https://doi.org/10.1098/rspb.2016.0446

Hematinrequiring plasmodial myxomyceteJournal of Bacteriology 84:1104–1110.https://doi.org/10.1128/jb.84.5.11041110.1962

Mode Selection in Compressible Active Flow NetworksPhysical Review Letters 119:028102.https://doi.org/10.1103/PhysRevLett.119.028102

Variability in behavior and the nervous systemEncyclopedia of Human Behavior 4:447–458.

Peristaltic transport and mixing of cytosol through the whole body of Physarum plasmodiumMathematical Medicine and Biology 29:263–281.https://doi.org/10.1093/imammb/dqr010

Physics and Chemistry of Protoplasmic StreamingAnnual Review of Plant Physiology 11:323–340.https://doi.org/10.1146/annurev.pp.11.060160.001543

Dynamic organization of Physarum plasmodiumCell Motility and the Cytoskeleton 10:107–116.https://doi.org/10.1002/cm.970100115

A novel growth mode of Physarum polycephalum during starvationJournal of Physics D 51:244002.https://doi.org/10.1088/13616463/aac2b0

Coordination of contractility, adhesion and flow in migrating Physarum amoebaeJournal of the Royal Society, Interface 12:20141359.https://doi.org/10.1098/rsif.2014.1359

Nonsteady peristaltic transport in finitelength tubesJournal of Fluid Mechanics 248:129–151.https://doi.org/10.1017/S0022112093000710

Similarity in Neuronal Firing Regimes across Mammalian SpeciesThe Journal of Neuroscience 36:5736–5747.https://doi.org/10.1523/JNEUROSCI.023016.2016

Physarum polycephalum —a new take on a classic model systemJournal of Physics D 50:413001.https://doi.org/10.1088/13616463/aa8699

Intracellular mechanochemical waves in an active poroelastic modelPhysical Review Letters 110:138102.https://doi.org/10.1103/PhysRevLett.110.138102

Decisionmaking without a brain: how an amoeboid organism solves the twoarmed banditJournal of The Royal Society Interface 13:20160030.https://doi.org/10.1098/rsif.2016.0030

Variability in neural activity and behaviorCurrent Opinion in Neurobiology 25:211–220.https://doi.org/10.1016/j.conb.2014.02.013

Periodic traction in migrating large amoeba of Physarum polycephalumJournal of The Royal Society Interface 12:20150099.https://doi.org/10.1098/rsif.2015.0099

Migratory behaviour of Physarum polycephalum microplasmodiaThe European Physical Journal Special Topics 224:1199–1214.https://doi.org/10.1140/epjst/e2015024552

Patterns of cell thickness oscillations during directional migration of Physarum polycephalumEuropean Biophysics Journal 44:349–358.https://doi.org/10.1007/s0024901510287

BookDevelopmental Biology of PhysarumCambridge, United Kingdom: Cambridge University Press.

Peristaltic pumping with long wavelengths at low Reynolds numberJournal of Fluid Mechanics 37:799–825.https://doi.org/10.1017/S0022112069000899

Dimensionality and dynamics in the behavior of C. elegansPLOS Computational Biology 4:e1000028.https://doi.org/10.1371/journal.pcbi.1000028

Energysaving with low dimensional network in Physarum plasmodiumJournal of Physics D 50:154003.https://doi.org/10.1088/13616463/aa635a

Rhythmicity, recurrence, and recovery of flagellar beatingPhysical Review Letters 113:238103.https://doi.org/10.1103/PhysRevLett.113.238103

Coordination of eukaryotic cilia and flagellaEssays in Biochemistry 62:829–838.https://doi.org/10.1042/EBC20180029

Behavioral responses to changing environmentsBehavioral Ecology 26:665–673.https://doi.org/10.1093/beheco/aru183

ATP and calciumcontrolled contraction in a saponin model of Physarum polycephalumCell Structure and Function 9:135–141.https://doi.org/10.1247/csf.9.135
Decision letter

Herbert LevineReviewing Editor; Northeastern University, United States

Naama BarkaiSenior Editor; Weizmann Institute of Science, Israel
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
Thank you for submitting your article "Emergence of behavior in a selforganized living matter network" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
The paper on "The emergence of behavior …" has been seen by three referees all of whom have worked in the field of quantitative measures of organismal behavior. Based on their reports and my own reading of the paper, it is clear that the paper is not acceptable in its current form. As will be explained in detail below, there are many questions regarding exactly what is being claimed and proven regarding the modal analysis of the data. As there are no significant qualms about the experimental data per se and there is universal agreement that the questions being posed are important and the Physarum setting is a reasonable choice for addressing those questions, I have decided that the authors should be given a chance to address these questions. Only if these questions can be resolved by what undoubtedly would be a major revision would the manuscript be reconsidered for publication.
Essential revisions:
Without a more detailed presentation of the meaning of the various results shown for the eigenvalue spectrum and the model structure, it is very hard to assess the extent to which the observed changes are actually correlated with important behavioral responses. This is partially because of vague and unsupported statements such as (line 161) "This observation finally underlines that the continuous spectrum of modes and its variability in activation is intrinsic to the organism's behavior", and (line 215) "Our findings suggest that a continuous spectrum of contraction modes allows the living matter network P. polycephalum to quickly transition between a multitude of behaviors using the superposition of multiple contraction patterns."
Specific points that need to be addressed in this regard include:
– The discussion of the eigenvalue spectrum is severally lacking. While it is true that the authors don't find a lowdimensional system (although I am not sure why they expected to), they do find a powerlaw spectrum. There is a vast literature of graphs that have this property (including many from biology and ecology) and the authors should connect with this work. It is also true that purely random networks also have a continuous eigenvalue spectrum, although it is not a power law it can look like one sometimes. How well can the authors quantitatively determine that the observed data is a power law and not more representative of a random process. Is there some specific control dataset that can be generated by randomization of the original data to demonstrate the biological meaning of the features they find?
– It is difficult to judge if the findings from the tube explain the results for the networks. This is due to the fact that different aspects of the PCA are shown for different cases. Some problems are (i) is the symmetry breaking seen in the tube also present in the network and how can we see it (e. g. in the shape of the modes) and (ii) do the numbers of significant modes change in the idle tube in a similar way before signal application vary in a similar fashion than in the network.
– The authors should carefully bring out the relation between qualitative change of dynamics and the behavior (i. e. response to external stimulus) more clearly to justify title of their study. The authors seem to have shown that the behavior is realized by a change in contraction dynamics (one form of selforganisation) in an otherwise relatively constant random network topology (another form of selforganisation). Given that the selforganisation in the topology of Physarum has been studied abundantly in previous work (e. g. Nakagaki et al. in Science …), the novelty here is the claim that the dynamics correspond to different "behaviors" such as moving towards or away from a stimulus or enhanced pumping of optical fluid.
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Emergence of behavior in a selforganized living matter network" for further consideration by eLife. Your revised article has been evaluated by Naama Barkai (Senior Editor) and a Reviewing Editor.
All reviewers and I agree that the manuscript has been significantly improved but there are some remaining issues that need to be addressed, as outlined below:
As can be seen from the detailed reviews given below, there is still a somewhat mixed opinion about exactly how informative the mode analysis is regarding behavior. There are also some residual technical questions that could be answered in a modest revision. Given these reports and a positive recommendation from an initial third reviewer, I expect that the paper will be acceptable for publication once these questions are properly addressed.
Reviewer #1:
Overall, the manuscript has improved in revision and I appreciate the efforts of the authors. There are, however, still rough patches and I list some specific comments below.
I am still left wanting for more understanding of the continuous spectrum. Perhaps control is just spatially local? How dependent is the finding of highdimensionality on the particular representation of the behavior local contractions0? Or on the method (PCA)?
Line 16
I recommend removing "Surprisingly" as it is not clear why we would actually expect a lowdimensional covariance space. I would also remove an additional occurrence of "surprisingly" in the last paragraph of the introduction.
Line 26
The choice of Marom et al. (2002) as a reference is odd as this refers to neural activity in cultured networks, that is those without a body and thus without behavior.
Line 37
I would consider adding more recent references from work on ciliated organisms, e.g. work from the KY Wan and R Goldstein groups.
Line 103 How do we know that phi_4 is really preferentially exciting thicker tubes? A single image could simply be a spatial accident. Is there a correlation over the whole dataset?
Line 129 The temporal correlation of the modes should be better described. A first glance this is very confusing as the modes from PCA are, of course, uncorrelated across the whole dataset. Thus these correlations are coming from the window of + 15 frames. This should be clarified in the main text and/or the caption for Figure 2.
Line 151
Why do we jump from a description of Figure 3A immediately to Figure 3E and then later come back to BD?
Reviewer #3:
While the revised manuscript is improved in several ways, I remain unconvinced that the author's PCAbased analysis has revealed substantially new understanding about the behavior of P. Polycephalum. My main criticisms are the following:
1. As the authors mention in the paper, there are many ways to get the observed continuous power spectrum. E.g., pure 1/f noise would suffice. As such, this seems like a negative result, and I am not sure it adds to our understanding of the system.
2. The PCA modes appear to be approximately combinations of 2D plane waves, as is expected for a general 2D system. Indeed, for most 2D constructions without specific enhanced structures, a Fourierlike decomposition can be shown to be optimal. Given the increasing wavenumber for increasing mode number (as the authors defined using the power/eigenvalues), I think the data indicates that spatial frequency would be more informative compared to PCA.
3. The authors most interesting finding appears in Figure 3. relating the amplitudes in different modes and the number of "significant" modes to the application of a stimulus. I believe the authors do show that the response *is not* a simple 2state traveling wave (peristaltic pump), but the authors don't then follow up with a description/understanding of what the response actually *is*. Apparently, the organism engages many more modes to produce the behavior. Why does this combination require so many modes and how does it work? If the response is actually linear in shape (or step shaped, or anything with many frequencies), a Fourier decomposition would involve many modes (actually an infinite number). I struggle to see how the mode decomposition as presented has taught us how the slime mold moves or responds to stimuli.
https://doi.org/10.7554/eLife.62863.sa1Author response
Essential revisions:
Without a more detailed presentation of the meaning of the various results shown for the eigenvalue spectrum and the model structure, it is very hard to assess the extent to which the observed changes are actually correlated with important behavioral responses. This is partially because of vague and unsupported statements such as (line 161) "This observation finally underlines that the continuous spectrum of modes and its variability in activation is intrinsic to the organism's behavior", and (line 215) "Our findings suggest that a continuous spectrum of contraction modes allows the living matter network P. polycephalum to quickly transition between a multitude of behaviors using the superposition of multiple contraction patterns."
We thank the reviewers for pointing out the need for a more detailed presentation of the results.
We have now significantly modified the manuscript to better highlight how the changes in contraction dynamics correlate with behavioural responses, in particular with the most robust behavioural response in Physarum: creation of locomotion fronts upon food encounter. We modified the highlighted sentences and provided direct references to the figures depicting the described results.
Specific points that need to be addressed in this regard include:
– The discussion of the eigenvalue spectrum is severally lacking. While it is true that the authors don’t find a lowdimensional system (although I am not sure why they expected to), they do find a powerlaw spectrum. There is a vast literature of graphs that have this property (including many from biology and ecology) and the authors should connect with this work. It is also true that purely random networks also have a continuous eigenvalue spectrum, although it is not a power law it can look like one sometimes. How well can the authors quantitatively determine that the observed data is a power law and not more representative of a random process. Is there some specific control dataset that can be generated by randomization of the original data to demonstrate the biological meaning of the features they find?
We thank the reviewers for pointing out the need for additional theoretical discussion of the eigenvalue spectrum to supplement the already existing experimental analysis of the impact that network topology has on spectrum.
First, we would like to note that while we agree that the observed eigenvalue spectrum suggests a power law, we refrain from claiming so in the article. As already mentioned by the reviewers, a power lawlike eigenvalue spectrum can have many origins and is notoriously difficult to delineate, see e.g. M. Stumpf and M. Porter, (2012) “Critical Truths About Power Laws”, Science, Vol. 335, Issue 6069, pp. 665666. The most important observation for us is that the spectrum is continuous and has no obvious cutoff value. To further corroborate that the observed structure of the eigenvalue spectrum is the result of network specific contraction dynamics, we provide two additional analyses:
We created a control data set by randomising the original data along the temporal axis. Specifically, we independently randomize the data in each column (for each spatial point) of the data matrix and then perform PCA on this randomized data. The spectrum of the randomized data, see Figure 1—figure supplement 1, shows no remnants of a power law. The comparison between original and control data shows that the original eigenvalue spectrum truly reflects contraction dynamics that underlie processes of biological meaning: sensing nutrients, feeding, and locomotion. We now depict the mode distribution of the randomized data in the manuscript as Figure 1B.
To further exclude network topology as the source of the continuous eigenvalue spectrum, we generated simulated data of a traveling wave and imposed it onto the network. The PCA analysis reveals only two modes, with their coefficients shifted in phase by 90 degrees to superimpose into a traveling wave, see Author response image 1. Thus, spatial network complexity is not the cause of the continuous eigenvalue spectrum.
– It is difficult to judge if the findings from the tube explain the results for the networks. This is due to the fact that different aspects of the PCA are shown for different cases. Some problems are (i) is the symmetry breaking seen in the tube also present in the network and how can we see it (e. g. in the shape of the modes) and (ii) do the numbers of significant modes change in the idle tube in a similar way before signal application vary in a similar fashion than in the network.
We thank the reviewers for pointing out that we did not fully demonstrate the analogy between the dynamics within a network and a single tube.
In point (i), the reviewers refer to the phenomenon of ‘symmetry breaking’ in the creation of locomotion fronts. While we agree that it is tempting to use this term we prefer to characterise the observed phenomenon as a reversal of locomotion direction. The analogy between single tube and network holds independent of semantics. Namely, prior to the stimulus, the network shows the formation (growth) of a locomotion front at the bottom left corner. After the application of the stimulus, at the opposite edge of the network, the locomotion direction is reversed. We find that reversal of locomotion direction is correlated with the activation of distinct contraction patterns, see Figure 3 (poststimulus) and Figure 3—figure supplement 2 (prestimulus). Modes 2 and 3 in Figure 3 B are the drivers of locomotion front generation poststimulus. We also note the activation of mode 1 roughly 20 minutes after the stimulus application (see Figure 3 D). Notably, mode 1 has an almost identical spatial structure to mode 1 prestimulus shown in Figure 3—figure supplement 2B. The reactivation of this mode indicates that this contraction pattern is intrinsic to the network and is not simply erased by the stimulus.
The contraction dynamics associated with locomotion front reversal is analogous in the tube and the network: poststimulus, the original contraction pattern is dominated by a new mode, but with the possibility to switch back to the original pattern. The onedimensional nature of the single tube additionally allows us to develop a mechanistic understanding of locomotion in terms of standing and traveling wave contraction patterns.
In point (ii) the reviewers are asking whether there is a similarity between the dynamics of the number of significant modes in the single tube and in the network prestimulus. In Author response image 2 we show the relevant curve of both systems over equal length time intervals before the stimulus. During these intervals, both the single tube and the network experience slow, directed growth. In both cases the number of significant modes shows a considerable variation over time.
To demonstrate the analogy in contraction dynamics between network and single tube we now include new Figure 3 with Figure 3supplement 2, respectively, in the revised manuscript and we explain these results in section “Stimulus response behavior is paired with activation of regular, largescale contraction patterns interspersed by manymode states”.
– The authors should carefully bring out the relation between qualitative change of dynamics and the behavior (i. e. response to external stimulus) more clearly to justify title of their study. The authors seem to have shown that the behavior is realized by a change in contraction dynamics (one form of selforganisation) in an otherwise relatively constant random network topology (another form of selforganisation). Given that the selforganisation in the topology of Physarum has been studied abundantly in previous work (e. g. Nakagaki et al. in Science …), the novelty here is the claim that the dynamics correspond to different “behaviors” such as moving towards or away from a stimulus or enhanced pumping of optical fluid.
The reviewers are correct: we use the decomposition of tube contractions to show that behaviour is the result of a rich repertoire of contraction states. For the first time, the attention is drawn away from the topological changes in the network and instead focused on the processes that drive them. In Figure R3 and R4 we correlate contraction dynamics and measured growth behaviours, finding distinct contraction patterns for different locomotion front states, before and after stimulus application. Notably, our analysis of network contractions shows different behaviours within an almost constant network topology. Specific behaviours are associated with a small number of dominant modes (see the different growth fronts in the network) while transitions between behaviours are marked by simultaneous activation of many modes. Our theoretical analysis using the spectrum of randomised data in Figure 1—figure supplement 1, further shows that not network complexity but contraction dynamics are key to behaviour.
We thank the reviewers for recognizing and summarizing the message of the article, which we now made clearer by including a quantitative presentation demonstrating the correspondence between contraction dynamics and behaviour for the network shown in updated Figure 3 with Figure 3supplement 2.
References:
1. Stumpf, M. and Porter, M., (2012) Critical Truths About Power Laws. Science, Vol. 335, Issue 6069, pp. 665666.
2. Berman GJ, Choi DM, Bialek W, Shaevitz JW (2014) Mapping the stereotyped behaviour of freely moving fruit flies. J R Soc Interface 11(99). Doi:10.1098/rsif.2014.0672.
3. Julien, J.D. and Alim, K. (2018) Oscillatory fluid flow drives scaling of contraction wave with system size. PNAS October 16, 2018 115 (42) 1061210617.
[Editors’ note: further revisions were suggested prior to acceptance, as described below.]
Reviewer #1:
Overall, the manuscript has improved in revision and I appreciate the efforts of the authors. There are, however, still rough patches and I list some specific comments below.
I am still left wanting for more understanding of the continuous spectrum. Perhaps control is just spatially local? How dependent is the finding of highdimensionality on the particular representation of the behavior local contractions0? Or on the method (PCA)?
We thank the reviewer for their thoughtful feedback. In our elaborate answer in the following we first present additional data on the dynamics of the continuous spectrum before we discuss the idea of localized control and what governs the highdimensionality of the spectrum, before we finally justify our choice of method.
To address the continuous nature of the spectrum, it is important to note that the contraction dynamics of P. polycephalum is not a sharp switching between a few discrete largescale contraction states. For the contraction dynamics to change smoothly as observed by us, a continuous spectrum is required allowing for contractions over a smooth spatial scale to be activated. In Figure 3—figure supplement 4, we show the instantaneous rank of the top 80 modes over time for the unstimulated network. The instantaneous ranking is based on the relative amplitudes (see Equation (2) in the article) and serves as a measure of activity for a mode. We observe that modes frequently change their instantaneous rank, moving across a wide range of ranks which is only possible in a continuous spectrum.
The concept of localised control, as suggested by the reviewer, appears to be a question regarding the underlying mechanism by which contractions are organised across the network. We emphasize that the spectrum we observe is the immediate result of the observed phenomenology of contractions and not of the mechanism by which contractions are coordinated. The continuous spectrum of longrange patterns that we observe may result from a number of different mechanisms, for instance flowbased transport of signaling molecules [2] or hydrodynamic coupling, yet it is beyond the scope of this work to discern between these.
To address the highdimensional nature of the contraction dynamics we start with the reviewer’s suggestion that local contractions are the cause of high dimensionality. To discuss this, we show in Figure 1B a modified version of a plot that was included in our previous response to the reviewers. In this figure we compare the eigenvalue spectrum of the original data to the spectrum of randomised data. As indicated in Figure 1B, the spectrum of the randomised data (gray), defines an upper noise bound (red), separating the spectrum into large eigenvalues and small eigenvalues which include noise. Before discussing the possibility of local contractions, it is useful to reiterate the interpretation of the large eigenvalues. These eigenvalues correspond to contraction modes with a large spatial scale (see e.g. Figure 3C). The patterns of these modes reflect the network’s locomotion behaviour in response to the stimulus and targeted activation of thick tubes in the network (also see [1] for a discussion of behaviour on longer time scales). Notably, the number of large eigenvalues is substantial (~80). Thus, the contraction dynamics is highdimensional even when one disregards the part of the spectrum that lies below the upper noise bound. We now point to the continuous spectrum of large eigenvalues in the revised manuscript. While it is conceivable that one can build a model of purely local contractions which yields a continuous eigenvalue spectrum, such a model would not lead to modes showing the largescale patterns compatible with network morphology and behaviour as we observe them. However, a contribution from local contractions to the spectrum that we observe is not excluded. For instance, it is possible that the small eigenvalues (lowranked modes) in Figure 1B, are associated with localised contractions, for instance a single tube in the network. The interpretation of such local contractions is however strongly dependent on the network’s specific morphology and local state. We emphasize again that the presence of local contractions would not change the fact that we consistently observe longrange patterns across the network.
Finally, we address the role of the choice of decomposition method for the finding of high dimensionality. In one dimension, PCA is equivalent to Fourier decomposition (regarding this point, see also our reply to reviewer #3’s second question further below). For the single tube we are thus effectively performing Fourier decomposition and find a continuous spectrum (see Figure5—figure supplement 1). In two dimensions there are other decomposition methods which differ from PCA by the assumptions made for the decomposition basis. Changes of basis typically do not affect the dimensionality, and PCA is in fact often used as a preanalysis step that feeds into another decomposition method to change the basis. Therefore, amongst the standard 2D decomposition methods the finding of high dimensionality holds qualitatively. We provide a discussion on the choice of method in the revised Appendix 2.
Line 16
I recommend removing “Surprisingly” as it is not clear why we would actually expect a lowdimensional covariance space. I would also remove an additional occurrence of “surprisingly” in the last paragraph of the introduction.
We have removed both occurrences of the word “surprisingly”. In the abstract (line 16), we have replaced it by “notably” to emphasise that it is central to our work.
Line 26
The choice of Marom et al. (2002) as a reference is odd as this refers to neural activity in cultured networks, that is those without a body and thus without behavior.
We are grateful to the reviewer for pointing this out. We agree that this choice of reference is not an optimal one given the context of the study. In the revised version, we removed this reference and replaced it with a reference [Mochizuki et al.,2016], providing an overview of neuronal activity in organisms exhibiting behaviour.
Line 37
I would consider adding more recent references from work on ciliated organisms, e.g. work from the KY Wan and R Goldstein groups.
We thank the reviewer for their suggestion. We have added [Wan and Goldstein, 2014] and [Wakefield et al., 2018] as new references.
Line 103 How do we know that phi_4 is really preferentially exciting thicker tubes? A single image could simply be a spatial accident. Is there a correlation over the whole dataset?
We thank the reviewer for their question. Visual inspection of brightfield frames indeed shows that thicker tubes are a dominant feature over long time intervals of the data set. Furthermore, principal components are ranked according to their variance (the eigenvalues) which is given by the time average of the square of the mode amplitude. A mode which is activated only over a short time interval is thus very unlikely to receive a ranking as high as mode phi_4.
Line 129 The temporal correlation of the modes should be better described. A first glance this is very confusing as the modes from PCA are, of course, uncorrelated across the whole dataset. Thus these correlations are coming from the window of + 15 frames. This should be clarified in the main text and/or the caption for Figure 2.
We thank the reviewer for making us aware of this potentially confusing point. The answer is that the PCA modes (for example the network modes shown in Figure 1C, or the tube modes shown in Figure 5C) are spatially uncorrelated, however their activation in time can certainly be correlated. This is what we quantify in Figure 2 and we choose windows of +15 frames to see how this correlation changes over time. Following the reviewer’s suggestion we have clarified this point in the text and the caption of Figure 2.
Line 151
Why do we jump from a description of Figure 3A immediately to Figure 3E and then later come back to BD?
We thank the reviewer for pointing this out and we agree with them that this is a clumsy presentation. We have also noticed the same issue in Figure 5. In the revised manuscript, we rearranged the order of subplots, adjusted the letter labels in Figure 3 and Figure 5, adapted the figure captions and ensured correct referencing across the article.
Reviewer #3:
While the revised manuscript is improved in several ways, I remain unconvinced that the author's PCAbased analysis has revealed substantially new understanding about the behavior of P. Polycephalum. My main criticisms are the following:
1. As the authors mention in the paper, there are many ways to get the observed continuous power spectrum. E.g., pure 1/f noise would suffice. As such, this seems like a negative result, and I am not sure it adds to our understanding of the system.
We thank the reviewer for emphasising the point that a continuous spectrum can have multiple sources, including certain types of noise. With reference to our reply to reviewer #1’s first question given above (see also Figure 1B), we believe that there is enough evidence that a substantial number of the highranked modes are representative of large scale contraction dynamics of the network. However, this does not exclude the possibility that part of the spectrum (in particular the lower ranked modes with smaller eigenvalues) correspond to noise and possibly some form of local contractions. In fact, the possibility of noise playing an active role in shaping the network’s behaviour to us appears like an interesting feature, worthy of a future study. Regarding this point, please also see our reply to your third point below.
2. The PCA modes appear to be approximately combinations of 2D plane waves, as is expected for a general 2D system. Indeed, for most 2D constructions without specific enhanced structures, a Fourierlike decomposition can be shown to be optimal. Given the increasing wavenumber for increasing mode number (as the authors defined using the power/eigenvalues), I think the data indicates that spatial frequency would be more informative compared to PCA.
We thank the reviewer for their comment regarding Fourier decomposition as an alternative to decomposition into principal components. While the reviewer is specifically asking for the 2D case, we would first like to emphasize that in 1D, Principal Component Analysis is equivalent to Fourier decomposition. This is indeed apparent in our PCA analysis of the single tube data set where the principal components shown in Figure 5C correspond precisely to half a period of a sine and cosine Fourier mode.
In 2D, the situation is more complicated for multiple reasons. First, applying Fourier decomposition makes sense when there is a clear periodic structure visible in the data. However, from the brightfield movies we do not see spatially periodic patterns that are longterm stable. Furthermore, there are in fact enhanced structures given by the network’s morphology, such that we cannot expect waves to propagate homogeneously across the network. However, for us the most important reason for employing PCA is the ability to resolve temporal dynamics. While we could in principle apply 2D Fourier decomposition (setting our previous reasons for not doing so aside for the moment), we would need to apply Fourier analysis separately to every frame in our data set. This would mean however, that we have no information about the temporal evolution of mode activation. The Fourier modes would be different from one frame to the next and the activation of large scale patterns over time would be obscured.
For the reasons that we listed above, Principal Component Analysis is a widely used approach to analysing 2D systems with temporal dynamics. This is even the case for systems that have a simple geometry such as C. elegans with clear Fourierlike modes, where however the most important aspect of the analysis is how the coupling of modes evolves over time [3]. We added an explanation of the advantage of PCA over Fourier decomposition to the revised Appendix 2.
3. The authors most interesting finding appears in Figure 3. relating the amplitudes in different modes and the number of "significant" modes to the application of a stimulus. I believe the authors do show that the response *is not* a simple 2state traveling wave (peristaltic pump), but the authors don't then follow up with a description/understanding of what the response actually *is*. Apparently, the organism engages many more modes to produce the behavior. Why does this combination require so many modes and how does it work? If the response is actually linear in shape (or step shaped, or anything with many frequencies), a Fourier decomposition would involve many modes (actually an infinite number). I struggle to see how the mode decomposition as presented has taught us how the slime mold moves or responds to stimuli.
We thank the reviewer for their question regarding the interpretation of contractions after stimulus application.
We briefly summarise the relevant observations from Figure 3. What we present in Figure 3D is that over the course of ~25min after stimulus application, a sequence of dominant contraction patterns is activated. These dominant patterns are given by the three highest ranked modes. Starting at minute 85, first mode 3 (orange line), then mode 2 (blue), followed by mode 1 (red), and finally again mode 2 are activated. The structure of these modes is shown in Figure 3C. The activation of these dominant patterns is interspersed by shorter time intervals where a large number of modes are roughly equally active. Additionally, in Figure 3B we quantify the behaviour of the network and find upwards locomotion behaviour. Finally, we also note that we make the analogous observations for the single tube presented in Figure 5. Thus, just like for the single tube, we find largescale contraction patterns also for the network after stimulus application. However, due to their 2D nature and the complex network morphology, these patterns are not just simple sine and cosine waves as explained in previous point. Nevertheless, activation of specific oscillation modes can be directly linked with behaviour. This is illustrated for the activation of mode 3 by the pink box extending across subfigures Figure 3B and D. A detailed quantification of the flow rates in the network associated with specific contraction patterns is hard in complex morphologies and well beyond the scope of our study. However, given the structure of the modes one can expect them to act like pumps generating mass redistribution.
The situation is similar for the states with many significant modes marking transitions from one dominant pattern to the next. In the single tube we are actually able to correlate the number of significant modes with the magnitude of the flow rate in the tube which governs mass redistribution, see Figure 4. A large number of significant modes is linked with a small flow rate and thus small mass redistribution. In this state, modes do not add up to give a regular largescale contraction pattern, but their superposition yields irregular patterns with a short length scale. While a similar quantification for irregular contraction states in the network is beyond the scope of this study, we believe that one of the strong points of our work is to reason by analogy with the single tube that the many mode states do not lead to mass redistribution. Additional evidence for this is given in Figure 2, where we show that when there are many significant modes, the modes are not strongly correlated in time, again suggesting that their superposition does not produce the regular largescale patterns required for a large flow rate.
We believe that the core message of our finding is that many modes are employed to produce a broad range of behaviours, where transition between behaviours is enabled by adjustment of mode amplitude. We hope that we have now illustrated this key finding by revising Figure 3 and Figure 5 as outlined in our reply to reviewer #1’s last comment above.
References:
1. Kramar, M. and Alim, K.. Encoding memory in tube diameter hierarchy of living flow network. Proceedings of the National Academy of Sciences. 2021; 118 (10).
2. Julien, J.D. and Alim, K. (2018) Oscillatory fluid flow drives scaling of contraction wave with system size. PNAS October 16, 2018 115 (42) 1061210617.
3. Stephens, GJ., JohnsonKerner, B., Bialek, W., Ryu, WS. Dimensionality and Dynamics in the Behaviour of C. elegans. PloS Computational Biology. 2008; 4(4):110.
https://doi.org/10.7554/eLife.62863.sa2Article and author information
Author details
Funding
Simons Foundation (400425)
 Philipp Fleig
IMPRS for Physics of Biological and Complex Systems
 Mirna Kramar
Max Planck Society
 Philipp Fleig
 Mirna Kramar
 Michael Wilczek
 Karen Alim
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Senior Editor
 Naama Barkai, Weizmann Institute of Science, Israel
Reviewing Editor
 Herbert Levine, Northeastern University, United States
Publication history
 Received: September 7, 2020
 Preprint posted: September 8, 2020 (view preprint)
 Accepted: December 20, 2021
 Version of Record published: January 21, 2022 (version 1)
Copyright
© 2022, Fleig et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,118
 Page views

 108
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Medicine
 Physics of Living Systems
Internephron interaction is fundamental for kidney function. Earlier studies have shown that nephrons signal to each other, synchronize over short distances, and potentially form large synchronized clusters. Such clusters would play an important role in renal autoregulation, but due to the technological limitations, their presence is yet to be confirmed. In the present study, we introduce an approach for highresolution laser speckle imaging of renal blood flow and apply it to estimate the frequency and phase differences in rat kidney microcirculation under different conditions. The analysis unveiled the spatial and temporal evolution of synchronized blood flow clusters of various sizes, including the formation of large (>90 vessels) and longlived clusters (>10 periods) locked at the frequency of the tubular glomerular feedback mechanism. Administration of vasoactive agents caused significant changes in the synchronization patterns and, thus, in nephrons’ cooperative dynamics. Specifically, infusion of vasoconstrictor angiotensin II promoted stronger synchronization, while acetylcholine caused complete desynchronization. The results confirm the presence of the local synchronization in the renal microcirculatory blood flow and that it changes depending on the condition of the vascular network and the blood pressure, which will have further implications for the role of such synchronization in pathologies development.

 Cell Biology
 Physics of Living Systems
In addition to diffusive signals, cells in tissue also communicate via long, thin cellular protrusions, such as airinemes in zebrafish. Before establishing communication, cellular protrusions must find their target cell. Here, we demonstrate that the shapes of airinemes in zebrafish are consistent with a finite persistent random walk model. The probability of contacting the target cell is maximized for a balance between ballistic search (straight) and diffusive search (highly curved, random). We find that the curvature of airinemes in zebrafish, extracted from livecell microscopy, is approximately the same value as the optimum in the simple persistent random walk model. We also explore the ability of the target cell to infer direction of the airineme’s source, finding that there is a theoretical tradeoff between search optimality and directional information. This provides a framework to characterize the shape, and performance objectives, of noncanonical cellular protrusions in general.