An Active Poroelastic Model for Mechanochemical Patterns in Protoplasmic Droplets of Physarum polycephalum

Motivated by recent experimental studies, we derive and analyze a two-dimensional model for the contraction patterns observed in protoplasmic droplets of Physarum polycephalum. The model couples a description of an active poroelastic two-phase medium with equations describing the spatiotemporal dynamics of the intracellular free calcium concentration. The poroelastic medium is assumed to consist of an active viscoelastic solid representing the cytoskeleton and a viscous fluid describing the cytosol. The equations for the poroelastic medium are obtained from continuum force balance and include the relevant mechanical fields and an incompressibility condition for the two-phase medium. The reaction-diffusion equations for the calcium dynamics in the protoplasm of Physarum are extended by advective transport due to the flow of the cytosol generated by mechanical stress. Moreover, we assume that the active tension in the solid cytoskeleton is regulated by the calcium concentration in the fluid phase at the same location, which introduces a mechanochemical coupling. A linear stability analysis of the homogeneous state without deformation and cytosolic flows exhibits an oscillatory Turing instability for a large enough mechanochemical coupling strength. Numerical simulations of the model equations reproduce a large variety of wave patterns, including traveling and standing waves, turbulent patterns, rotating spirals and antiphase oscillations in line with experimental observations of contraction patterns in the protoplasmic droplets.


