Oscillatory motion of a droplet in an active poroelastic two-phase model
Dirk Alexander Kulawiak, Jakob Löber, Markus Bär, Harald Engel
OOscillatory motion of a droplet in an activeporoelastic two-phase model
Dirk Alexander Kulawiak , Jakob L¨ober , Markus B¨ar ,and Harald Engel TU Berlin - Institut f¨ur Theoretische Physik, Hardenberstr. 36, 10623Berlin, Germany; Max-Planck-Institut f¨ur Physik komplexer Systeme,N¨othnitzer Straße 38, 01187 Dresden, Germany; Physikalisch-TechnischeBundesanstalt, Abbestrasse 2-12, 10587 Berlin, GermanyE-mail: [email protected]
July 2018
Abstract.
The onset of self-organized droplet motion is studied in aporoelastic two-phase model with free boundaries and substrate friction. In themodel, an active, gel-like phase and a passive, fluid-like phase interpenetrate onsmall length scales. A feedback loop between a chemical regulator, mechanicaldeformations, and induced fluid flow gives rise to oscillatory and irregulardroplet motion accompanied by spatio-temporal contraction patterns inside thedroplet. By numerical simulations in one spatial dimension, we cover extendedparameter regimes of active tension and substrate friction, and reproduceexperimentally observed oscillation periods and amplitudes. In line with recentexperiments, the model predicts alternating forward and backward fluid flowat the boundaries with reversed flow in the center. Our model is a firststep towards a more detailed model of moving microplasmodia of Physarumpolycephalum.
1. Introduction
Dynamic processes in cells, and cell motility in particular, are intriguingexamples of large-scale spatio-temporal order in systems far from thermodynamicequilibrium [1, 2, 3]. Here, the continuous turnover of ATP by molecularmotors [4] provides the energy to drive mechano-chemical contraction-expansionpatterns and, ultimately, locomotion. Biological examples of these phenomenaare reviewed and discussed in [5, 6, 7, 8, 9].A well-studied model organism exhibiting a huge variety of spatio-temporalmechano-chemical patterns with and without locomotion is the true slime-mold
Physarum polycephalum [10, 11, 12]. Physarum is unicellular, but a cell containsmultiple nuclei and can grow to the size of several square meters [13]. Physarummicroplasmodia are an artificial form of Physarum with a size between 100 µ mand 1 mm that do not occur in nature [14, 15]. They are composed of a gel-like a r X i v : . [ c ond - m a t . s o f t ] A ug ctive Poroelastic Two-Phase Model peristaltic and amphistaltic [16, 21]. In both modes, microplasmodia alternate between forwardand backward motion with a well-defined period. The forward motion is largerthan the backward motion, resulting in a net displacement within each period.In the more frequently observed peristaltic mode, motion is driven by mechano-chemical waves originating at the tail and traveling towards the front. In theamphistaltic mode, front and tail contract in anti-phase oscillations.Common models for the cytoskeleton are based on active fluid and gel models[22, 23, 24, 25, 26]. In contrast, some models for the crawling type of amoeboidcell motility [27, 28] neglect intracellular flows. As opposed to simple fluids andsolids, which are governed by a single momentum balance equation, poroelasticmedia belong to the class of two-fluid models. These are characterized byindividual momentum balance equations for each of the constitutive phases. Sucha description is useful if two phases with largely different rheological propertiesinterpenetrate on relatively small length scales, such as groundwater permeatingporous rock [29], the superposition of normal and inviscid superfluid helium inhelium II [30], or cytosol pervading the cytoskeleton.Poroelastic two-phase models have been used successfully [31, 32, 33]as ingredients in detailed models to replicate the pattern found in restingmicroplasmodia of Physarum [14]. Our work is based on the simple generic modelof a poroelastic active droplet introduced by Radszuweit et al. in [34]. Therein, afeedback loop between an advected chemical regulator and active stress gives riseto self-organized spatio-temporal contraction patterns. This occurs even withoutthe inclusion of a nonlinear reaction-diffusion kinetics for the chemical regulator,that were part of the detailed Physarum models mentioned above.Appropriate boundary conditions must be introduced to close theseporoelastic models and all earlier approaches utilized fixed boundaries to studyresting microplasmodia. While these fixed boundaries are simpler to implementin numerical simulations and allow us to study mechanical deformations in thebulk, they preclude the possibility of deformations of the droplet boundary andtherefore motion of the droplet as a whole.Here, we introduce free boundary conditions that allow for deformationsand motion of the droplet boundary. Considering a free boundary problemcomplicates numerical simulations and can significantly change the solution ofa given problem, especially for fluid dynamics [35, 36].Previous work has shown that our model is able to exhibit self-organizedspatially non-symmetric deformations [34]. In order to describe the impact of ctive Poroelastic Two-Phase Model
2. Model
In our poroelastic two-phase model, the homogeneous isotropic droplet consistsof an active, gel-like phase and a passive, fluid-like phase that interpenetrate atrelatively small length scales [37, 38]. The passive phase flows with velocity v.The active gel phase is a visco-elastic solid with mechanical displacements u andvelocity ˙ u .Both phases individually satisfy a momentum balance equation expressedwith stress tensors, where σ g and σ f denote the stress in the gel and in the fluid,respectively. The total stress is given by σ = ρ g σ g + ρ f σ f , with ρ g ( ρ f ) denotingthe volume fraction of the gel (fluid) phase. Deformations of the medium resultin an exchange of volume between the respective fractions. The time evolutionfor the gel fraction is given by ˆ ρ g = ρ g (1 − ∂ x u ). Here, ρ g is the initial, spatiallyconstant gel fraction. However, due to our small strains approximation, only theconstant term ρ g enters in our model equations. Furthermore, we assume thatthere are no other phases present, and the volume fractions satisfy ρ g + ρ f = 1at all times [39].A droplet occupies a one dimensional, time dependent domain B withboundaries denoted by ∂ B . Assuming that B is infinitely large in the y-direction,the boundary is straight, and we omit terms that depend on interface tensionor bending. Free boundary conditions in x-direction enable the boundary todeform and move in response to bulk flow and deformation [40, 41, 42]. Weassume that the droplet is surrounded by an inviscid fluid described with stresstensor σ out = − p out . The exact value of the outside hydrostatic pressure p out ctive Poroelastic Two-Phase Model p out = 0. At the droplet’s boundary, the total stress has to be continuous acrossthe interface. This gives the first boundary condition σ − p (cid:12)(cid:12)(cid:12) ∂ B = σ out (cid:12)(cid:12)(cid:12) ∂ B = 0 , (1)where the subscript ∂ B denotes evaluation at the boundary of domain B . Becauseof the two momentum balance relations in the poroelastic model, we need asecond boundary condition. Assuming there is no polymerization of actin at andno permeation of the fluid phase through the boundary, the velocity of gel andfluid must match. This gives rise to the additional boundary condition˙ u (cid:12)(cid:12)(cid:12) ∂ B = v (cid:12)(cid:12)(cid:12) ∂ B . (2)Note, that the free boundary conditions require the evaluation of stress tensorsand flow fields at the boundary, whose position itself must be determined inthe course of solving the evolution equations. We circumvent this problem bytransforming the system to a co-moving frame of reference. In general, continuummechanics allows us to use different coordinate frames to formulate the modelequations [43].We distinguish between the lab frame (LF) with spatial coordinates X andthe gel’s body reference frame (BRF) with material coordinates x . The materialdisplacement field u connects both frames by u ( x, t ) = X ( x, t ) − x . Note thatthe domain as well as its boundary, which is time-dependent in the LF, becomesstationary in the BRF. The gel velocity ˙ u is given by the material time derivative˙ u = ∂ t u + ( ∂ x u ) ˙ x . By definition, the gel is fixed in its BRF ( ˙ x = 0), and thematerial time derivative simplifies to ˙ u = ∂ t u . On the downside, transformingstress tensors given by linear constitutive laws from the LF to the BRF gives riseto many geometric nonlinearities. We simplify by linearizing in the strains, i.e.,assuming | ∂ x u | (cid:28)
1. However, note that we do not assume the displacements u to be small. The displacements may, for example, grow linearly in time withoutbounds for a droplet moving with constant center of mass velocity. See [34, 31]for details on the transformation from the LF to the BRF and Fig. 6 for a visualcomparison of a quantity plotted in the BRF and the LF.The fluid phase is modeled as a passive viscous liquid with stress σ f = η f ∂ x v ,where η f is the viscosity. The stress σ g = σ ve + σ act of the gel phase is decomposedin a passive part, σ ve , and an active part, σ act . We assume that the passive part isa viscoelastic Kelvin-Voigt solid with σ ve = E∂ x u + η g ∂ x ˙ u [44]. Here, E is Young’smodulus and η g the viscosity. Recently, the effects of alternative viscoelasticmodels for the passive gel stress were investigated in [45]. A poroelastic modelwith nonlinear elasticity was introduced in [46].Hence, we can write momentum balances in the BRF for both phases. TheReynolds numbers that arise from flows in the medium are assumed to be smalland Re (cid:28)
1. Thus, inertia effects can be neglected, and the intra-droplet flow is ctive Poroelastic Two-Phase Model ρ g ∂ x ( σ g − p ) + f g + f fric = 0 , (3) ρ f ∂ x ( σ f − p ) + f f = 0 , (4)where p is the pressure stemming from the incompressibility of the mediumexpressed as ∂ x ( ρ g ˙ u + ρ f v ) = 0 . (5)The friction between both phases is given by Darcy’s law together with Newton’sthird law f g = − f f = ρ g ρ f β ( v − ˙ u ). We assume a linear substrate friction force f fric = − ρ g γ ˙ u for friction between gel and substrate, and no friction between fluidand substrate.The active stress is assumed to be governed by the concentration c of achemical regulator species σ act = T ( c ) = T − ξ c c . (6)Here, T is a homogeneous stress that is inhibited by the regulator c , and ξ describes the strength of this active stress. This dependency is in contrast withthe assumption of an activating regulator species in the simple model of a one-component active fluid [25] and in line with observation regarding the effect ofcalcium in Physarum microplasmodia [47].The regulator c is dissolved in the fluid and advected with the fluid flow v in the LF. Furthermore, the regulator is diffusing with a diffusion coefficient D c . Transforming the advection-diffusion equation from the LF to the BRF, andlinearizing in the gel strains ∂ x u , yields an advection-diffusion equation with therelative velocity of the fluid to the gel v − ˙ u as the advection velocity, ∂ t c + ∂ x [( v − ˙ u ) c ] = D c ∂ xx c. (7)We assume that no regulator molecules can cross the droplet’s membrane,resulting in a no-flux boundary condition for the variable c . Thus, the totalamount of the regulator is conserved. Nevertheless, there can be local differencesin the regulator concentration, which in turn drive mechanical deformations via aspatially varying active stress T ( c ). When necessary, a nonlinear reaction kineticsfor the regulator species can be taken into account, for details see [31, 33].In summary, the model equations are given by ρ f η f ∂ xx v + ρ g η g ∂ xx ˙ u + ρ g E∂ xx u − ρ g γ ˙ u − ∂ x p = − ρ g ∂ x T ( c ) (8) η f ∂ xx v − ρ g β ( v − ˙ u ) − ∂ x p = 0 (9) ∂ x ( ρ g ˙ u + ρ f v ) = 0 (10) ∂ t c + ∂ x [( v − ˙ u ) c ] − D c ∂ xx c = 0 . (11)Note, that all patterns in our model emerge through self-organization.Furthermore, the only nonlinear terms are the advection term for the regulator ctive Poroelastic Two-Phase Model ∂ x T ( c ). We introduce the dimensionless P´ecletnumber Pe = ξ/ ( D c β ) as a measure for the ratio of diffusive to advective timescales to characterize the strength of the active tension [34].
3. Results
Equations 8-11 are solved with parameter values adopted from [34] and listedin Tab. 1 unless stated otherwise. The initial condition is the weakly perturbedhomogeneous steady state (HSS) with u = ˙ u = v = 0 and c = c = 1. Due to theincompressibility of the medium (Eq . L is constant, andthe location of the droplet’s left boundary is used as a measure of its position.We observe fluid flow coupled to concentration and deformation patters both forresting and moving boundaries. Depending on parameter values of the frictioncoefficient γ and the P´eclet number Pe, the droplet’s position over time remainsfixed or undergoes regular or irregular oscillations. However, in all cases, thetemporal average of the droplet’s position vanishes, i.e. there is no net motion.To visualize the dynamics of regulator and gel flow, we show space-time plots ofthese quantities in the BRF. Time in s − − P o s i t i o n i n µ m Time in s x i n µ m Regulator c
Figure 1. Regular oscillations of the droplet’s position overtime (top) are accompanied by spatially antisymmetric regulatoroscillations (bottom).
The active tension (Pe = 6) is slightly above thecritical value where the HSS destabilizes. After initialization, the restingdroplet exhibits spatially symmetric regulator oscillations with period T s = 24 sand slowly growing amplitude. After ≈ T as = 98 s. The regulator concentration is plotted in a body referenceframe co-moving with the gel phase. Symmetric and antisymmetric spatio-temporal oscillations
The case of a regular oscillation in Fig. 1 exemplifies how the different parts of themodel interact and motion arises. Shortly after initialization, the concentration ctive Poroelastic Two-Phase Model Symmetric oscillations, t = 0 v drop Antisymmetric oscillations, t = 0 x in µm Symmetric oscillations, t = 0.5 T s x in µm v drop Antisymmetric oscillations, t = 0.5 T as -5.00.05.0 0.02.04.06.0-5.00.05.0 -6.0-4.0-2.00.0 R e g u l a t o r c A c t i v e t e n s i o n ∂ x T i n − kg µ m s Figure 2. Snapshots of spatially symmetric (left) and antisymmetric(right) regulator oscillations.
The regulator concentration c is depicted inblue, the active tension ∂ x T ∼ ∂ x c/ (1 + c ) in red, and black arrows indicatethe advection velocity v − ˙ u . The length of the black arrows indicates theamplitude of the advective flow. Spatially symmetric oscillations have a periodof T s = 24 s and c oscillates between high concentration at the center (top left)and high concentration at the boundaries (bottom left). No movement occurs.In the second case, c oscillates between two configurations antisymmetric toeach other with a period of T as = 98 s (right). The droplet is moving withvelocity v drop towards the direction of high regulator concentration as indicatedby the green arrows. oscillates in a spatially symmetric manner with period T s = 24 s (Fig. 1, bottom)and a constant position of the droplet (Fig. 1, top). The blue lines in Fig. 2 showsnapshots of the regulator concentration as a function of space. The regulatordistribution changes periodically from a high concentration at the center and lowvalues at the boundaries (top left panel) to a low concentration at the centerand high values at the boundaries (bottom left panel). In the top left panel, theactive tension ∂ x T ∼ ∂ x c/ (1+ c ) (red line, Eq .
6) generates a symmetric advectionflow from the boundaries towards the center (black arrows). This results in evenmore regulator piling up at the center, thus causing an even stronger flow. Theflow deforms the droplet and elastic tension in the gel builds up. Over time, thiselastic tension increases and counteracts to the active tension. The advection ofregulator weakens, and at some point regulator diffusion takes over. The regulatorconcentration at the center starts to decrease, diminishing active tension andthereby advection even further. When the elastic tension overcomes the activetension the direction of advection changes. The concentration starts to pile up atthe boundaries (bottom left panel), and the process repeats. During the wholeoscillation cycle both the profiles of gel ˙ u and fluid v are symmetric in space.Hence, spatially symmetric oscillations result in immobile droplets.At t ≈ ctive Poroelastic Two-Phase Model ≈
200 s, the amplitude of regulator oscillations saturates,and the droplet performs periodic motion in phase with the regulator oscillationsand a period of T as = 98 s. Fig. 2 (right) shows two snapshots of antisymmetricoscillations at different times. When the droplet’s position is at its maximumdisplacement to the right, the regulator concentration changes from having itsmaximum at the right boundary to being higher on the left boundary. When themaximum of the concentration switches sides, the active tension changes sign,causing an advective flow to the left. A rapid movement to the left starts witha maximum droplet velocity of v drop = 2 . µ m / s. While this movement takesplace, the concentration piles up on the left boundary (top right) and the elastictension inside the droplet builds up. Similar to the symmetric case, this elastictension is acting in opposition to the active tension yielding a lower advection anddroplet velocity. At some point, advection becomes weaker than diffusion. Then,the concentration has reached its maximum and begins to decrease. While theactive tension is diminishing, the elastic stress inside the gel phase still buildsup. Once it overcomes the active tension, the tension in the droplet startsto decay. This causes a change of direction of the advective flow, leading toa rise of concentration on the right, and the process repeats (bottom right).This antisymmetric oscillation results in a periodic movement with a maximumdisplacement of u max ≈ µ m but not in net motion. Varying active tension strength and substrate friction excites different modes ofmotion
Depending on the strength of the active tension as measured by the dimensionlessP´eclet number Pe and the friction between droplet and substrate, the droplet’sposition is stationary or oscillates with one or more frequencies. We decomposethe position over time data in frequency components and use the number ofexcited Fourier modes K to characterize a droplet’s movement. We group theresults into four different categories: “HSS” if the HSS is stable against smallperturbations, “non-moving”, when the HSS is unstable, but the droplet doesnot move ( K = 0), “regular” in the case of movement with less than 5 modes(1 ≤ K ≤
4) and “irregular” if K ≥
5. Fig. 3 gives an overview which modeof movement emerges (top) and how the droplet speed develops (bottom) whenchanging active tension strength Pe and the substrate friction γ . ctive Poroelastic Two-Phase Model − − − − − γ - Substrate friction in kg/s P e - A c t i v e t e n s i o n HSS Non moving Regular Irregular − − − − γ - Substrate friction in kg/s M e a n s p ee d i n µ m / s Pe = = = Figure 3. Mode of motion (top) and speed (bottom) for differentstrengths of substrate friction γ and active tension as measured bythe dimensionless P´eclet number Pe = ξ/ ( D c β ). Top: We group theresults into four different categories: “HSS” if the HSS is stable against smallperturbations, “non-moving”, when the HSS is unstable, but the droplet doesnot move ( K = 0), “regular” in the case of periodic movement with less than5 Fourier modes (1 ≤ K ≤
4) and “irregular” if K ≥
5. The thick black linedenotes the critical P´eclet number Pe cr that results from the linear stabilityanalysis. Bottom: The mean droplet speed increases for larger P´eclet numbersand shows a maximum for medium values of γ . In agreement with the linear stability analysis from [34], the HSS is stableagainst small perturbations if Pe is below a critical P´eclet number Pe cr . Abovethis critical value, the HSS destabilizes and non-linear numerical simulations haveto be carried out to determine the resulting patterns. The different categoriesof movement often arise in bands that appear or disappear when the parameterschange.The linear stability analysis predicts an oscillatory instability to long-wavelength modes where initially a symmetric standing wave-like pattern withtwo nodes appears from a superposition of left and right traveling waves. In thelong term, however, an antisymmetric standing wave pattern with one node ofdouble wavelength emerges. An adaptation (coarsening) of the pattern towardslarger wavelength is often found in such systems. Overall this coarsening is acrucial point, because we show that the change in symmetry of the pattern iswhat leads to notable motion of the boundary.If γ is too large, the friction inhibits any pattern formation. The HSS isstable, and small perturbations decay. A higher substrate friction shifts Pe cr tohigher values, and for γ = 10 − kg / s the HSS is always stable. If γ is too low( < − kg / s), there is effectively no friction between droplet and substrate. Inthis case, our equations do not possess full rank any more and we omit showingthese results.In the region in-between, we can observe regular and irregular oscillatorymotion as well as non-moving droplets. For a high substrate friction, thedroplet performs regular oscillations. Beginning with γ = 5 × − kg / s, irregular ctive Poroelastic Two-Phase Model ≥ γ this band of irregular solutions isshifted to a lower value of Pe. With γ = 10 − kg / s non-moving solutions occurfor intermediate Pe. In contrast to the case of transient symmetric oscillationsin Fig. 1, here the symmetric oscillations remain stable on the timescale of thesimulation length (at least 10 000 s).For γ = 10 − kg / s the critical P´eclet number is Pe cr = 5 .
5. For Pe > . ≥ .
5, the droplet moves again and the number ofexcited Fourier modes K rises fast. With a further increase of Pe the droplet’smotion becomes irregular. Between Pe = 8 . . . γ . For Pe = 6 . γ ≤ − kg / s and γ ≥ × − kg / s).In-between it performs regular oscillations with a mean speed of 2 µ m / s. Forhigher values of Pe, the droplet’s speed is shifted to higher values with a maximumspeed of about 3 . µ m / s for Pe = 8 . . µ m / s for Pe = 12 .
5. Just asfor a lower active tension, droplets are at rest for a strong friction and the speedapproaches an almost constant value for a weaker friction.
Time in s − − P o s i t i o n i n µ m
300 350 400 450
Time in s x i n µ m Regulator c Figure 4. Irregular oscillations of the droplet’s position over time(top) appear together with asymmetric regulator profiles (bottom).
With an active tension (Pe = 9 .
5) higher than in Fig. 1 the droplet’s movementand the regulator dynamics become irregular and the spectrum of its trajectoryis continuous (data not shown). ctive Poroelastic Two-Phase Model Comparison of our model with experiments on Physarum microplasmodia
In the following, we compare the predictions of our model with some experimentalresults recently reported by [16, 21]. In experiments, directed motion of Physarummicroplasmodia was found to be accompanied by oscillations with a periodbetween 85 s and 110 s. During the peristaltic mode of motion, microplasmodiaundergo a forward displacement of d F ≈ µ m, followed by a backwarddisplacement of d B ≈ µ m-20 µ m. As displayed in Fig. 3, we find regularoscillations for large parameter regimes in our simulations. As shown in Fig. 1,once the dominant pattern with regular, antisymmetric regulator oscillations hasemerged, the droplet’s position over time oscillates with a period of T = 98 s,which in line with the measurements, and undergoes displacements of d ≈ µ m.However, we always find d = d F = d B such that no net motion occurs. x i n µ m Gel velocity ˙ u in µm/s − Figure 5. Gel velocity in simulation (top) and experiment (bottom).
At a fixed position, the flow direction alternates between forward and backwardwith a period of about 100 s in experiment as well as in simulation. Flow atthe center is opposite and of weaker magnitude than flow at the boundaries.Experimental backward flow is of weaker magnitude than the forward flow,resulting in a a net motion. Simulated forward and backward flows have equalmagnitude and thus cancel exactly. Bottom figure taken from [16]. Copyright:Journal of Physics D: Applied Physics by IOP Publishing.
Additionally, we compare experimentally observed flow patterns within themicroplasmodia with our simulations. The space-time plot in Fig. 5 shows thegel flow (top) and experimentally obtained ectoplasmic flow (bottom) for theperistaltic mode of motion from [16]. At a fixed position in our simulations, theflow alternates periodically between forward and backward flow. While the flowdirections at front and back are equal, with a higher magnitude at the back,the center part is flowing towards the opposite direction with weaker magnitudethan flow at the boundaries. Experimentally, flows towards the front have alarger magnitude than backward flows, resulting in a net propagation velocityof the entire microplasmodia of v exp ≈ . . µ m / s. In simulations, forwardand backward flows have equal magnitude and thus cancel within one oscillationperiod.Spatio-temporal measurements of the regulator dynamics ( Ca ) reveal thatmicroplasmodia motility is accompanied by calcium waves. These resembletraveling waves in the peristaltic mode and standing waves in the amphistaltic ctive Poroelastic Two-Phase Model
4. Discussion
To model motile cells, continuum mechanical models must be supplemented withfree boundary conditions. Here, we extended the poroelastic model with rigidboundaries from [31, 32, 33, 34] to the case of free boundaries and included linearfriction with the substrate. In this minimal model, we explore the conditions forself-organized motion of an active poroelastic droplet. This model is a first steptowards a more detailed description of moving Physarum microplasmodia.We observed different modes of motion ranging from resting droplets todroplets performing regular and irregular oscillations. The symmetry breakingfrom a standing wave with mirror symmetry and wavelength equal to the systemlength to an asymmetric standing wave with half a wavelength in the systemafter a long transient is necessary for motion of the boundary (Fig. 1 and Fig. 2).In [34], only symmetric oscillations were reported, and there were no transitionsto asymmetric oscillations. In addition, we found that the droplet’s speed has apeak for intermediate values of γ (Fig. 3).In our earlier work [34], traveling wave patterns appeared that were reflectedat the boundaries and moved back and forth. The chaotic pattern is anintermediate state between the region of stable standing waves (symmetric orasymmetric) and the traveling domain pattern. The different patterns havemostly physiological significance as the quantitative and qualitative aspects ofdroplet motion in the model change.In the parameter plane spanned by the friction coefficient ( γ ) and thestrength of active tension (Pe), we identified parameters that reproduceexperimentally observed oscillation periods of about 100 s [16, 21]. Additionally,as shown in Fig. 5, the simulated flow patterns are qualitative in line withexperimentally measured flow patterns.The obvious question arises why, for all cases of periodic and evenirregular motion, the time averaged position vanishes, and no net motionoccurs. An essential ingredient to achieve directed motion is a mechanismwhich breaks the front-back symmetry and thus establishes a polarity [49]. Theregulator distributions produced by pure diffusion-advection dynamics may lookasymmetric at certain instants in time, but the long time averaged distributionis always symmetric. While this is expected for a regular oscillation such asFig. 1, it is surprising for the irregular dynamics as shown in Fig. 4. The absenceof a time-averaged asymmetry in the regulator dynamics for the irregular caseindicates that an additional mechanism to establish a polarity is required to model ctive Poroelastic Two-Phase Model c is to include reaction kinetics for the regulator as donein models for resting Physarum droplets earlier [31, 33]. There, unidirectionaltraveling mechano-chemical waves were reported in contrast to the back andforth moving waves found in our study here and previously in [34]. The reactionkinetics allows for a temporal variation of the total concentration of the regulatorencoded in the variable c . Additionally, the total amount of the regulator is notconserved anymore.Moreover, the following argument shows that net motion is in generalimpossible for a substrate friction γ constant in space. For constant mass density,the velocity ¯ v of the droplet’s center of mass is obtained by spatially averagingthe gel velocity¯ v = 1 V (cid:90) B ˙ udx, (12)where V denotes the constant volume of the droplet (length in one spatialdimension). Adding the force balances for gel and fluid phases, Eq . . v = 1 V (cid:90) B γ ∂ x ( σ − p ) dx. (13)For γ constant in space, we may use the Gauss theorem to transform the volumeintegral to a surface integral with normal vector n ,¯ v = 1 γV (cid:73) ∂ B ( σ − p ) ndS. (14)Together with the free boundary condition Eq .
1, we immediately obtain avanishing center of mass velocity ¯ v = 0. This is confirmed by a direct calculationof ¯ v in numerical simulations as given by Eq . ctive Poroelastic Two-Phase Model Appendix A: Numerical Details
We solve the equations of motion on an one-dimensional Chebyshev-Lobatto grid[56] of size L with N points. We utilize no-flux boundary conditions for theconcentrations c and free boundary conditions for the mechanical equations. Weformulate our model in the gel’s body reference frame (BRF). For details aboutthe derivation of the model refer to [31, 33]We split the full equations from from Eq . .
11 into a mechanical and anadvection-diffusion part and solve each part separately. We use pseudo-spectralmethods (Chebyshev) for the discretized spatial derivatives and the Euler methodfor time-stepping.For the mechanical part, we introduce U ≡ ∂ t u = u − u t ∆ t , where t denotesthe current time-step and variables without explicit time dependency are at time t + ∆ t . Then, we arrive at u − U ∆ t = u t (15) ρ f η f ∂ xx v + ρ g η g ∂ xx U + ρ g E∂ xx u − γρ g U − ∂ x p = − ρ g ∂ x T ( c t ) (16) η f ∂ xx v − ρ g β ( v − U ) − ∂ x p = 0 (17) ∂ x ( ρ g U + ρ f v ) = 0 . (18)Then, we solve the advection-diffusion part semi-implicitly with ∂ t c = c − c t ∆ t . Thisapproach yields c + ∆ t∂ x (cid:0) w t c (cid:1) − ∆ tD c ∂ xx c = c t , (19)where w = v − U is the fluid velocity in the gel’s BRF. This yields the linearequation (cid:0) − ∆ tD c ∂ xx + ∆ t∂ x w t (cid:1) c = c t . (20)We solve our equations using python [57] with the iterative gmres solver fromscipy and an ILU preconditioner. Appendix B: Regulator concentration in body reference and labframe
We solve our model equations in the gel’s BRF and the resulting quantities aredefined in this frame. However, we as observers are located in the lab frame (LF).Fig. 6 shows how a regular regulator oscillation (compare with Fig. 1 in the main ctive Poroelastic Two-Phase Model u with X = x + u ( x ), where x theposition in the BRF and X is the position in the LF. x i n µ m Body reference frame
Time in s − − X i n µ m Lab frame R e g u l a t o r c Figure 6. Spatially antisymmetric regulator oscillation in bodyreference (top) and lab frame (LF) (bottom).
In the LF, the droplet ismoving into the direction where c has a local maximum. The regular regulatoroscillation yields a periodic movement with a period of about 98 s. The concentration in Fig. 6 exhibits spatially antisymmetric oscillations inthe BRF and switches between a state with a high value of regulator at the leftboundary and a low value at the right boundary and the reversed state. In the LF,this yields a periodic movement of the droplet. As long as the concentration hasa local maximum at a certain boundary the droplet is moving into this direction. ctive Poroelastic Two-Phase Model Appendix C: Parameters
Table 1. Model parameters
Par Description Value Units N Number of grid points 120 -∆ t Numerical time step 0 .
001 s D c Regulator diffusion 200 µ m s − L Length 125 µ m ρ g Gel fraction 0 . ρ f Fluid fraction 0 . η g Viscosity gel 10 − µ m s η f Viscosity fluid 2 × − µ m s β Friction between bothphases 10 − µ m s E Young modulus 0 . kg µ m s Pe Active tension 6 - γ Substrate friction 10 − kg s − These parameters are used throughout this work any deviation is explicitlymarked. These parameters are based on typical estimates for eukaryotic cells.Taken from [34].
Acknowledgments
We thank Markus Radszuweit for helpful discussions about his previous workon the model. Computational resources were provided by the Institut f¨urTheoretische Physik at the TU Berlin. DAK was funded by the German ScienceFoundation (DFG) within the GRK 1558.
5. References [1] Murray J D 2007
Mathematical Biology 1: An Introduction [2] Winfree A T 2001
The Geometry of Biological Time (New York, NY: SpringerNew York) ISBN 978-0-387-98992-1 URL http://link.springer.com/10.1007/978-1-4757-3484-3 [3] Prigogine I 1978
Science [4] Kolomeisky A B and Fisher M E 2007
Annual Review of Physical Chemistry [5] Gross P, Kumar K V and Grill S W 2017 Annual Review of Biophysics [6] Nishikawa M, Naganathan S R, J¨ulicher F and Grill S W 2017 eLife URL https://elifesciences.org/articles/19595 ctive Poroelastic Two-Phase Model [7] Howard J, Grill S W and Bois J S 2011 Nature Reviews Molecular Cell Biology [8] Mayer M, Depken M, Bois J S, J¨ulicher F and Grill S W 2010 Nature [9] Kumar K V, Bois J S, J¨ulicher F and Grill S W 2014
Physical Review Letters
URL https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.112.208101 [10] Oettmeier C, Brix K and Doebereiner H G 2017
Journal of Physics D: Applied Physics
URL http://iopscience.iop.org/article/10.1088/1361-6463/aa8699 [11] Aldrich H C and Daniel J W (eds) 1982
Cell Biology of Physarum and Didymium (NewYork: Academic Press) ISBN 978-0-12-049601-3 URL [12] Teplov V A 2017
Journal of Physics D: Applied Physics ISSN 0022-3727, 1361-6463URL http://iopscience.iop.org/article/10.1088/1361-6463/aa6727/meta [13] Fessel A, Oettmeier C and D¨obereiner H G 2015
Nano Communication Networks [14] Takagi S and Ueda T 2008 Physica D: Nonlinear Phenomena [15] Takagi S and Ueda T 2010
Physica D: Nonlinear Phenomena [16] Zhang S, Guy R, Lasheras J C and del ´Alamo J C 2017
Journal of Physics D:Applied Physics
URL http://iopscience.iop.org/article/10.1088/1361-6463/aa68be/meta [17] Oettmeier C, Lee J and D¨obereiner H G 2018
Journal of Physics D: Applied Physics http://stacks.iop.org/0022-3727/51/i=13/a=134006?key=crossref.6d5c80eecb2c462fdc81a5ff46884ab5 [18] Rieu J P, Delanoe-Ayari H, Takagi S, Tanaka Y and Nakagaki T 2015 Journal of The RoyalSociety Interface URL http://rsif.royalsocietypublishing.org/cgi/doi/10.1098/rsif.2015.0099 [19] Rodiek B and Hauser M J B 2015
The European Physical Journal Special Topics http://link.springer.com/10.1140/epjst/e2015-02455-2 [20] Rodiek B, Takagi S, Ueda T and Hauser M J B 2015
European Biophysics Journal http://link.springer.com/10.1007/s00249-015-1028-7 [21] Lewis O L, Zhang S, Guy R D and del Alamo J C 2015 Journal of The RoyalSociety Interface URL http://rsif.royalsocietypublishing.org/cgi/doi/10.1098/rsif.2014.1359 [22] J¨ulicher F, Kruse K, Prost J F and Joanny J 2007
Physics Reports [23] Joanny J F and Prost J 2009
HFSP Journal [24] K¨opf M and Pismen L 2013 Physica D: Nonlinear Phenomena http://linkinghub.elsevier.com/retrieve/pii/S0167278913001577 [25] Bois J S, J¨ulicher F and Grill S W 2011
Physical Review Letters
URL http://link.aps.org/doi/10.1103/PhysRevLett.106.028103 [26] Strychalski W, Copos C A, Lewis O L and Guy R D 2015
Journal of Computa-tional Physics http://linkinghub.elsevier.com/retrieve/pii/S0021999114006780 [27] Kulawiak D A, Camley B A and Rappel W J 2016
PLOS Computational Biology URL http://dx.plos.org/10.1371/journal.pcbi.1005239 [28] L¨ober J, Ziebert F and Aranson I S 2014
Soft Matter http://xlink.rsc.org/?DOI=C3SM51597D [29] Coussy O 2004 Poromechanics ctive Poroelastic Two-Phase Model NJ: Wiley) ISBN 978-0-470-84920-0 URL [30] Landau L D and Lifschitz E M 1987
Fluid Mechanics [31] Radszuweit M, Engel H and B¨ar M 2014
PLoS ONE URL http://dx.plos.org/10.1371/journal.pone.0099220 [32] Radszuweit M, Engel H and B¨ar M 2010
The European Physical Journal SpecialTopics [33] Alonso S, Strachauer U, Radszuweit M, B¨ar M and Hauser M J 2016
Physica D: NonlinearPhenomena [34] Radszuweit M, Alonso S, Engel H and B¨ar M 2013
Physical Review Letters
URL http://link.aps.org/doi/10.1103/PhysRevLett.110.138102 [35] Balossino R, Pennati G, Migliavacca F, Formaggia L and Dubini G 2006 13[36] Friedman A 2015
Philosophical Transactions of the Royal Society A: Mathematical,Physical and Engineering Sciences http://rsta.royalsocietypublishing.org/lookup/doi/10.1098/rsta.2014.0368 [37] Alt W and Dembo M 1999
Mathematical Biosciences [38] Joanny J F, J¨ulicher F, Kruse K and Prost J F 2007
New Journal of Physics http://iopscience.iop.org/article/10.1088/1367-2630/9/11/422/meta [39] Dembo M and Harlow F 1986 Biophysical Journal [40] Larripa K and Mogilner A 2006 Physica A: Statistical Mechanics and its Appli-cations [41] Nickaeen M, Novak I L, Pulford S, Rumack A, Brandon J, Slepchenko B M and Mogilner A2017
PLOS Computational Biology URL http://dx.plos.org/10.1371/journal.pcbi.1005862 [42] K¨opf M H and Pismen L M 2013
Soft Matter http://pubs.rsc.org/-/content/articlelanding/2013/sm/c3sm26955h [43] Chaves E W V 2013 Notes on Continuum Mechanics [44] Banks H T, Hu S and Kenz Z R 2011
Advances in Applied Mathematics and Mechanics http://journals.cambridge.org/abstract_S207007330000151X [45] Alonso S, Radszuweit M, Engel H and B¨ar M 2017 Journal of Physics D: Applied Physics URL http://iopscience.iop.org/article/10.1088/1361-6463/aa8a1d [46] Taber L, Shi Y, Yang L and Bayly P 2011
Journal of mechanics of materials and structures https://msp.org/jomms/2011/6-1/p35.xhtml [47] Yoshimoto Y, Matsumura F and Kamiya N 1981 Cytoskeleton http://onlinelibrary.wiley.com/doi/10.1002/cm.970010404/full [48] Dupont G, Falcke M, Kirk V and Sneyd J 2016 Models of Calcium Signalling (Cham: Springer International Publishing) ISBN 978-3-319-29645-6 URL http://link.springer.com/10.1007/978-3-319-29647-0 [49] Rappel W J and Edelstein-Keshet L 2017
Current Opinion in Systems Biology [50] Lewis O L and Guy R D 2017
Journal of Physics D: Applied Physics URL http://iopscience.iop.org/article/10.1088/1361-6463/aa76c3/meta [51] Weber C, Rycroft C H and Mahadevan L 2017 arXiv preprint arXiv:1710.03633
URL https://arxiv.org/abs/1710.03633 ctive Poroelastic Two-Phase Model [52] Lu X, Ren L, Gao Q, Zhao Y, Wang S, Yang J and Epstein I R 2013 ChemicalCommunications http://xlink.rsc.org/?DOI=c3cc44480e [53] Ren L, She W, Gao Q, Pan C, Ji C and Epstein I R 2016 Angewandte Chemie InternationalEdition http://doi.wiley.com/10.1002/anie.201608367 [54] Ren L, Wang M, Pan C, Gao Q, Liu Y and Epstein I R 2017 Proceedings of the NationalAcademy of Sciences [55] Epstein I R and Gao Q 2017
Chemistry - A European Journal http://doi.wiley.com/10.1002/chem.201700725 [56] Trefethen L N 2000 Spectral Methods in MATLAB ( Software, environments, tools no 10)(Philadelphia, Pa: SIAM) ISBN 978-0-89871-465-4[57] Oliphant T E 2007
Computing in Science & Engineering10–20 ISSN 1521-9615 URL