Peristalsis by pulses of activity

Peristalsis is a series of actively generated wave-like muscle contractions and relaxations which propagate along the body of an organism. Smooth muscle tissues develop such contractions to produce a peristaltic wave in the digestive tract Sinnott et al. (2017); Brandstaeter et al. (2019). Crawling by peristalsis enables animals like snails, earthworms, slugs, and terrestrial planarians, to advance in narrow spaces  Gray and Lissmann (1938); Chapman (1958); Jones (1978); Alexander (1982); Quillin (1999); Tyrakowski et al. (2012), moreover, based on geometrical symmetries only, peristaltic waves were shown to be an optimal motility strategy in such systems Agostinelli et al. (2018).

In this Letter we develop a prototypical model of a peristalsis in a segmented limbless organism. We assume that it crawls along a flat surface by extending its forward end and then bringing up its rear end. To achieve this goal the organism generates a solitary wave which travels from the front to the rear. The space-time distribution of activity in such living systems is known to be highly adaptive  Boyle et al. (2012a) and the mechanism of this adaptability has recently become a subject of great interest in robotics Daltorio et al. (2013); Boyle et al. (2012b).

Peristaltic waves are also of general physical interest as elementary nonlinear excitations of active matter. Propagating active pulses reminiscent of peristaltic waves are ubiquitous in nature from shimmering in honeybees Kastberger et al. (2014) to Mexican waves on stadiums Cartwright (2006). Comparable phenomena in the form of propagating activity bands are also observed in flocking colonies of swarming robots and other similar systems Solon et al. (2015); Ngamsaad and Suantai (2018). Some of these behaviors can be quantified using models of excitable media Dudchenko and Guria (2012) or models involving some kind of globally synchronized CPG (central pattern generator) Ijspeert (2008). However, such models have been questioned in other cases clearly dominated by mechanical sensory feedback and neuromechanical proprioception Paoletti and Mahadevan (2014). Given the distinctly mechanical functionality of peristalsis, we forgo the reaction-diffusion framework Miller and Dunkel (2020) and neglect the possible role of CPG, and assume instead that physical forces not only drive the associated localized waves of activity but also secure the signaling pathways controlling, for instance, the necessary internal delays.

As a toy model, capturing only the main effects, we consider a mass-spring chain capable of generating active stresses. It is implied that behind such activity is an endogenous machinery of the type involved in muscle tetanization and we assume that the associated energy flow through the system can be adequately represented by a non-constitutive component of stress. We show that the ensuing, apparently purely mechanical, model can generate directional peristaltic locomotion without relying on externally coordinated actuators or digital controllers. Our intentionally minimalistic approach emulates (and can be extended towards) more comprehensive continuum theories of active media with internally generated active stresses known for both fluids  Prost et al. (2015); Jülicher et al. (2018) and solids Hawkins and Liverpool (2014); Maitra and Ramaswamy (2019); Moshe et al. (2018); Scheibner et al. (2019). While these models directly account for energy consumption and energy dissipation, our approximate model neglects both.

Passive solitary waves have been long employed in the actuator-driven soft robotics imitating peristalsis Raney et al. (2016); Nadkarni et al. (2014); Deng et al. (2020). Instead, here we rely on self-driven active solitary waves and show that peristalsis can be modeled by a train of such waves propagating through otherwise passive matter. Rather remarkably, our analysis reveals the existence of a critical motility regime in such systems where active pulses assume realistic rectangular shape with variable width. This ensures broad repertoire of responses and we argue that so-interpreted criticality may be a characteristic feature of the physiological peristalsis.

Bodies of annelid animals are usually divided into a series of metameres, the segments that are fundamentally similar in muscular structure and functionality Tanaka et al. (2012); Quillin (1999). To model such organisms we first neglect friction Kandhari et al. (2019); Keller and Falkovitz (1983); Wadepuhl and Beyn (1989), and represent them schematically as a chain of springs connected in series. The dynamics of such system is described by the FPU equations Keller and Falkovitz (1983)

