Coordination in Coupled Arrays of Stiff Filaments—Modelling and Simulation

: We present the mechanical model of an array of elastic ﬁlaments and simulate the response with different mechanical couplings. This class of systems is inspired by robust and elegant solutions for locomotion mechanics that have emerged in several small-scale biological entities in the form of beating protrusions, such as cellular cilia and eukaryotic ﬂagella. The collective dynamics of cilia arrays reveals important features such as array alignments, two-phase asymmetric beating of individual ﬁlaments, and the emergence of metachronal coordination, which make them suitable for bio-inspired terrestrial and aquatic locomotion. The model presented here is the basis for further developments towards the design of terrestrial and aquatic locomotion systems for general purpose robotic devices.


Introduction
Several small-scale biological entities have evolved highly adaptable, low specialization locomotion mechanics systems in the form of beating protrusions, such as cellular cilia and eukaryotic flagella [1,2]. The fundamental structure of many eukaryotic cilia and flagella is the axoneme, which is a regular cylindrical arrangement of elastic filaments connected by motor proteins, which act as active elements from the perspective of actuation and induce beating patterns at the core of locomotion actions for small-scale motility [3,4].
The biomechanics of axonemal beating protrusions is extensively studied in the context of solid-fluid interaction, with the investigation often focused on the emergence of beating patterns induced by the coupling with an external flow [2,3,[5][6][7]. Fluid interactions and induced beating evolution in small-scale protrusions are of particular interest in mathematical biology, since they are at the core of micro-swimming and also relate to important physiological and developmental phenomena such as mucous flow control and protection against infections, as well as embryonic development [8,9]. From the perspective of locomotion in fluids, the collective dynamics of cilia arrays has important implications as it reveals important aspects such as array alignments, two-phase asymmetric beating of individual filaments, and the emergence of metachronal coordination [10][11][12]. Metachronal coordination in arrays of coupled protrusions consists of synchronization with constant phase shift between contiguous elements, resulting in travelling waves through the array, ultimately associated with locomotion mechanics.
Recent work in [13] proposed that emerging metachronal waves in arrays of stiff filaments are also fundamentally related to mechanical coupling through the base in addition to hydrodynamic coupling, suggesting that the emergence of metachronal waves can be achieved by pure basal mechanical coupling. This has important implications for terrestrial locomotion, since from a design point of view, arrays of flexible protrusions are a robust solution as they can adapt and morph to a variety of terrains, and from the dynamics point of view, metachronal coordination and travelling waves are suitable to induce motion when the protrusions engage with a substrate/terrain. Once the dynamics of this class of systems can be predicted, important optimization considerations for robotics applications can be addressed, such as minimal actuation constraints and induced coupling topologies to minimize energy consumption and maximize autonomy, for example.
Here, we present the mechanical model of an array of elastic filaments connected through an elastic medium (base) and consider different basal elastic couplings, namely forces and couples. Each filament is considered as a stiff polymer with an active internal force, consistent with axonemal structures. This class of mechanical objects is suitably modelled by using constrained Euler elastica, resulting in coupled nonlinear PDEs describing the evolution of a kinematic descriptor for the geometric shape, and a Lagrange multiplier that enforces inextensibility. The response is simulated for different coupling schemes, hinting at possible designs for walker-type locomotion mechanics based on their collective dynamics. Further developments are focused on general purpose robotic devices, which are attractive solutions for a variety of applications that include environmental monitoring and protection, civilian and military defence, and generally hazardous missions [14][15][16][17][18].

Stiff Filaments with Internal Action
A cross-sectional image and a schematic of an axonemal structure are shown in Figure 1, reproduced from [19] with the author's permission. Axonemes are the fundamental structures at the core of locomotion of small-scale structures like cells and bacteria, in which it is often achieved through the beating of cilia and flagella protrusions. Axonemal beating is the complex result of the co-action of internal electrochemical motors and of external forces resulting from hydrodynamic interactions and from forces transmitted through the base [3]. The resulting nonlinear mechanics allows for the emergency of a variety of possible patterns, which confer robustness and adaptability to this type of structure with respect to the locomotion function and suggest the possibility that the solution obtained at the evolutionary scale could be fundamental and therefore desirable in engineered systems designed to engage unknown environments.  [19] with the author's permission). A simplified mechanical model for an axonemal beating structure in two dimensions is schematised in Figure 2, in which two filaments are separated by approximately a constant distance 2a [19]. The arc length s spans the neutral axis of the system with position vector r(s), and points on the boundary filaments are described by the position vectors r 1 (s) = r(s) + an(s) and r 2 (s) = r(s) − an(s), where n is the unit normal to the neutral axis. The bending of the system is driven by an internal system of shear forces f (s) [4]. In biological active filaments, the internal forces f (s) result from the concurrent effect of dynein motors (active) and passive elastic connection across microtubule doublets. The bending stiffness of the filaments is related to an internal bending moment that works on the curvature. In an engineered realization of the same systems, passive and active links can be reproduced respectively by distributed motors and by viscoelastic links. At every point, the pair of equal magnitude and opposite sign forces f is equivalent to an internal couple (bending moment) of magnitude a f . The internal forces induce relative sliding of the filaments, which is related to the slope of the neutral axis by [19]: Therefore: where κ is the curvature of the neutral axis. The equations of motions can be derived from an action principle [20]. Let r(s, t) be generalized coordinates describing the neutral axis of the system. We have: where ψ(s, t) is the local tangent angle to the neutral axis, with respect to a global coordinate frame. If L and R are respectively the Lagrangian and the Rayleigh dissipation function of the system, then the equations of motion are: The Rayleigh dissipation function is quadratic in the generalized velocities: with c being a constant dissipation constitutive parameter. We model this system as a nonlinear rod, governed by curvature with a distributed system of internal shearing forces f (s, t). The potential energy of the system is given by [20]: where λ is a multiplier that enforces the local inextensibility condition r · r = 0. For small deformations, the curvature is the second derivative of the transverse displacement, and the potential energy reduces to the one of an inextensible Euler beam. Since s is an arc length, the following intrinsic relations define the tangent and the normal vectors to the parametric curve r: The kinetic energy of the system is given by: where the superimposed dot means the time derivative. The variation of the potential energy gives: where f = f is the magnitude of f , and: By using the relation δκ = δ(n · r ) = n · δr [1] Appendix A and by two integrations by parts, we obtain: Boundary terms conjugate to δr are boundary bending moments, and boundary terms conjugate to δr are external forces, with explicit tangential and normal components. By introducing: we can see that this quantity is the integrated tangential component of the forces acting on the filament, and therefore, it acts as a physical tension. This gives also an interpretation to λ, which is the portion of the physical tension necessary to enforce the inextensibility condition. The equations of motion for a single filament are: with boundary conditions: where F 0 , F , M 0 , M are external forces and moments applied at the boundaries. In order to obtain evolution equations for ψ and τ, we first consider the following: Therefore: where ν = k M κ − a f . By carrying out the spatial derivatives: Therefore, by projecting along n and t, we obtain the governing equations: In order to close the system, we need a constitutive model for f , which is in general comprised of a passive and an active component. Different constitutive models have been considered in the literature [1,4] to fit observations about flagella beating, specifically focusing on passive elastic coupling and on the active action related to dynein motors. Here, we do not specify the active component in order to isolate the effect of the basal coupling, and we consider the passive component to be described by the linear relation: with the elastic constant k f characterized as: where E f and A f are respectively Young's modulus of the filaments' material and cross-sectional area. At the boundary s = 0, we consider an elastic coupling with the body of the robotic device, which will be discussed in the following subsection. For the purpose of this study, we are interested in the dynamics induced by the basal coupling, and therefore, we consider the extremes at s = to be free. The boundary conditions will be formalized below.

