Model of metameric locomotion in active directional filaments

Locomotion in segmented animals, such as annelids and myriapods (centipedes and millipedes), is generated by a coordinated movement known as metameric locomotion, which can be also implemented in robots designed to perform specific tasks. We introduce a theoretical model, based on an active directional motion of the head segment and a passive trailing of the rest of the body segments, in order to formalize and study the metameric locomotion. The model is specifically formulated as a steered Ornstein-Uhlenbeck curvature process, preserving the continuity of the curvature along the whole body filament, and thus supersedes the simple active Brownian model, which would be inapplicable in this case. We obtain the probability density by analytically solving the Fokker-Planck equation pertinent to the model. We also calculate explicitly the correlators, such as the mean-square orientational fluctuations, the orientational correlation function and the mean-square separation between the head and tail segments, both analytically either via the Fokker-Planck equation or directly by either solving analytically or implementing it numerically from the Langevin equations. The analytical and numerical results coincide. Our theoretical model can help understand the locomotion of metameric animals and instruct the design of metameric robots.

Introduction. -Active particles are systems far from equilibrium that utilize energy sources to perform nonthermal motions [1][2][3][4]. Examples of active particles include both living systems such as cytoskeletal filaments [5][6][7], microorganism colonies [8,9], and bird flocks [10], as well as synthetic systems such as self-propelled colloids [11] and driven granular matter [12]. Among these active systems one can also include animals with metameric structure that enables them to perform a special type of gait referred to as metameric locomotion. Limbless and manylegged invertebrates can generate a coordinated gait via a collective movement of segments or appendages, not dissimilar to a metachronal wave in cilia arrays [13], that efficiently propels the body forward. Inspired by the crawling creatures, robots have been developed that simulate the metameric locomotion in order to perform tasks that otherwise traditional robots would accomplish only inefficiently [14,15]. While presently the studies of metameric locomotion focus mostly on the experimental observation of living worms and millipedes [16][17][18][19], as well as related biology inspired robot engineering problems [15,[20][21][22][23][24][25], theoretical models highlighting the dynamical properties of the motion trajectory are lacking. In what follows we formulate an active directional filament model to simulate the metameric locomotion and obtain analytical results describing its properties, including the probability density and physical correlators related to the trajectory of the motion.
We start from the continuity constraint for the filament tracing the body frame as one of the defining features of the model describing metameric locomotion. The shape of the moving filament should be not only continuous in position, but also with continuous orientation and curvature fields, as there are indications that nematodes and other metameric animals trace out locomotion trajectories with continuous curvature [17], likely related to the fact that the directional steering itself is continuous, possibly corresponding to lowest energy consumption. The continuous curvature constraint already disqualifies variants of the acp-1 tive Brownian model [26,27], where the orientation of the active particles changes with time continuously but is not differentiable, corresponding to a discontinuous curvature field. We thus resort to a higher order model dynamics, with continuous position, orientation as well as curvature fields. Without any sensorial feedback the metameric locomotion can be presumed to persist as rectilinear, however, the constant influx of environment information and the inaccuracy in the coordinated response can be listed as sources of the unavoidable curvature noise. We therefore in addition assume that the curvature can be reasonably described as undergoing an Ornstein-Uhlenbeck process, possibly steered by external forces. In addition, due to the connectivity of the metameric body the locomotion shows a trailing property, where distal body sections follow the trajectory of the head. For the same reason the body filament can be considered as inextensible globally as well as locally, leading to a filament velocity with a constant norm. All the listed constraints of the metameric locomotion model can be subsumed by describing it as an active directional filament motion, generated by an active head with curvature noise and a passive body trailing the trajectory of the active head.
The model constraints laid out lead to a set of equations of a higher order Langevin form that allows us to obtain the probability density of and the correlators characterizing the metameric locomotion analytically. Two paths are then taken to solve the model: (i) obtain the correlators, i.e., the mean-square orientational fluctuations, the orientational correlation function and the mean-square separation between the head and tail of the body filament, directly from the solution of the Langevin equations; (ii) obtain the correlators via the solution of the corresponding Fokker-Planck equation, which gives the probability density of the head particle, from which one can calculate the probability density for the whole moving filament. Since the solution of the Fokker-Planck equation carries the complete information regarding the metameric locomotion, we will mainly pursue the calculation of the second path, but also establish that the results from the Fokker-Planck method and those stemming directly from the Langevin equations coincide. Moreover, we will numerically integrate the Langevin equations and obtain the running averages of the trajectories and show that they agree well with the analytical results, consequently validating both the analytical as well as the numerical results.
Model. -We now propose the analytical form of the active directional filament model for metameric locomotion. The active filament lives on a two-dimensional infinite flat plane. It consists of an active head and a passive trailing body. The motion of the head maintains a velocity with constant norm v 0 . The length of the filament is L. Denote the position of the body segment at an arclength distance l from the head as rl(t), where t is time and the dimensionless arclength is defined asl = l/L with 0 <l ≤ 1. Then the position of head and tail are r 0 (t) and r 1 (t), respectively. Denote the time for the filament to move a distance equal to its body length as T = L/v 0 . Since the rest of the body segments are trailing the trajectory of the head, this implies that for t > T In the main text, we will focus on the large time regime t > T . The marginal case with t ≤ T will be trivial to discuss after we obtain the results of the case with t > T [28]. As long as no ambiguity rises, we will omit the subscript of the variables pertinent to the active head for simplicity. The motion of the active head with a constant velocity norm follows the dynamical equations [29] where t(t) and n(t) are, respectively, the tangential and normal vectors of the head motion, while κ(t) is the instantaneous curvature.
In two-dimensional space, the tangential vector can be rewritten as t = (cos θ(t), sin θ(t)) with θ(t) being the orientational angle of the tangential vector of the head segment, see Fig. 1. Then it follows that By assumption the curvature of the head segment undergoes a steered Ornstein-Uhlenbeck process where β > 0 is curvature decay constant and ξ(t) is white noise with g being the constant noise amplitude. For generality, we also added an external force term, −∇U , derived from a position and orientation dependent steering potential U (r, t), with γ being a positive coupling coefficient. From the above Langevin equations, we can derive the corresponding Fokker-Planck equation. Note that in the p-2 absence of the steering potential the correlators related to position r(t) can be calculated from the orientation data, see Sec. Correlators, so it suffices to solve the Fokker-Planck equation for the probability density as a function of only orientation and curvature. Assume the probability of finding at time t the active head in the orientation interval [θ, θ + dθ] and curvature interval [κ, κ + dκ] is P 0 (θ, κ, t|θ 0 , κ 0 , 0) dθ dκ, where θ and κ are, respectively, the orientation and curvature at time t, and θ 0 and κ 0 are, respectively, the orientation and curvature at time t = 0. Then we end up with the following Fokker-Planck equation which coincides with a two-dimensional Ornstein-Uhlenbeck process and is thus analytically solvable [28]. This is of course only true when there is no steering potential or the steering potential depends only on orientation [28].
Probability density. -Assume the active head initially has orientation and curvature distribution where θ(t) and κ(t) are, respectively, the mean orientation and curvature dependent on time and σ(t) = σ θθ σ θκ σ κθ σκκ , which is the variance matrix with elements Then the solution of Eq. (7) is given by [28,30] Clearly, the probability density depends explicitly on the orientation, curvature and time. It is noteworthy that the range of orientation in the probability density function is (−∞, ∞) instead of being wrapped into the interval [0, 2π). The reason is that the state of a filament is not a periodic function of the orientations of individual segments. Specifically, different from a rod, the filament can coil into spiral structures when the orientation of the head varies to exceed the range [0, 2π). We now remark on the behavior of the probability density Eq. (10) at large time limit. When t → ∞, σ θκ → gv 0 /(2β 2 ) and σ κκ → g/(2β). The mean curvature vanishes while the mean orientation remains finite. The orientational component of the variance matrix diverges linearly with increasing time while other components remains finite. It indicates a localized distribution of curvature around zero and an extended distribution of orientation around a finite value at large time.
The probability density of the whole filament, trailing the head motion, can be readily obtained from the probability density of the active head. The probability density of finding a segment of a single filament at time t with orientation θ and curvature κ is Recall that the passive segments of the filamentous body follow the trajectory of the active head. The state of the passive segments therefore lags behind the state of the head, and we conclude that The above integral can be performed numerically. Let us define the dimensionless variables by using the following characteristic quantities: we define the characteristic time t c = 1/β, length l c = v 0 /β, curvature κ c = (g/β) 1/2 and orientational angle θ c = (gv 2 0 /β 3 ) 1/2 , yielding the dimensionless timet = t/t c , lengthl = l/l c , curvatureκ = κ/κ c , orientational angleθ = θ/θ c and noiseξ = ξt c /κ c . Shown in Fig. 2 are the probability densities of finding a segment of a single filament dependent on orientation and curvature of the segment at different times with varying initial orientations and curvatures. The timet = 4 in Fig. 2(a1, b1), whilet = 10 in Fig. 2(a2, b2). The initial orientationθ 0 = 0 and curvatureκ 0 = 0 in Fig. 2(a1, a2), whileθ 0 = 1 andκ 0 = 1 in Fig. 2(b1, b2). The length of the chain is fixed to beL = 1 in all the subfigures. The probability density is partially heterogeneous and localized at timet = 4 [ Fig. 2 (a1)] and becomes completely heterogeneous at a larger timet = 10 [ Fig. 2 (a2)]. At large times, due to the diverging orientational component and finite curvature component in the variance matrix, the distribution of orientation is extended, while the distribution of curvature remains localized. Different initial orientation and curvature [ Fig. 2(b1, b2)] affect only the positions of the mean orientation and curvature. The shapes of the distributions are, however, not altered.
Correlators. -Now that we have analytically solved the Fokker-Planck equation and obtained the probability density, we can calculate correlators such as mean-square orientational fluctuations, orientational correlation function and mean-square separation. The correlators depend p-3 in general on two parameters: time as well as position along the filament.
We first discuss generally the correlators of the active head at time t 1 and t 2 . Since the state of the distal segments of the filament lags behind the state of the active head, the correlators between the head and tail at time t 2 are just a special case with The mean-square separation is related to the orientational correlation function by while the orientational correlation function can be calculated from the mean-square orientational fluctuations. Note first that t(t 1 ) · t(t 2 ) ξ = cos(θ 2 − θ 1 ) ξ and θ 2 − θ 1 = v 0 t2 t1 κ(t ′ )dt ′ , where θ 2 ≡ θ(t 2 ) and θ 1 ≡ θ(t 1 ). Since κ undergoes the Ornstein-Uhlenbeck process, θ 2 − θ 1 is a Gaussian variable whose characteristic function can be readily obtained, and the orientational correlation function t(t 1 ) · t(t 2 ) ξ follows as just the real part of the characteristic function of θ(t 2 ) − θ(t 1 ), that is The calculations of the mean-square separation and the orientational correlation function then reduce to the calculations of the mean and the mean-square orientational fluctuations, which can be obtained either from the Fokker-Planck probability density or directly from the Langevin equations.
For the curvature, the following result containing scaled time-transformed Wiener process is used To integrate the orientational angle we use the discretized Euler version of dynamics just as for the case of the position vector dynamics Shown in Fig. 3 are the 3D plots of the correlators, i.e., the mean-square orientational fluctuations [ Fig. 3(a)], the orientational correlation function [ Fig. 3(b)], and the mean-square separation [ Fig. 3(c)] between the head segment and the trailing tail segment of the filamentous body. Fig. 4 displays the cross sections of the 3D plot in Fig. 3 at fixed lengthL = 1 [ Fig. 4(a, b, c)] and at fixed time t = 10 [ Fig. 4(d, e, f) with the analytical ones, consequently validating both the analytical as well as the numerical results. Note that the mean-square orientational fluctuations Eq. (18) does not depend on initial orientation. The initial curvature only affects the correlators at small times. At large times, as shown in Fig. 4(a, b, c), the initial curvature is irrelevant because the exponential factor in the term containing the initial curvature in Eq. (18) decays to zero at large times. θ c = (gv 2 0 /β 3 ) 1/2 characterizes the orientational fluctuap-5 tions during the characteristic time t c = 1/β. Large curvature diffusion constant g/2 and norm of active velocity v 0 both promote the orientational fluctuations, while higher ability to maintain a straight body (large β) suppresses the orientational fluctuations. Consequently, higher θ c results in higher mean-square orientational fluctuations, lower orientational correlation, and smaller mean-square separation between the head segment and the tail segment of the filamentous body. From Fig. 4(f), one can see that the mean-square separation depends quadratically on the length at small length and linearly at large length. At medium lengths, there is a crossover from the ballistic to the diffusive behavior of the mean-square separation.
Conclusions. -In summary, we defined an active directional filament model to describe the metameric locomotion of both segmented animals as well as segmented robots. As the curvature of these segmented bodies should be continuous, we designed a model with this property explicitly enforced, so that by necessity it generalizes the active Brownian motion model, which does not preserve a continuous curvature. The crawling creatures and wobbling robots are modeled as active directional filaments, being composed of otherwise identical body segments with the sole proviso, that only the head segment is active with its curvature undergoing a steered active Ornstein-Uhlenbeck stochastic process. The rest of the segmented body is assumed to passively trail after the active head so that its dynamical state simply lags by various amounts from the state of the head segment.
We obtain the probability density of the active head by solving the appropriate Fokker-Planck equation and from this we calculate the probability density for the whole filament by integrating the probability density of the active head along the rest of the segmented filamentous body. The probability density of the whole filament is found to be partially heterogeneous at small time scales but becomes completely heterogeneous at large time scales with extended distribution of orientation and localized distribution of curvature. The various correlators of the dynamical state of the segmented filament are calculated by three different methods, i.e., analytically from the Fokker-Planck equation, analytically from the Langevin equations, as well as numerically by integrating the Langevin equations.
The results obtained by three different methods coincide. Initial orientation and curvature are shown to be irrelevant at large time in the expression for the meansquare orientational fluctuations, the orientational correlation function and the mean-square head-to-tail separation. A characteristic orientation emerges, depending on the activity, diffusion of curvature and curvature time decay constant, which determines the mean-square orientational fluctuations and consequently the orientational correlation function and the mean-square head-to-tail separation. Our theoretical model sheds new light on the metameric locomotion of segmented animals but can also be useful in the design of metameric, segmented robots [14,15,21,25].

Model of metameric locomotion in active directional filaments
Supplementary Material

I. DERIVATION OF PROBABILITY DENSITY FROM FOKKER-PLANCK EQUATION
Here we solve the probability density P 0 (θ, κ, t|θ 0 , κ 0 , 0) fulfilling the Fokker-Planck equation However, we will not solve the above Fokker-Planck equation directly. We first observe the Langevin equations, from which the Fokker-Planck equation is derived where Rearrangement of Eqs. (S3) into matrix form gives us d dt which is just two-dimensional Ornstein-Uhlenbeck process. Denote We have the Fokker-Planck equation for two-dimensional Ornstein-Uhlenbeck process The solution of Eq. (S7) is well known [1] and given by where We now use spectral decomposition to write the above equations explicitly. We express B in biorthogonal vectors (repeated indices are summed) Then (S11) Finally we obtain explicitly the mean orientation and curvature and the variance matrix which completes the solution of the Fokker-Planck equation Eq. (S1). Now that we have the joint probability density dependent on both orientation and curvature, we can readily obtain the marginal probability density dependent separately only on orientation and curvature. The marginal distribution of multivariate normal distribution can be obtained by just dropping irrelevant parts, which can be easily proved. Hence we obtain the following two marginal probability distributions where (S16)

II. DERIVATION OF PROBABILITY DENSITY WITH ORIENTATIONALLY STEERING FIELD
The presence of orientationally steering field can profoundly alter the behavior of the probability density of the active head at large time. Assume the motion of the active head is orientationally steered by a harmonic potential where ǫ is a positive constant and θ p represents the preferred orientation. Then the Langenvin equations become Denote λ 1,2 = 1 2 (β ± β 2 − 4ǫv 0 ). The solution of corresponding Fokker-Planck equation is like Eq. (S8) where θ ′ (t) and κ ′ (t) are, respectively, the mean orientation and curvature and σ ′ is the variance matrix with elements being Note the real parts of λ 1,2 are always positive. Hence when t → ∞, The presence of orientationally steering field renders the mean orientation vanish and the variance of orientation finite at large time compared to the case without orientationally steering field. It is noteworthy that the transition of orientational variance from diverging to finite is continuous at ǫ = 0.

III. FILAMENT PROBABILITY DENSITY AND CORRELATORS AT SMALL TIME
In the main text, we focus on the large time regime with t > T , where T is the time for the filament to move a distance equal to its body length. We now discuss the marginal case with t ≤ T . A metameric animal may drop accidentally into a new environment or a metameric robot can be set manually to start moving. When t ≤ T , part of the distal segments remain at the initial state, having not become the lagged state of the active head yet.
The mean and the mean-square orientational fluctuations between the head and tail (θ(t) − θ 1 (0)) ξ and (θ(t) − θ 1 (0)) 2 ξ are again trivially related to θ(t) − θ 0 ξ and (θ(t) − θ 0 ) 2 ξ . Therefore, the calculations of the mean-square orientational fluctuations, the orientational correlation function and the mean-square separation between the head and tail at time t ≤ T are all reduced to the calculations of θ(t) − θ 0 ξ and (θ(t) − θ 0 ) 2 ξ . From the main text, we know which completes the calculations of the correlators in the marginal case with t ≤ T .