where $\varepsilon_{n}(t)=(u_{n+1}(t)-u_{n}(t))/a$ is the strain in a spring whose ends undergo displacements $u_{n}(t)$and $a$ is the equilibrium length of a spring. The inertial term, allowing the system to overcome the discreteness-induced trapping, proved to be important in ultra-soft robotics Fang et al. (2015); Jiang and Xu (2017); Tanaka et al. (2012); Deng et al. (2020). In physiological setting the apparent mass density $\rho$ can be viewed as a parameter introducing an activity-related time delay in the response of stretch receptors Badoual et al. (2002); Serra-Picamal et al. (2012).

We assume that the constitutive relation for the stress $\sigma$ has two branches: passive and active, see Fig.1(a). For simplicity, the soft elastic response along the passive branch is considered linear $\sigma=E\varepsilon,$ with elastic modulus $E$. To describe the active branch (analog of muscle tetanus), we write $\sigma=\sigma_{a}+E\varepsilon,$ where $\sigma_{a}>0$ is a constant active stress; in more detailed models it can be a variable with sigmoidal response taking a value zero in the passive phase Paoletti and Mahadevan (2014). We further assume that switching from passive to active response takes place when the ’unjamming’ threshold strain $\varepsilon_{c}$ is reached Serra-Picamal et al. (2012). We neglect the possibility of hysteresis and assume that unloading below the threshold $\varepsilon_{c}$ brings the system back into the passive state, see Fig.1(a). To nondimensionalize the system (1) we normalize length by the system size $L$, time by $L/c$, where $c=\sqrt{E/\rho}$, and stress by $E$.

It will be convenient to first deal with a continuum approximation of the discrete problem. To this end we introduce the continuous strain field $\varepsilon(x,t)$, where $\varepsilon(nh,t)=\varepsilon_{n}(t)$ and assume that $h=a/L\ll 1$; we will also use the convenient rescaling $\tilde{u}=u/\varepsilon_{c},\,\tilde{\varepsilon}=\varepsilon/\varepsilon_{c}$ and $\tilde{\sigma}=\sigma/\varepsilon_{c}$. If we now Padé-approximate the nonlocal operator in the right hand side of (1) and leave only the lowest order terms, we obtain Rosenau (1986); Truskinovsky and Vainchtein (2006); Gómez et al. (2012)

To generate a solitary wave solution of (2) we impose a traveling wave anzats $\varepsilon(x,t)=\varepsilon(\eta)$ where $\eta=(x-Vt)/h$ and $V$ is the dimensionless velocity of the pulse.

If we center the active pulse performing local extension at $\eta=0$ and denote its width by $2d$ we can write the associated stress distribution in the form $\tilde{\sigma}(\eta)=\tilde{\varepsilon}(\eta)+\tilde{\sigma}_{a}\text{rect}(\eta/(2d)),$ where $\text{rect}(x)=H(x+1/2)-H(x-1/2)$ and $H(x)$ is the Heaviside function. We require that $\tilde{\varepsilon}(\eta)\to 0$ as $\eta\to\pm\infty$ and impose at $\eta=\pm d$ the matching conditions $\llbracket\tilde{\varepsilon}\rrbracket=\llbracket d\tilde{\varepsilon}/d\eta\rrbracket=0$ and set $\tilde{\varepsilon}(\pm d)=1$. Then, integrating (2) twice and applying the boundary/matching conditions, we obtain the equation