Weak Form and Finite Element Formulation
The dynamics of the system of coupled stiff filaments is simulated by using finite elements implemented in a code written in Mathematica c 12. The finite element solution is based on the following weak or variational form of the system of governing equations.
For i = 1, . . . , N, define the arc-length s i ∈ [0, ] spanning the neutral axis of the ith filament. For s i ∈ (0, ), that is in the open interval spanning the neutral axis, the plane dynamics of the ith filament is governed by the following set of equations adapted from (20): Therefore, there are in total 2N fields (ψ i , τ i ) to which we associate the 2N test function fields (ψ i ,τ i ), which can be interpreted as the variations of the states in the sense of the calculus of variations, or otherwise, they can be interpreted as test functions in the framework of numerical approximation methods. The test functions satisfy homogeneous boundary conditions assigned to the trial solutions ψ i and τ i [21]. Integration with respect to s i over the domain [0, ] after multiplying the expressions in (23) respectively byψ i andτ i gives: Integration by parts of different terms (see Appendix A) gives: From the boundary terms, we establish the following variational duality relations: Since ψ i = κ i is the curvature, by using the physical quantities introduced in Section 2, the duality relations can be rewritten as: which shows that the dual of the tension τ is the tension per unit length 2κM − τ , where M = k M κ − aF is the moment, which is also the integral of the normal force k M κ − a f .