Introduction
The true slime mold Physarum polycephalum is an extensively studied system in biophysics. The plasmodial stage is of particular interest, since it exhibits, despite the relatively simple organization of this unicellular organism, seemingly ''intelligent'' physiological processes [1]. In this context the term ''intelligent'' means that, given an external stimulus, the plasmodium optimizes its cell shape, vein network and growth with respect to transport efficiency as well as robustness with respect to link deletion and avoidance of unfavorable conditions [2,3]. Recent experiments along these lines show that plasmodia were able to reproduce public transport networks on the scale of a petri dish [4] and to ''solve'' maze problems such as finding the shortest path between two food sources placed at the exits of a labyrinth [5]. Several groups have also investigated the topology and dynamical evolution of the vein network in large Physarum plasmodia with graph theoretical and statistical physics tools [6][7][8]. A second remarkable phenomenon is the synchronization of the contraction patterns in the tubular vein network that generates shuttle streaming to distribute nutrients efficiently throughout the organism [9]. From the perspective of biophysics it is natural to consider these phenomena in the framework of self-organized complex systems [10]. For the formulation of mathematical models a basic understanding of chemical and mechanical processes in the protoplasm is needed.
A first model for strand contraction combined the viscoelastic properties of the ectoplasmic wall with a reaction kinetics that regulates the contractile tension of the actomyosin system [11,12]. Later, several models in the form of reaction-diffusion [13] and reaction-diffusion-advection equations [14,15] were formulated that use homogenized quantities, for instance the average strand thickness. These models describe Physarum protoplasm as an oscillatory medium and treat the mechanical feedback in a simplified, qualitative way. More realistic models consider, instead, a two-phase description that distinguishes a fluid sol ( = cytosol) and a solid gel ( = cytoskeleton) phase. Some of these models account for sol-gel transformations and were used to explain flowchannel formation [16] and front dynamics [17].
Experiments with microplasmodia, i.e. small plasmodia of sizes ranging from 100 mm to several millimeters, provide a possibility to study internal amoeboid dynamics of Physarum without the pronounced vein structures usually present in Physarum cells of larger size. Such microplasmodia are produced by extracting cytosol from a Physarum vein and placing it on a substrate. Given a sufficient amount of cytosol, protoplasmic droplets will reorganize and form a new independent cellular entity. During the first hours of this process such cells show a surprising wealth of spatiotemporal mechanical contraction patterns [18,19] (and U. Strachauer & M.J.B. Hauser, unpublished data). The fact that the cell morphology does not change dramatically and that the cell does not migrate during the first hours, permits observation of the mechanical deformation patterns and waves in a quasi-stationary setting. The observed patterns include spirals, traveling and standing waves as well as antiphase oscillations (see Fig. 1).
Various patterns were reproduced previously by a qualitative particle-based model [20]. However, this description provided no information about the mechanical quantities that are essential to understand the intracellular deformation waves and patterns seen in the experiments.
In a more general context, studying the spatiotemporal instabilities and the related symmetry breaking in intracellular processes has become important to understand many biological processes. In a pioneering paper [21], Turing suggested that the interplay of reactions and diffusion processes provides a fundamental mechanism for morphogenesis. Later, models for intracellular pattern formation that include mechanical forces and the resulting advection processes have been suggested [22,23]. Active gel models describe the cytoskeleton as an active viscous fluid [24]. In contrast, experiments on inhomogeneous hydration in cells, where large pressure gradients in the cell are observed [25] indicate that the cytoplasm can behave like a porous elastic sponge-like solid (cytoskeleton) penetrated by a viscous fluid phase (cytosol) [26,27]. Moreover, several multiphase flow models have been proposed as appropriate descriptions of cytoplasmic dynamics [28,29].
In this paper, we derive and investigate a poroelastic two-phase model of the cytoplasm assuming a viscoelastic solid phase (cytoskeleton) and a fluid phase (cytosol). Furthermore, we incorporate an active tension in the solid phase which is regulated by the concentrations of free calcium ions in the fluid phase, that are in turn advected by the cytosolic flow field. To account for the calcium oscillations observed in experiments with Physarum, a simple one-dimensional active poroelastic model derived earlier [30] is extended to two dimensions and supplemented by a coupling to an oscillatory reaction-diffusion dynamics of the intracellular calcium concentration [31]. Unlike the regulating species in [30], the total free calcium concentration in the model presented here is not conserved anymore and varies strongly in time. As a result, the model displays a short-wavelength instability as opposed to the long-wavelength instability found in [30]. The choice of the poroelastic approach is motivated by the fact that the typical oscillation period of 1 -2 minutes connected with the spatiotemporal contraction patterns discussed above is considerable shorter than the experimentally observed time of 3 -6 minutes at which the cytoskeleton starts to exhibit fluid behavior [32]. Hence, the resulting model describes the cytoskeleton as an active viscoelastic solid coupled to a passive fluid in contrast to earlier works that had modeled the cytoskeleton itself as an active fluid [22] addressing long time scales, for which fluidization of the cytoskeleton has already occurred.
The inclusion of the calcium oscillator is necessary because it is known to be essential in the regulation of the contractile actomyosin system of Physarum polycephalum [33]. The mechanical and biochemical model parameters are mostly taken from the experimental literature on Physarum and therefore allow quantitative comparisons between model results and experiments. Altogether, in this article we derive and analyze a model for the intracellular dynamics of protoplasmic droplets that treats the cellular mechanics in the framework of a continuous two-phase active poroelastic model coupled to an oscillatory biochemical medium. Active mechanical stresses induce internal pressure gradients that generate cytosolic flow. To account for the latter we extend the equations for the internal calcium oscillator proposed in [31] to a reaction-diffusion-advection (RDA) model closing the mechano-chemical feedback loop. In the methods section we introduce and derive the mathematical mode with a description divided into a mechanical and a biochemical part. Subsequently, the choice of the model parameters is discussed and the numerical methods used to discretize and solve the resulting partial differential equations (PDEs) are described. The next section contains the results obtained by linear stability analysis at the homogeneous steady state (HSS) and a two-parameter phase diagram with numerical simulations. We present also a selection of qualitatively different contraction patterns obtained from simulations of our model, compare them to earlier experimental findings and demonstrate that the variety of patterns found in the experiments with Physarum droplets is reproduced successfully. The paper is concluded with a discussion of the presented model, its limitations and possible extensions.

Model: Mechanical part
Physarum protoplasm contains a rudimentary form of an actomyosin system that is also present in cells of higher vertebrates. In contrast to muscle cells the actin filaments in Physarum are randomly oriented in the cortex [34,35]. In our model, we assume that the cytoplasm contains a solid filamentous phase (gel phase) that has viscoelastic properties and exerts contractile tension on the system. This is illustrated schematically in Fig. 2. The derivation given below is analogous to a recently published generic model for active poroelastic media [30]. Unlike the regulating species in Ref. [30], the total free calcium concentration in the model presented here is not conserved anymore and varies strongly in time. The fluid part of the cytoplasm is modeled as a passive fluid (sol phase) that permeates the cytoskeleton [25]. The velocity field in the sol phase will be expressed by the variable v. Typical Reynolds numbers that arise from the cytoplasmic flow are small (Re%1 ) and inertia is negligible. The relative volume fraction of solid material is denominated as r gel and the fraction of  [18,19]. The color represents the local phase of oscillation obtained by a Fourier transformation of the spatiotemporal height data: a) standing wave, b) many irregular spirals, c) traveling wave, d) antiphase patterns, and e) single spiral. doi:10.1371/journal.pone.0099220.g001 the fluid material as r sol with the additional constraint r sol zr gel~1 .
We define a body-reference coordinate system x and a displacement field u that gives the deviation of the deformed coordinates X: u(x,t)~X(x,t){x at a time t. The gel velocity is the substantial time derivative _ u u~L t uz(=u) : _ x x of the displacement field. Since the gel is fixed in the reference coordinate system ( _ x x~0), the substantial time derivative _ u u is identical to the partial time derivative L t u.
To determine the flow and displacement field we consider forcebalance equations of the form = : (r gel (s pass: gel zT a 1))zf gel~0 ð1Þ = : (r sol s pass: where the passive sol and gel stresses s pass: sol=gel are determined by linear constitutive laws and f sol=gel are force densities. The active stress generated in the bulk of the gel phase is expected to be isotropic and hence, given by a multiple of the unit tensor: s act: gel~T a 1. We assume only small deformations and thus restrict ourselves to linear elastic theory. We consider the gel phase as a porous viscoelastic active material [30] that is able to exert contractile stresses by interaction of the myosin-motor system with the actin filaments. The filament orientation in the cortex of Physarum that mainly determines the active and passive properties of the medium is random [35]. Hence, we use an isotropic constitutive law for the elastic stress-strain relation. It has been suggested to consider the cytoplasm as an incompressible medium [36]. In the three-dimensional bulk the total mass flux must be zero. For small strains this can be expressed as [29]: In the following we assume constant sol and gel fractions r 0 sol and r 0 gel throughout the medium. This is justified, because we consider only small deformations. As a result the transport of cytosol and the related potential inhomogeneities of the fields r sol and r gel lead only to second-order corrections in the mechanical equations (for details see Text S1 in File S2).
We include a hydrostatic pressure p into the stress tensors of sol and gel that originates from the incompressibility of the material expressed by Eq. (3) and neglect the osmotic pressures caused by a difference in the chemical potential (see Text S1 in File S2). We assume a Kelvin-Voigt viscoelastic constitutive law (see e.g. [37]) for the gel phase. Using Darcy's law w~{ k g =p the following relation for the drag force is obtained where the parameter b is the ratio between the dynamic viscosity g of the cytosol and the permeability k of the porous medium.
Instead of the fluid velocity v in the laboratory frame, we need to consider the velocity w~v{ _ u u in the body-reference frame introduced above. The sol phase is considered as a passive Newtonian fluid. Hence, only viscous stresses are included for the sol phase. As a result, Eq. (2) corresponds to the Brinkman equation [38]. With the usual expression w~{ k g =p for Darcy's law and Eq. (2) one can relate the coefficient b in our model to the viscosity g of the cytosol and the permeability k of the medium: br 0 gel~g =k. We divide the overall passive stresses in a dissipative and nondissipative (elastic) part: The factors g shear sol=gel in front of the symmetric traceless parts denote the shear viscosities and g bulk sol=gel the bulk viscosities for sol and gel phase, respectively. In the elastic stress G is the shear modulus and K the compression modulus. The mechanical force balance equations in the final form (for details, see Text S1 in File S2) together with the incompressibility condition read: Model: Chemical part Calcium ions play a key role in the regulation of contraction in Physarum polycephalum. To describe the calcium kinetics, we use the biochemically realistic model presented by Smith and Saldana [31]. This model describes an oscillation mechanism that is driven by a phosphorylation-dephosphorylation cycle of myosin light chain kinases (MLCKs). A crucial feature of the Smith-Saldana model are autonomous calcium oscillations that occur already in the absence of mechanical feedback. Another entirely different mechanism was introduced to explain oscillations of plasmodial strands of Physarum, where the necessary feedback for the oscillations is provided by mechano-sensitive channels. This feedback acts upon a non-oscillatory calcium reaction kinetics and is therefore required to obtain the oscillations [12,39]. This model, however, predicts that the free calcium concentration and the mechanical tension oscillate in phase, whereas experiments exhibit an antiphase oscillation with a phase shift of p between these two quantities [33]. Moreover, experiments in the homogenate of Physarum plasmodium wherein deformations and mechanical stresses are not possible yield calcium oscillations giving further evidence for an autonomous calcium oscillator [40]. The Smith-Saldana model [31] can be reduced to two ordinary differential equations involving the free calcium concentration n c and the fraction of phosphorylated myosin-light-chain kinases w: These equations involve the myosin-bound calcium concentration and the calcium-dependent phosphorylation rate of the MLCK The chemical parameters introduced in Ref. [31] include the concentration of myosin N M , the equilibrium calcium concentration N c , a pumping and leaking rate k V and k L of calcium vacuoles, and the dephosphorylation rate of the MLCK k E . Furthermore, we have the affinity K a,b of calcium to the LC1 binding site on the myosin light chain (MLC). The value of this affinity depends on whether the LC2 binding site is phosphorylated (b) or not (a). The AC-cAMP-PKA signal cascade is modeled by Eq. (9) using a Hill-equation with coefficients b, K ? , and k 0 Q . For a detailed explanation of the terms in Eq. (7), we refer the reader to the publications [31,41]. The model is completed by a third equation that relates the activated fraction of myosin to the active tension T a . In contrast to the original work of Saldana and Smith, we introduce a relaxation equation analogous to models of cardiac myocytes [42]: where t T is the relaxation time the tension needs to approach its equilibrium value and F T represents the mechanochemical coupling strength. According to the model in Ref. [31] the fraction of activated myosin motors (available for binding to actin) is given by In this equation two additional parameters, the phosphorylation rate k P and dephosphorylation rate k D of the LC1 binding site come into play. The relaxation time constant t T is introduced to obtain a phase shift between calcium concentration and active tension that is p. Earlier, a direct proportionality was used [31]. But this work could not reproduce the phase shift between calcium oscillation and active stress that was observed in experiments. Thus, we have adjusted t T to obtain the correct phase relation.
The next step in the derivation of the model is a spatial extension of Eq. (7) to a reaction-diffusion-advection system. The only species in this model that is transported by diffusion is the free calcium n c . Altogether, the following equation is obtained where w is the fluid velocity in the body-reference frame introduced above and D c denotes the diffusion constant of the free calcium in the cytosol.

Summary of the model
Collecting all pieces of the model described so far, we obtain the following two-dimensional system of partial differential equations: 0~= : Here, we have introduced the shear and bulk viscosities g shear sol=gel and g bulk sol=gel of sol and gel phase and the linear elastic shear and compression modulus of the gel phase K and G, respectively. Note that the sol fraction and the free calcium concentration are given in the body-reference frame. The pressure p formally appears as a Lagrangian multiplier, due to the incompressibility condition (like in the incompressible Navier-Stokes equation).
The above equations are defined in a circular domain that mimics the geometry of the Physarum droplets in experiments [18,19] (and U. Strachauer & M.J.B. Hauser, unpublished data). Zero displacement and zero cytosolic flow are imposed at the boundaries. Assuming, a non-permeable membrane of the droplet, no-flux boundary conditions =n c (x,t) : n(x,t)~0 are employed for the calcium concentration n c , where n is the normal vector at the boundary. As initial conditions we apply a small random noise j with zero mean as a perturbation to the homogeneous steady state solution.
To compare with the experimentally measured height profiles of the droplet, one needs to estimate the local height field H(x,t) that does not appear explicitly in the two-dimensional model. We relate the relative height deviation h(x,t) :~(H(x,t){H 0 )=H 0 (where H 0 is the height in the reference configuration) to the divergence of the displacement field computed in the two-dimensional model: This approximation is based on the idea that deformations are locally isotropic [35]. Though Eqs. (13)-(18) may appear quite complex at first glance, their essential aspect is the combination of the mechanical dynamics of a poroelastic medium and the calcium oscillator.

Parameters
A summary of the parameters used in the model is given in Tables 1 and 2. The values for the chemical oscillator are taken from Ref. [31]. For the affinity K a that appears in the functions f c and h in Eq. (7) two different values are considered. This parameter represents the affinity of the myosin light chain kinase to calcium ions in the unphosphorylated state. Autonomous calcium oscillations are obtained for the parameter K a~2 :3 mM {1 , whereas a stationary calcium concentration is found for K a~2 :0 mM {1 in the absence of mechanical feedback. The effect of the parameter K a on the autonomous calcium oscillator model for Physarum was studied systematically in Refs. [31,41]. For the diffusion coefficient of the free calcium we use a typical value of D c~0 :03 mm 2 =min for small ions in cytoplasm [43].
In Physarum, the percentage of actin in the cytoplasm is estimated to be in the range 15{25% [44]. Since other proteins also contribute to the solid gel phase, we set r 0 gel~0 :25 and r 0 sol~0 :75. The Young modulus E of a Physarum strand was determined to be about 10 kPa [45]. For sponge-like materials, measurements show that the Poisson ratio is very low: n&0 [46] implying G&K in two dimensions. Therefore, we have set the parameters G~K~8:9 kPa.
A typical value for the generated tension in plasmodial strands is T a &20 kPa [47]. According to Eq. (16), the active tension T a relaxes to an equilibrium value of T ? a~F T h. Because h max v0:1(see Ref. [48]), values up to F T~3 50 kPa for F T have been considered. Moreover, we have varied F T in our study to illustrate the influence of the strength of mechanochemical coupling on the pattern dynamics.
For the dynamic sol viscosity g shear sol in Physarum, values in the range of 0:1{0:5 Pas where measured [49]. Alternatively, one can calculate the sol viscosity from velocity profiles measured in Physarum [50]. With this method we obtain a dynamic sol viscosity of around 10 Pas. In the model, the value of the sol viscosity is then set to g shear sol~1 Pa s. In a study of a composite network containing actin filaments and microtubuli [51] the frequency dependent complex dynamic shear modulus G(v)~G'(v)ziG''(v) was measured. For a viscoelastic material described by the Kelvin-Voigt model, one identifies G''(v)~vg shear gel . With a typical frequency v~2pmin {1 of oscillations in Physarum a value of g shear gel &1Pas is obtained. For this value of gel viscosity the timescale, on which viscous gel stress is relevant is of the order of 10 {4 s. Assuming that the relevant timescale here is about 1 min, viscous stress in the gel can be neglected. However, due to a stabilizing effect on the numerics, we keep this term. We neglect the influence from the bulk viscosities and set g bulk sol=gel~0 . Since the permeability k!' 2 pore , the drag coefficient is b!g sol =' 2 pore , where ' pore is the average pore size. Porous structures in Physarum exhibit several spatial scales: The actin network pore size is ' pore~0 :25 mm [34] and the membrane structure of invaginations possesses pores with a typical diameter of about 10 mm [35]. If one assumes that larger porous structures determine the permeability, a drag coefficient of b&10 4 kg=(mm 3 min) is obtained. Since, however, the porous structure of the Physarum cytoskeleton may vary substantially over time, we

Linear stability analysis
Linear stability analysis is used here to investigate the spatiotemporal instability near the homogeneous steady state (HSS) without deformation and without any fluid motion. This HSS is given by u ? :v ? :0, p ?~c onst: and the steady state values of chemical components and active tension that follow from the implicit equations f c (n ? c ,w ? )~0, f w (n ? c ,w ? )~0 and T ? a~F T h(n ? c ,w ? ). The growth rate and the frequency (for imaginary eigenvalues) of a small perturbation (dn c , dw, dT a , du, dv, dp)e iqxzlt to the HSS is given by the dispersion relation l(q).

Numerical integration
The integration of the full PDE model given by Eqs. (13) - (18) was done with a hybrid method consisting of a finite element (FEM) and finite volume discretization scheme (FVM). This part was implemented as a stand-alone software written in C. Two-dimensional meshes, however, were generated with the free software Triangle [52]. With the operator-splitting technique the time integration step was divided into substeps that are carried out with different types of solvers: a fourth order Runge-Kutta method for the nonlinear reaction part, a linear FEM with implicit Euler stepping for the parabolic and elliptic parts and a FVM using the dual mesh of the triangulation (Voronoi diagram) for the advection step. (Operator splitting leads to subproblems of the form L t n c~D Dn c z:::(parabolic PDE), a linear elastic problem (mDz(mzl)== : )u~:::(elliptic PDE) and an advection problem L t n c z= : (n c w)~::: (hyperbolic PDE).) The number of nodes for the mesh discretizing a disc with radius r~1 mm was N~5218. This results in a typical mesh resolution of Dx&0:024 mm, whereas wavelengths obtained in the simulations do not go below l min &0:6 mm. The ratio Dx=l min v1=40 is assumed to provide a sufficient accuracy in order to resolve the profiles of the wave patterns and allow us to qualitatively compare our simulation results with experimental data. The time step size Dt was chosen to be 0:01 min that is about one order of magnitude smaller than the fastest time scale in the reaction kinetics given by Eqs. (16) -(18). The typical simulation length was 100 min.

Linear stability analysis and dispersion relation
Here, the mechanochemical coupling strength F T and the drag coefficient b are varied to reveal their influence on the stability of the HSS and the shape of the dispersion curve l(q~DqD) . We have carried out numerical simulations where other model parameters were varied (results not shown) and found that indeed the drag coefficient b and the coupling strength F T lead to the largest qualitative changes. In contrast, changes in sol viscosity by one order of magnitude had almost no effect on the simulation results. In Fig. 3 the branch of eigenvalues with the largest real part is shown for two scenarios (a and b) with different affinity K a and for different values of F T . The real and imaginary parts of the  dispersion relation are displayed in separate plots. Fig. 3 a) shows the case where the HSS is stable in absence of mechanochemical coupling (F T~0 ). If F T is increased above a critical value, the HSS exhibits a wave instability (oscillatory Turing instability). In Fig. 3 b) we consider a case where the HSS is already unstable against oscillations for F T~0 , i.e. Re(l(0))w0. For larger wavenumbers q of the perturbation the imaginary part of the growth rate increases indicating that the frequency of waves is larger than the frequency of homogeneous oscillations.
In the case shown in Fig. 3 b), there are two mechanisms that destabilize the HSS: the homogeneous oscillatory mode with q~0 that originates from the calcium kinetics described by Eqs. (7) and a wave-like perturbation at finite wavenumber (q=0) induced by the mechanical feedback. Above a critical mechanochemical coupling strength F T wF thr T w0 the fastest growing mode (mode with largest real part of the eigenvalue) has a finite wavelength q=0 and a nonzero imaginary part indicating wave dynamics.
In Fig. 4 a) the influence of variation of the drag coefficient b on the dispersion relation for fixed F T is shown. The wavelength L c~2 p=q c of the perturbation with largest growth rate increases with growing b until the maximum with finite wavelength in the dispersion curve disappears. This discontinuity is visible in Fig. 4 b) where L c is plotted versus the drag coefficient b for different values of F T . This figure also shows that the fastest growing wavelength L c decreases with the mechanochemical coupling strength. Nevertheless, the wavelength depends only weakly on b over a large range of values.
The linear stability analysis gives already an important insight into the essential physics of the model defined by Eqs. (13) -(18): The mechanical part leads to a fastest growing mode with nonzero wavenumber while the calcium oscillator can switch the longwavelength instability observed in the purely mechanical system [30] to a short wavelength instability by removing the conservation constraint for the integrated concentration of free calcium.

Numerical simulations
For different values of the coupling strength F T and the drag coefficient b, we present numerical simulations of the dynamics of full nonlinear model equations. The corresponding phase diagram is shown in Fig. 5. Therein, the solid black line separates the phase plane according to the shape of the dispersion curve obtained from the linear stability analysis of the HSS: It shows the threshold value F thr T as a function of the parameter b. When F T vF thr T there is no discrete wavenumber q k w0, for which l(q k )wl(q 0~0 ). In the opposite case F §F thr T there exists at least one q k w0 with l(q k )wl(q 0 ). The prediction of homogeneous oscillations for F T vF thr T from linear stability analysis is confirmed by numerical  Tables 1  and 2. A video file displaying the spatiotemporal dynamics including the transient initial phase corresponding to subfigure c) is included in the supporting material (see Movie S1 in File S1). doi:10.1371/journal.pone.0099220.g006  Tables 1  and 2. A video file corresponding to subfigure c) is included in the supporting material (see Movie S2 in File S1). doi:10.1371/journal.pone.0099220.g007 simulations for bv10 5 kg=(mm 3 min) (see Fig. 5, blue discs). However, for larger values of b one finds a considerable discrepancy between linear stability analysis and simulation results. For homogeneous oscillations in the calcium concentration the flow and deformation fields both vanish.
To address the question how much influence the advective coupling relative to diffusion has, one can consider a Péclet number that is given by the ratio of diffusive to advective time scales (see also Ref. [22] for a motivation of this definition) where h max &0:01 is the typical amplitude for the oscillations in the variable h appearing in Eq. (10) for K a~2 :3 mM {1 . In Fig. 5 the blue line corresponds to Pe~1. Note, that above the line Pew1, there are no homogeneous oscillations and different types of patterns prevail. In some cases, simple patterns like rotating spirals are traveling waves are also found for Pev1. Altogether, this consideration shows that the mechanochemical coupling has to be strong enough to overcome the homogenizing effect of diffusion for mechanochemical waves and patterns to emerge. Above F thr T , traveling, standing and spiral waves occur in the simulations. In Fig. 6 a traveling wave is depicted by snapshots and space-time plots. In addition to the free calcium concentration n c the plots show the relative height field h and the sol velocity v, since these quantities are most likely to be measured in experiments. We get an average velocity of the wave front that is about 0:086 mm=s.
For values of b around 5 : 10 4 kg=(mm 3 min) standing waves are observed for the relative height field h, while for the concentration field n c a traveling wave appears (see Fig. 7).
Traveling and standing waves often coexist with single rotating spiral patterns. The pattern selection depends on the initial conditions, in particular on the number of phase singularities. Two counter-rotating spirals annihilate and give way to a traveling wave, whereas a single rotating spiral is stable. A single rotating spiral is presented in Fig. 8. No distinct direction of propagation or rotation is preferred reflecting the symmetries of the system. The coexistence region of traveling waves and spirals is marked by violet squares in the phase diagram in Fig. 5. From the space-time plots in Fig. 8 c) one finds a wave speed of 0:047 mm=s.
For larger mechanochemical coupling strength F T , there is no coexistence of spirals with traveling waves: a rotating wave is the only attractor (see Fig. 5, green triangles). For large enough values of the drag coefficient b a mode with largest possible wavenumber L 1~2 L determines the emerging patterns. This is confirmed by the numerical simulations. For bv10 5 kg=(mm 3 min) periodic patterns are obtained, even for large coupling F T (see Fig. 5,  brown squares). These patterns have a characteristic wavelength of the same order as the system size. The simplest pattern is an antiphase oscillations in the form of a radial wave that is reflected at the boundaries (see Fig. 9) that has the symmetry as the experimentally found pattern in Fig. 1d.
For large drag coefficients (see Fig. 10 and phase diagram in Fig. 5, pink circles), irregular patterns like the experimental pattern in Fig. 1b with a wavelength significantly smaller than the system size are obtained. For the wave segments in these spatiotemporal patterns we get a typical velocity of 0:03 mm=s that is much slower than that of traveling or spiral wave.

Discussion
In this article we have combined a recently published novel mechanical continuum model of the cytoplasm as a poroelastic active medium [30] with the Smith-Saldana model for calcium oscillations in Physarum protoplasm [31]. Upon increase of the mechanochemical coupling strength the homogeneous steady state in this model is destabilized by an oscillatory Turing instability with finite wave number. Homogeneous oscillations of the calcium concentration are replaced by spatiotemporal patterns and waves connected with local deformation and fluid motion. Practically all spatiotemporal contraction patterns observed experimentally in Refs. [18,19] including rotating spirals, traveling and standing waves, antiphase oscillation and irregular, chaotic waves could be reproduced by numerical simulations of this model. In contrast to earlier models (see e.g. Ref. [20]) for Physarum protoplasm, a closed set of mechanical force-balance equations is derived that allows for explicit computations of pressure and flow fields. The cytoplasm is treated as a two-phase material, consisting of a passive fluid sol phase and an active solid viscoelastic gel phase. The basic idea was already sketched in [41], where however only the case was considered, where the mechanochemical coupling is approximated by a global coupling in the dynamics for calcium. Here, we have instead analyzed and simulated the interplay of mechanical deformation of the cytoskeleton, fluid flow in the cytosol and chemical ( = calcium) concentration and the associated spatiotemporal dynamics.
Unlike other models that treat the cytoskeleton as a viscoelastic fluid (see e.g. Ref. [53]), we consider the cytoskeleton to be a viscoelastic solid described the Kelvin-Voigt model. This is a valid approximation in Physarum, since the relaxation time of elastic tension is larger by about a factor of three than the calcium oscillation period in the system [32]. The porosity of the solid gel phase allows a flow of cytosol that carries calcium and regulates the tension generation. This is represented by an advective transport of calcium in addition to diffusion and introduces a nonlinear feedback mechanism that couples flows and deformations to a nonlinear reaction kinetics. Related models of active gels and fluids [22,24,54] consider the transport of motors, whereas we have considered the transport of calcium that acts as a motorregulating species in Physarum.
We have limited our considerations to small displacements u and deformations D=uD. Thus, the equations are formulated up to linear order in the displacement field and its gradient.
A linear stability analysis of a homogeneous steady state with constant calcium concentration and zero deformation and flow has been carried out for the presented model. The resulting dispersion relations reveal that the mechanical feedback provides a new mechanism of an oscillatory Turing instability with nonzero wavenumber if the mechanochemical coupling strength is positive F T w0. In numerical simulations, homogeneous oscillations get destabilized for sufficiently large coupling strength F T . Rotating spirals and traveling waves obtained in the simulations are common patterns found in experiments [18,19] (and U. Strachauer & M.J.B. Hauser, unpublished data). The wave speed we obtain for traveling and spiral waves agree with the findings in [18]. Note, however, that the antiphase and irregular patterns (Figs. 9 and 10) were obtained at large coupling strength F T yielding also large deformations. Thus, these results violate the small deformation assumption of our model and should be used only for qualitative comparisons. For a quantitative treatment of large deformations, an extension of the above model to a nonlinear elasticity formulation [55] becomes necessary.
The mechanism of the experimentally observed transition between different patterns on a timescale much larger than the typical calcium oscillation period is not fully understood. The phase diagram in Fig. 5 shows that a minor change in the coupling strength F T or drag coefficient b can result in a different type of pattern. Hence, a variation of one of these parameters over the time of the experiment may be responsible for qualitative changes.
For the sake of simplicity, we choose a simple Kelvin-Voigt model as a starting point and found that it is sufficient to reproduce the experimental findings on deformation patterns in Physarum. Future work, e. g. on the modeling of moving droplets of Physarum, may require a more complicated model that exhibits In order to model the stage of migration of a Physarum droplet free-boundary conditions have to be imposed and an interaction with the substrate must be considered (see e.g. [36]). Such an extension of the model may eventually help to understand the interplay of chemical and mechanical processes in the selforganized amoeboid movement of Physarum droplets.

Supporting Information
Comments about the sol continuity equation, the osmotic pressure, as well as a derivation of the force balance equations from an extremal principle can be found in the Supporting Information (Text S1 in File S2) provided with this article. Furthermore, we include movies (see Movies S1-S5 in File S1) with the spatiotemporal dynamics corresponding to Figs. 6 -10.

Supporting Information
File S1 Supporting movies. Movie S1. Spatiotemporal dynamics of traveling wave pattern: Relative height field h in color and flow field v given by arrows corresponding to Fig. 6. The simulation time is 100min beginning with the initial state. Movie S2. Spatiotemporal dynamics of standing wave pattern: Relative height field h in color and flow field v given by arrows corresponding to Fig. 7. The simulation time is 100min beginning with the initial state. Movie S3. Spatiotemporal dynamics of spiral wave pattern: Relative height field h in color and flow field v given by arrows corresponding to Fig. 8. The simulation time is 100min beginning with the initial state. Movie S4. Spatiotemporal dynamics of radial wave: Relative height field h in color and flow field v given by arrows corresponding to Fig. 9. The simulation time is 100min beginning with the initial state. Movie S5. Spatiotemporal dynamics of irregular wave pattern: Relative height field h in color and flow field v given by arrows corresponding to Fig. 10. The simulation time is 100min beginning with the initial state. (TAR) File S2 Supporting information. Text S1. This next contains some more detailed information about the derivation of the mechanical model and a comment about the negligence of the osmotic swelling pressure. (PDF)