where the right hand side implicitly depends on $\tilde{\varepsilon}$. Similar solitary wave solution, describing local contraction , can be obtained if we set $\tilde{\sigma}(\eta)=\tilde{\varepsilon}(\eta)+\tilde{\sigma}_{a}(1-\text{rect}(\eta/(2d)))$ and require that $\tilde{\varepsilon}\to\lambda$ when $\eta\to\pm\infty$ where $\lambda=\tilde{\sigma}_{a}/(V^{2}-1)$. The extension pulses exist in the range $1<V<V_{*}$ where $V_{*}\equiv\sqrt{\tilde{\sigma}_{a}/2+1}$, so they are supersonic which does not mean that they are fast given that the underlying elastic medium is almost an acoustic vacuum Nesterenko (2013). The contraction pulses exist in the range $V_{*}<V<V_{**}$ where $V_{**}\equiv\sqrt{\tilde{\sigma}_{a}+1}>V_{*}$.

The phase portrait of the system (3) at $1<V<V_{*}$ is shown in Fig. 1(b). Two non-degenerate saddle points at $\tilde{\varepsilon}_{\pm}$ lie on the same Rayleigh line $V^{2}=(\tilde{\sigma}(\tilde{\varepsilon})-\tilde{\sigma}(\tilde{\varepsilon}_{+}))/(\tilde{\varepsilon}-\tilde{\varepsilon}_{+})$, see Fig. 1(a). Solitary waves describing the extension pulses, correspond to homoclinic trajectories starting and ending at $\tilde{\varepsilon}_{+}$. Periodic trains of such pulses correspond to closed trajectories encircling the degenerate center at $\tilde{\varepsilon}=1$. As $V\to V_{*}$ homoclinic trajectories become heteroclinic and the solitary waves turn into kinks; at $V=V_{*}$ we have $S_{1}=S_{2}$ in Fig. 1(a) and therefore the associated macroscopic discontinuity is dissipation free Truskinovskii (1982). The structure of these solutions is similar to the one in flocking models Caussin et al. (2014); Solon et al. (2015) modulo the fact that here we omit the explicit description of inflow and outflow of energy. Note that the contraction pulses, which exist in the complimentary range of parameters $V_{*}<V<V_{**}$, correspond to homoclinic trajectories starting and ending at $\tilde{\varepsilon}_{-}$.

In view of the piecewise linear nature of our model, the homoclinic solution of (3) can be written explicitly

where $z=V/\sqrt{12(V^{2}-1)}$ and $\tanh(d/z)=-(1-\lambda)^{-1}.$ Defining the amplitude of the pulse as $A=\max\tilde{\varepsilon}(\eta)-\min\tilde{\varepsilon}(\eta)$, we find for extension pulses that $A=\lambda+(1-\lambda)/\cosh(d/z)$. The corresponding explicit solution for contraction pulses is presented in (SOM, ).

The functions $d(V)$ and $A(V)$ for both types of pulses are shown in Fig. 2(a,b). The two families are separated by the critical value of parameter $V=V_{*}$ where the solitary waves take the form of infinitely separated kinks. At this point the parameter $d$, playing the role of the correlation length, diverges even though the pulse amplitude remains finite (taking the value $A=2$). In the limits $V\to 1$ and $V\to V_{**}$ we obtain sonic waves in passive and active states, respectively; note that the passive limit is singular, see (SOM, ) for more details. The typical functions $\tilde{\varepsilon}(\eta)$ for different values of $V$ are shown in the insets in Fig. 2. We emphasize that only the near-critical pulses have a physiologically realistic rectangular form.

Using the relation $\tilde{\varepsilon}=-Vd\tilde{u}/d\eta$, we can obtain the amplitude of the displacement increment culminating the passing of a pulse, for instance, in the case of a stretching pulse $\Delta\tilde{u}=2d\lambda$. As a rough description of a peristaltic wave train, represented by a succession of $N$ such pulses, we can write $\Delta\tilde{u}=2dN\lambda$.