Boundary Conditions and Coupling
At the free end of each filament, all the dynamic quantities are zero, and therefore, the boundary terms in (25) evaluated at s i = are set to zero. On the other hand, the coupling between filaments through the base is expressed through the boundary terms evaluated at s i = 0.
The idea of basal coupling is schematised in Figure 3, which shows an array of stiff filaments anchored to a base that provides the coupling. The array includes N filaments labelled with the integer i = 1, 2, . . . , N, each described by a curvilinear abscissa s i ∈ [0, ], where is the common length. In Figure 3, the field ψ i (0, t) is the rotation of the ith filament at the point of attachment to the base. To model rigid coupling, one would impose kinematic constraints of the type ∑ N j=1 a j ψ j (0, t) = 0 or ∑ N j=1 b j κ j (0, t) = 0 for some set of coefficients a j or b j . On the other hand, force coupling through the base is appropriately modelled through boundary forces and moments that are exerted on each filament as specified by the boundary terms at s i = 0 in (25). The appropriate force terms are identified by the duality relations (27). Based on these relations, we define the generalized normal force, the generalized moment, and generalized tangential force respectively as the duals of the rotation ψ, the curvature κ, and the tension τ: Note that the just introduced generalized forces and moments involve derivatives of the original physical quantities introduced in Section 2 since the model (23) is obtained from (14) by taking a spatial derivative of the position vector r to express the system in terms of ψ and τ.
There are multiple possibilities to express the coupling among filaments through the generalized dynamic quantities just defined. Here, we consider the relatively simple possibility of modelling it as constitutive relations between generalized boundary forces and moments and boundary kinematic quantities. For rotation coupling, the boundary terms (in (25)) duals ofψ i are then given by: whereμ i (t) is an actuation for the ith filament, α µ is a constitutive parameter that plays the role of coupling strength, and σ ij (t) ∈ [0, 1] is the coupling between filament j and filament i. From the perspective of robotic applications, the actuation termsμ i act as inputs, and for the purpose of energy consumption, one would want to minimize the number of such terms, ideally having only one actuation with action propagated through the base to the non-actuated filaments. Moreover, more sophisticated analysis could involve the optimization through the coupling network topology induced by the distribution of coefficients σ ij (t), which play the role of weighted network links between nodes constituted by the extremes of the filaments attached to the base. The coupling in (31) can accommodate finite difference operators and therefore common lattice network structures. For example, for: the coupling has the structure of a second order central difference ψ i+1 (0, t) − 2ψ i (0, t) + ψ i−1 (0, t), which, upon introduction of the normalized coupling strengthk M = h 2 k M , where h is the distance between two filament attachment points on the base, allows studying a system with many filaments described by the limit: where Ψ is the field describing the rotation of the distribution of filaments at the interface with the base and S is a curvilinear abscissa. Similarly, for curvature coupling, the boundary terms duals ofκ i in (25) are given by: whereν i is an actuating generalized moment. If the rotation at the basis is prescribed, thenψ i (0) = 0, and similarly for the curvature.

Finite Element Discretization
To obtain numerical solutions for (25) where φ ψ1 , φ ψ2 , . . . , φ ψn , φ κ1 , . . . , φ ψn are Lagrange polynomials in H 2 [21], satisfying 0 (φ ψj ) 2 ds < ∞ and 0 (φ κj ) 2 ds < ∞ for all j, so that the integrals in (25) involving second derivatives ψ i are well defined. On the other hand, φ τ1 , . . . , φ τn is a set of Lagrange polynomials in H 1 [21], satisfying so that ψ i (s ik ) = a ik , ψ i (s ik ) = b ik , and τ i (s ik ) = c ik . Example plots of the adopted basis functions are given in the Appendix B. By substituting the projections (35) into the weak form (35), we obtain the finite element discrete form: where q = (a 1 , b 1 , . . . , a N , b N , c 1 , . . . , c N ) is the collection of all unknowns for the coupled system, of dimension 3nN −n, wheren is the total number of kinematic constraints, that is the number of constraints on ψ i , κ i , and τ i , i = 1, . . . , N. The 3nN −n × 3nN −n matrix M and the 3nN −n nonlinear operator C have the block structure: with the details given in Appendix C. The matrix K c has the block structure: and details can be found in Appendix C. The 3nN −n input vector U = (U 1 , . . . , U N ) depends on the actuation actions in (31) and (34), with entries given by: The operator F, mapping to a 3nN −n vector, accounts for the nonlinear dynamics of the filaments. It has the block structure: with entries determined by the integral terms of the weak form (25) as detailed in Appendix C.