To construct the actual periodic solution we can use the representation $\tilde{\sigma}(\eta)=\tilde{\varepsilon}(\eta)+\tilde{\sigma}_{a}\sum_{j}\text{rect}(\eta-jD/(2d_{p})),$ where it has been assumed that each pulse has a half width $d_{p}$ and the whole active lattice has the period $2D$. The matching conditions are now $\llbracket\tilde{\varepsilon}\rrbracket=\llbracket d\tilde{\varepsilon}/d\eta\rrbracket=0$ and $\tilde{\tilde{\varepsilon}}(jD\pm d_{p})=1$. While the whole solution can be again written explicitly, see SOM , here we only present the transcendental equation for the correlation length $d_{p}$ which is now regularized by the fixed system size $D$: $(1-\lambda)\tanh(d_{p}/z)=\tanh((d_{p}-D)/z).$

The obtained family of solutions is parameterized by $V$ and incorporates both stretching (active) and contraction (passive) pulses. To distinguish between the two it is convenient to re-define the half-width as $d=\min{(d_{p},D-d_{p})}$ and the amplitude as $A=|\tilde{\varepsilon}(0)-\tilde{\varepsilon}(D)|$. The resulting functions $d(V)$ and $A(V)$ are shown in Fig. 3.

The obtained periodic solutions can be used to model peristaltic locomotion. Suppose that the organism generates a periodic train of stretching pulses that propagate rearward with a velocity $V$ so that when a pulse reaches the tail, another one is initiated at its head with a fixed delay controlled by the parameter $D$. Observations show that when such pulse moves towards the tail, the latter fattens and gets anchored due to local increase of friction. With such an anchor present, the locomotion naturally occurs into the direction opposite to the direction of the pulse. We mimic such motility pattern in Fig. 4 (left), where the pulse is taken from the range $V<V_{*}$.

Suppose that the motion of the animal is of stick-slip-type with each advance corresponding to passing of a single pulse producing the forward displacement of the head $\Delta\tilde{u}=2d_{p}\lambda$. Since the next pulse arrives after the time $2D/V$, the mean translational velocity of the system is $\overline{v}=\Delta\tilde{u}V/(2D)$. Note that the simultaneous propagation of multiple peristaltic pulses along the animal body is also a possibility which we do not consider here.

The branch of contraction pulses corresponding to $V>V_{*}$ also generates a motility pattern shown in Fig. 4 (right), In this case the ’fattening’ and the resulting anchoring takes place around the pulse. Such regimes, driven by passive pulses in otherwise active medium, are not realistic due to the necessity to maintain active state throughout the whole body of the organism.

The behavior of the function $\overline{v}(V)$ for both kinds of motility ($V\lessgtr V_{*}$ ) is shown in Fig. 4 (center). We emphasize that around the critical point $V=V_{*}$, the macroscopic velocity $\overline{v}$ behaves singularly: it varies over a broad range around a single value of the control parameter. The corresponding individual pulses take the form of elongated rectangles which is one of the most characteristic features of peristaltic waves. Since the width of these rectangular pulses, and therefore the resulting motility velocity, can vary significantly, being positioned near such critical point, can help an organism to adapt its responses. The implied anomalous sensitivity to controls can facilitate optimal behavior in complex physiological conditions and would then be highly functional.

Given that the biological systems exhibiting peristalsis are often segmented, the question arises whether our oversimplified (quasi) continuum model (2) adequately represents the dynamics of its discrete prototype (1). To answer this question we now briefly consider the traveling wave solutions of the original discrete system.

We maintain the same normalization and use again the ansatz $\tilde{\varepsilon}_{n}(t)=\tilde{\varepsilon}(\eta)$, where $\eta=(nh-Vt)/h$. The discrete strain field $\tilde{\varepsilon}(\eta)$ satisfies the equation

with $\tilde{\sigma}(\eta)=\tilde{\varepsilon}(\eta)+\tilde{\sigma}_{a}H(-\eta)$. Solitary wave can be built directly from kinks which we therefore consider first.

Suppose that a switching point where $\tilde{\varepsilon}(0)=1$ is placed at $\eta=0$. Application of the Fourier transform to (5) allows one to obtain an explicit solution:

where $L(k)=\omega^{2}(k)-(kV)^{2}$, $\omega^{2}(k)=4\sin^{2}\left(k/2\right)$ and $\tilde{\varepsilon}_{+}=\tilde{\varepsilon}(\infty)$. Computing the integral in (6) we obtain the representation of this solution in the form of infinite series

where $Z^{\pm}=\left\{k\,:\,L(k)=0,\,\pm\mathrm{Im}\,k>0\right\}$ and $\tilde{\varepsilon}_{-}=\tilde{\varepsilon}_{+}+\lambda$. This solution exists only for the critical value of velocity $V=V_{*}$, see SOM for details, which suggests that admissible kinks in this model must necessarily satisfy the dynamic Maxwell condition.

To construct solutions representing stretching pulses we must solve (5) with the stress given by the same rect function as in the continuum case and use the same matching condition $\tilde{\varepsilon}(\pm d)=1$. The linearity of the system at fixed $d$ hints that the solitary wave solution can be obtained as a linear combination of two kinks $\tilde{\varepsilon}(\eta)=\tilde{\varepsilon}_{k}(\eta-d)-\tilde{\varepsilon}_{k}(\eta+d).$ The nonlinear relation $d=d(V)$ can be then found from the matching condition $\tilde{\varepsilon}_{k}(0)-\tilde{\varepsilon}_{k}(2d)=1$. By inverting the equation $(e^{-ik}-1)\hat{u}(k)=\hat{\varepsilon}(k)$ in the Fourier space we can also find the displacement field $u(\eta)$, see SOM , and compute the total displacement which is again $\Delta\tilde{u}=2d\lambda.$ The typical dependencies $d(V),$ and $A(V)$ for stretching pulses are illustrated in Fig. 5, where they are compared with the corresponding results of the continuum theory with $h=1$; the comparison confirms the overall adequacy of the continuum approximation despite the choice of a finite value for our ’small’ parameter.

We tested the numerical stability of the obtained discrete pulses by solving a range of initial value problems. An example of a stable propagation of a rectangular stretching pulse is shown in Fig. 6(a) for the chain with $1000$ springs. We used initial conditions $\tilde{u}_{n}(0)=0,\,\dot{\tilde{u}}_{1}(0)=v_{0}$ and $\dot{\tilde{u}}_{n}(0)=0,\,n>1$ and assumed free ends. The time evolution of the displacement field, shown in Fig. 6(b), illustrates the creation of the main driver of peristaltic motility, the finite displacement behind the propagating pulse.

To conclude, we developed a model of a dynamic passive-to-active transformation taking place in the front of a steadily moving pulse with the corresponding reverse transformation taking place in its rear. The resulting solitary wave solutions for the active medium were extended as periodic trains and used to model the peristaltic mode of self propulsion. We found that at the critical value of parameter the model generates singular regimes with diverging effective correlation length and argued that such criticality may be functional. Our results can be used for biomimetic reproduction of worm-like motion in applications ranging from endoscopic diagnostics Stefanini et al. (2006) to pipeline inspection Yamashita et al. (2011). While in the existing robotic systems activity is usually imitated by globally synchronized distributed actuators Agostinelli et al. (2018); Miller and Dunkel (2020); Jiang and Xu (2017) our focus on local mechanical feedback in muscle-type soft deformable materials opens new avenues in the modeling of peristalsis, see also Pehlevan et al. (2016); Umedachi et al. (2016); DeSimone (2020). The proposed prototypical model can serve only as a proof of concept and future work allowing one to make quantitative predictions, should incorporate energy supply and dissipation and account for realistic 3D geometry.

The authors thank G. Mishuris, P. Recho, A. Vainchtein for helpful discussions and acknowledge the support of the French Agence Nationale de la Recherche under the grant ANR-17-CE08-0047-02.