Simulation Results and Discussion
Finite element simulation results were obtained for the following set of material and geometric parameters: = 1 m, ρ = 1 kg m −1 , ζ = 0.5 kg m −1 s −1 , k M = 10 kg × m 3 /s 2 , α ν = 20 N m −1 , α µ = 20 N, and α τ = 5 m −1 . The final time for numerical integration was t f = 20 s, and a system of four coupled filaments was considered.
Each filament was discretized by five uniformly placed nodes in [0, ], and the resulting nonlinear system of ODEs (37) was integrated with the Mathematica built-in nonlinear solver NDSolve, which uses multi-step methods with local error convergence control and adaptive time steps. The absolute and relative error control precision parameters in NDSolve are respectively AccuracyGoal and PrecisionGoal. The two parameters are set by default to half the machine precision, which in the case of the results calculated here was 16 digits.
The boundary conditions τ i (L, t) = 0, i = 1, . . . , N are essential or kinematic since the tension τ is a state variable and are enforced by removing the N degrees of freedom c in (t), i = 1, . . . , N. For the boundaries at s i = 0, we considered two types of coupling: 1. Rotation coupling, defined by the boundary conditions: and generalized moment µ i given by (31).
2. Curvature coupling, defined by the boundary conditions: and generalized normal force ν i given by (34).
In both cases, the kinematic boundary conditions remove N degrees of freedom. Therefore, the discrete system of nonlinear ODEs (37) has dimension N(3n − 2) = 13N for n = 5 nodes. In both cases, the set of boundary conditions at s i = 0 is completed with the coupling in the dual of the tension: Coupling coefficients σ ij are set as in (32) to induce direct coupling among contiguous filaments. Future research developments include modelling with evolving coupling coefficients. Actuation inputs were set to:μ with the same setting adopted forν i . All initial conditions were set to zero. The internal forces f i were also set to zero for the purpose of simulation to isolate the effect of the basal coupling. The plots in Figure 4 show the evolution of the four tip displacements r i ( , t), which exhibited an oscillatory behaviour induced by the actuation. For the same set of parameters, there were clear differences in the responses induced by the two types of couplings, where the transmissibility of oscillatory patterns was visibly more attenuated in the case of rotation coupling. The evolution of the tangential tension at the attachment with the base are shown in Figure 5, where oscillatory patterns emerge due to the actuation in terms of the generalized transverse force or generalized moment.  Three snapshots of the configurations of the array of four filaments are shown in Figures 6 and 7, respectively for rotation and for curvature coupling. Considering that the two coupling conditions were simulated with the same set of material and geometric parameters, it is interesting to note how different beating patterns can be generated by controlling the actuation and the coupling. In this specific case, rotation coupling resulted in nearly rigid body rotation of the filaments, whereas with curvature coupling, there was much more pronounced bending. This has implications for the design, hinting at the possibility of inducing different modes of operation by controlling the input and the type of coupling, without necessarily changing the material and geometric properties of the system.

Conclusions
A model for an array of stiff filaments with basal coupling was presented. The filaments' model was based on some features of arrays of cilia ubiquitously observed in small-scale mobile biological structures. These structures evolved desirable properties for robotic locomotion adapted to uncertain, unstructured environments, specifically in terms of robustness and versatility to different types of terrain and substrates.
The nonlinear model of the system was numerically solved with finite elements, and simulation results were compared for two types of basal coupling. Oscillatory beating patterns were induced with minimal actuation of one filament through coupling. For the same set of parameters, different coupling induced beating patterns with considerably different shapes, hinting at the possibility of controlling the operation modes through coupling modes.
Future work will address quantitatively the wave propagation associated with collective beating and synchronous and metachronal modes induced by basal coupling.
Funding: This research received no external funding.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A
Different integrations by parts are performed to reduce the maximum order of the spatial derivation to two, so that the finite element formulation involves polynomial basis functions of minimal order of regularity [21]:

Appendix B
For a mesh of five uniformly spaced nodes in the domain [0, 1], the following plots show the basis functions φ ψ and ψ κ and their derivatives. The functions are constructed by imposing the constraints (36). Basis functions φ τ are plotted in Figure A3. for i = 1, . . . , N, where: