## Abstract

In this paper, we address the identifiability of constitutive parameters of passive or active micro-swimmers. We first present a general framework for describing fibres or micro-swimmers using a bead-model description. Using a kinematic constraint formulation to describe fibres, flagellum or cilia, we find explicit linear relationship between elastic constitutive parameters and generalized velocities from computing contact forces. This linear formulation then permits one to address explicitly identifiability conditions and solve for parameter identification. We show that both active forcing and passive parameters are both identifiable independently but not simultaneously. We also provide unbiased estimators for generalized elastic parameters in the presence of Langevin-like forcing with Gaussian noise using a Bayesian approach. These theoretical results are illustrated in various configurations showing the efficiency of the proposed approach for direct parameter identification. The convergence of the proposed estimators is successfully tested numerically.

### 1. Introduction

Using kinematics to reveal internal properties of deformable objects moving in fluids is a difficult and long-standing topic in mechanics, e.g. [1]. One can, for example, think about the motion of a deformable cell in shear flow where various dynamical modes are observed depending on elastic capillary number [2, 3], internal properties of red-blood cells [4] or capsules [5, 6] in a creeping flow. Inferring internal mechanical properties has indeed already being pursued for flexible capsules from their deformation within a flow in a number of contributions [7–9] where the combination of experimental measurements and direct simulations can lead to estimating some internal parameters of the membrane. In this case, the constitutive laws of the membrane deformation are generally nonlinear but they are supposed homogeneous. In this study, we restrict our interest to flexible filaments in flows at low Reynolds number, where the hydrodynamics is linear [10]. Passive and active fibres [11–17], flagellum or cilia [18–21] have also been analysed in order to infer the internal parameters of their mechanical modelling. Gadêlha *et al.* [19] have used flagellum buckling experiments [20] to infer passive properties, including the passive sliding resistance, of a flagellum. Riedel-Kruse *et al.* [21] have fitted beating shapes with model prediction using a linearized one-dimensional wave propagation model for the beating propagation along an active flagellum in order to obtain either passive or active parameters. In the cited contributions, the shape is obtained from numerically solving a coupled fluid–structure problem, with mechanical properties being given. Internal parameter are then obtained from a systematic exploration of parameter space so as to reach a satisfying shape compared with experiments.

Whether such empirical parameter identification might be sufficient to provide interesting and insightful estimation of fibre or membrane mechanical properties could be a matter of debate. In many cases, a very precise estimation might not be so crucial and the mechanical properties might be approximately known from other experimental observations (e.g. osmotic swelling [5]).

Nevertheless, one key question which ought to be known is the uniqueness of the shape–motion/parameter relation. More specifically, are the shape and motion of a flexible object unique when varying internal parameters? This is obviously a crucial question as a non-unicity scenario would possibly lead to wrong estimation for the internal parameters. To our knowledge, in the case of nonlinear constitutive law for a flexible object (having large deformation) in a Stokes flow, this is an open question. This question is known as the identifiability condition for internal parameters. Even in the case of *linear* constitutive laws, but arbitrary deformation, identifiability conditions for deformable membrane or fibres are not known. The goal of this contribution is to clarify internal parameter identification in the simplest case of elastic flexible fibres in Stokes flow. Our aim is (i) to provide clear conditions for which the internal elastic properties of a passive flexible fibre or an active micro-swimmer can be identified from the observation of their kinematics within the framework of a bead-model (BM) description, (ii) to propose an operational formulation for reliable identification with and without additional noise, and (iii) to validate this operational formulation in various configurations in order to illustrate its interest. As the elasto-hydrodynamic fluid–structure interaction problem associated with (passive or active) flexible objects is *a priori* nonlinear, there is not much use of Stokes linearity for parameter identification. However, when a discrete description of deformable filaments by a collection of rigid bodies is considered to describe complex flexible objects namely BMs in Stokes flows [11, 16, 22] some simplifications arise. When considering a collection of non-overlapping beads we show that the internal forces linearly depend on passive and active forcing, producing a linear elasto-hydrodynamic problem. However, the obtained results are general for extensible BMs with suitable linearizations, or with pre-averaging, for which these models can be recast as Gauss–Markov processes and would provide equally valid likelihoods.

A detailed description of BMs is presented in §2a while focusing on kinematic constraints associated with the fibre inextensibility. More precisely §2b(ii) describes how contact forces are found from using Lagrangian multipliers, so that internal parameters and kinematic observations are linearly related. This point opens many issues that we subsequently analyse and illustrate for parameter identification in §2b. First, we provide clear and simple identifiability conditions associated with active or passive internal parameters showing that identifying concomitantly the active and passive internal properties of micro-swimmers is not possible. By contrast, the identification of mechanical properties assuming the active part is known, or reciprocally, finding the active part knowing the passive properties is possible, and can be addressed with the direct inversion of a simple linear problem without the need for guess values for the parameters. Section 3 analyses identifiability conditions and parameter estimation either in the deterministic case in §3a or the stochastic one in §3b. Identification in the stochastic case is possible using a Bayesian approach (as in [23]) and provides unbiased estimators the convergence of which can be estimated. Finally, §4 provides the numerical test and validation of the approach in various configuration associated with passive in §4a,b, actuated in §4c or active fibres in §4d.

### 2. Models and parameters

#### (a) Object discretization with the bead model

##### (i) Bead models.

To model the locomotion of a micro-swimmer from describing its structural, active and elastic properties, we hereby use a BM. This class of models is interesting as they permit describing a complex deformable object as a flexible assembly of simple rigid ones, such as ${N}_{b}$ bonded spheres. BMs have been used to provide a sensible description of an assembly of spheres connected together to model complex objects such as polymers, active fibres [11–14, 16, 17] or micro-swimmers. The connection between the spheres (beads or blobs) can be ensured in different ways. To account for the resistance to stretching and compression, linear springs are very often used [11, 24, 25]. Nevertheless, most flexible objects are very weakly extensible, the description of which results in stiff bonds distressfully restrictive for time integration. Joint models have been proposed to handle inextensible assembly of objects with kinematic constraints on the beads' motion. In the most formulation, the beads are constrained such that the contact point ${c}_{i}$ for each pair of beads remains the same [18, 26]. To allow the object to bend, a gap between each pair of beads is necessary. One problem arising in such a joint model is that in the case of strong bending, some overlapping between the objects can occur, so that a repulsive force is used to prevent it. Such repulsive forces are generally stiff and again very restrictive for time integration. An alternative joint model has been introduced by Yamamoto & Matsuoka [27, 28]. To avoid artificial gaps (and thus repulsive forces), they iteratively solved a nonlinear system to satisfy a no-slip condition between touching beads.

##### (ii) The gear model.

Recently, a linear formulation of kinematic constraints using Lagrange multiplier formalism has been proposed for BMs: the ‘gear bead model’ (GBM) [29]. In this model, there is neither gap nor repulsive force between the beads to allow bending and no need for numerical parameters to be tuned. Kinematic conditions are associated with rolling contact connections at each contact point between neighboured beads, providing constraints for the local velocity and angular momentum of each bead. For example, at the contact point ${c}_{i}$ between two beads *i* and $i+1$, the contact point velocity ${\mathbf{V}}_{{c}_{i}}$ has to match exactly with each bead rigid body velocity so that

*a*is the bead radius, ${\mathbf{t}}_{i}$ is the unit vector connecting their centre of mass and ${\mathbf{v}}_{i},{\mathit{\omega}}_{i}$ are the translation and rotational velocity of bead

*i*. Also, quantitative agreement with previous models has been obtained for both slender objects (passive fibre, actuated filaments) and non-slender swimmers (

*Caenorhabditis elegans*), allowing its use in a wide variety of contexts [29].

#### (b) Theoretical framework of parameter estimation

##### (i) Bead model formulation and context.

We consider the following quantities to describe flexible objects of complex shape; ${\mathbf{r}}_{i}$ and ${\mathbf{p}}_{i}$ denote the position and orientation of each bead $i\in \{1,{N}_{b}\}$ as in [29]. The translation velocities ${\mathbf{v}}_{i}$ are collected into the vector $\mathbf{V}\equiv [{\mathbf{v}}_{1},{\mathbf{v}}_{2},\dots {\mathbf{v}}_{{N}_{b}}]$, and the angular velocities ${\mathit{\omega}}_{i}$ into $\mathit{\Omega}\equiv [{\mathit{\omega}}_{1},{\mathit{\omega}}_{2},\dots {\mathit{\omega}}_{{N}_{b}}]$. Similarly, $\mathbf{F}\equiv [{\mathbf{f}}_{1},{\mathbf{f}}_{2},\dots {\mathbf{f}}_{{N}_{b}}]$ collects the forces ${\mathbf{f}}_{i}$ and $\mathit{T}\equiv [{\mathit{t}}_{1},{\mathit{t}}_{2},\dots ,{\mathit{t}}_{{N}_{b}}]$ the torques ${\mathit{t}}_{i}$, for which upper-indexes (‘*b*’ for bending or ‘*a*’ for active) are also added when needed for the sake of clarity so as to distinguish those torques from tangent vectors. Generalized velocities $\mathcal{V}\equiv [\left(\begin{array}{c}{\mathbf{v}}_{1}\\ {\mathit{\omega}}_{1}\end{array}\right),\hspace{0.17em}\left(\begin{array}{c}{\mathbf{v}}_{2}\\ {\mathit{\omega}}_{2}\end{array}\right)\dots \left(\begin{array}{c}{\mathbf{v}}_{{N}_{b}}\\ {\mathit{\omega}}_{{N}_{b}}\end{array}\right)]$ and forces $\mathfrak{F}\equiv [\left(\begin{array}{c}{\mathbf{f}}_{1}\\ {\mathit{t}}_{1}\end{array}\right),\hspace{0.17em}\left(\begin{array}{c}{\mathbf{f}}_{2}\\ {\mathit{t}}_{2}\end{array}\right)\dots \left(\begin{array}{c}{\mathbf{f}}_{{N}_{b}}\\ {\mathit{t}}_{{N}_{b}}\end{array}\right)]$ are also used in the following for the sake of compactness. We consider ${N}_{\mathrm{c}}$ linear non-holonomic kinematic constraints associated with generalized velocities acting on the ${N}_{b}$ beads of the form

*J*symbol), and might stand for a generic relation for linearized constraints. $\mathbf{J}$ is a ${N}_{\mathrm{c}}\times 6{N}_{b}$ sparse matrix and $\mathbb{B}$ is a vector of ${N}_{\mathrm{c}}$ components. It is worth mentioning that the rolling contact condition of the gear model [29] is a special case of (2.2) for which $\mathbb{B}=\mathbf{0}$ while $\mathbf{J}$ embeds all constraints (2.1). The equation (2.2) is linear because $\mathbf{J}$ and $\mathbb{B}$ might either depend on time

*t*, positions ${\mathbf{r}}_{i}$ and orientations ${\mathbf{p}}_{i}$ but not on $\mathcal{V}$. Vector $\mathbb{B}$ is thus useful to include the presence of additional effects located at some specific bead position of the fibre such as forced actuation at one end, as exemplified in [29]. This is why in the following we name $\mathbb{B}$ an actuation vector. To include constraints in the dynamics of articulated systems one should search for additional forces which permit satisfying these constraints (2.2), i.e. the so-called

*contact*forces. These

*contact*forces are given by

*λ*is the Lagrange multipliers vector. The entire framework presented in the paper is valid for both two- and three-dimensional dynamic configurations.

##### (ii) Contact forces and linearity of the kinematic constraints of bead model.

Hydrodynamic interactions are provided through the solution of the mobility problem which relates forces and torques to the velocity and angular velocity of the beads. This formalism is particularly suitable for low Reynolds number flows for which total forces and torques on each bead are zero (as the flow is inertia-less). The total generalized force, i.e. the sum of generalized hydrodynamic forces ${\mathfrak{F}}_{h}$ and non-hydrodynamic contributions $\mathfrak{F}$, is zero:

Hence, the contact forces can be found explicitly and depend linearly on other forces ${\mathfrak{F}}^{\mathrm{\prime}}$ and ambient flow

*effective*mobility matrix taking into account the action of contact forces associated with kinematic constraints. Also, note that if those contact forces (as would be the case for extensible fibres) are zero and kinematic constraints are also zero $\mathfrak{M}$ reduces to the standard mobility matrix. Hence the following results are also of general scope for active extensible polymer chains, for one to include extensible forces, as well as the necessary time-stepping constraints associated with stiff bonds.

In §4, we illustrate the formulation in the case of GBM associated with kinematic constraints with sliding contacts for which in most cases (except in §4c) actuation vector $\mathbb{B}=\mathbf{0}$ in (2.2), and furthermore with no external ambient flow, so that ${\mathcal{V}}^{\mathrm{\infty}}=\mathbf{0}$, and thus $\mathbb{K}=\mathbf{0}$.

In the case (most often encountered) where the elastic forces have linear dependence with constitutive parameters ${\mathit{\Theta}}_{\mathrm{e}}$, a vector composed of bending modulus for flexural deformation and shear modulus for torsion for each bead connection, we can write

Relation (2.13) is a major result of this paper as it gives an explicit, algebraic formulation for the linearity between kinematics and internal parameters. As kinematics is supposed to be experimentally known, this formulation allows for a direct determination of the identifiability conditions since, in the linear case, this is a straightforward issue.

### 3. Identifiability conditions

#### (a) Deterministic case

Let us suppose that generalized velocity $\mathcal{V}$(t) are the observable (e.g. which could be obtained from some particle tracking/image segmentation post-processing of experimental images). From (2.9), they depend nonlinearly on particle position and orientations $\mathcal{V}$(t) through the *effective* mobility matrix $\mathfrak{M}(\mathcal{V}(t))$. They also depend linearly on static unknown parameter ${\mathit{\Theta}}_{\mathrm{e}}$, and dynamic unknown active force ${\mathfrak{F}}_{\mathrm{a}}(t)$. Let us also suppose that ${\mathrm{{\rm Y}}}_{\mathrm{e}}$ in (2.12) is also known at each time step from linear constitutive properties of the deformable object. The inverse of matrix $\mathfrak{M}$ is an effective resistance matrix. In the following, two viewpoints are possible to address identifiablility. One possible approach is to consider kinematic observable $\mathcal{V}$ and try to match the model prediction in the kinematic variable. A supplementary approach is to consider ‘kinematically deduced observable forces’ from using (2.9)

##### (i) Identification of passive elastic parameters.

Let us start with the first case, suppose $\mathbf{y}(t)$ is known from the observation $\mathcal{V}(t)$ and solving (3.1). Let us suppose there is no active force ${\mathfrak{F}}_{\mathrm{a}}(t)=0$, but an external flow so that $\mathbb{K}\ne 0$ in (2.9) and (3.1). Then, at each time step, one is able to define the following functional

Equivalently, if one considers the kinematic observable, the minimization of

In both cases, a direct $6{N}_{b}\times 6{N}_{b}$ linear system has to be solved, which we never found poorly conditioned (but if the number of bead does not exceed 100, a direct inversion might be pursued). However, force formulation, from (3.1), requires the action of the effective ‘resistance’ matrix ${\mathfrak{M}}^{-1}$ onto the generalized velocity vector, which might be more costly (especially if one uses more precise mobility matrix than far-field approximation [35]).

##### (ii) Identification of active forces.

Identifying active forces knowing passive ones is trivial in force formulation, for which the functional minimization

##### (iii) Identifiability of both active force and passive elastic parameters.

Let us consider a first ‘naive’ approach for concomitantly identifying both active and passive forces over a single period of time. In this case, we can formulate parameter identification as the solution of the following minimization problem for the force formulation

*t*(a continuous formulation is also possible but adds unnecessary complexity), and where it is meaningful to suppose some periodicity in the active component, so that $t\in [0,T]$. In the following, we denote each field evaluated at discrete time step ${t}_{i}=iT/N$ with $i\in [1,N]$ with index

*i*, e.g. ${\mathfrak{F}}_{\mathrm{a}}({t}_{i})\equiv {\mathfrak{F}}_{\mathrm{a}}(i)$. The unknown parameters are thus the collection vector ${\mathfrak{F}}_{\mathrm{a}}(1),{\mathfrak{F}}_{\mathrm{a}}(2),\dots {\mathfrak{F}}_{\mathrm{a}}(N)$ and ${\mathit{\Theta}}_{\mathrm{e}}$. As the functional (3.10) is quadratic, its Jacobian ${J}_{\mathrm{F}}\equiv {\mathrm{\nabla}}_{{\mathfrak{F}}_{\mathrm{a}}(t),{\mathit{\Theta}}_{\mathrm{e}}}F$ provides the linear system to solve for, in order to find the unknown parameters, formally

*d*a continuous $d-1$ parameter family of possible sum can be found. In other words, this system is always undetermined. Similarly for the kinematic formulation

##### (iv) Identifiability of periodic active forcing.

However, if the active component is strictly periodic, instead of trying to identify at various discrete time step (phase) ${t}_{i}$ within one period, one could try to evaluate at one given phase, e.g. ${t}_{i}$, over several periods, e.g. ${t}_{i}+T$, ${t}_{i}+2T$, etc. For simplicity, let us hereby simply consider two distinct periods with two distinct evaluations, one at ${t}_{i}$, and the other one at ${t}_{i}+T$. As, ${\mathfrak{F}}_{\mathrm{a}}({t}_{i})={\mathfrak{F}}_{\mathrm{a}}({t}_{i}+T)$, it is easy to realize that

*T*-periodic. One can then take the difference between the first and the second line of (3.14) to eliminate the active force so as to identify passive properties using the ‘force formulation’ (3.3). To be explicit, denoting

*T*. Hence, even if this section has shown that both active and passive components are not identifiable at the same time, one can nevertheless find strategies to circumvent the indeterminacy in the special case of a periodic active force within a variable flow, not having the same period.

#### (b) Stochastic case

In the stochastic framework, we keep with the kinematic formulation as the noise is added onto the kinematics. Hence, let us consider the Langevin-like equations for the bead motion

*k*of vector $\mathcal{W}(n)$ follows ${\mathcal{W}}_{k}(n)\sim \mathcal{N}(0,{\sigma}^{2})$. The noise can either be attributed to some experimental uncertainty in the measurements, or some intrinsic random contributions of the background.

In a stochastic framework, $\mathcal{V}$ are the observables, bead positions and orientations are the variables and ${\mathfrak{F}}_{\mathrm{a}}(t)$ and/or ${\mathit{\Theta}}_{\mathrm{e}}$ are the parameters to be estimated. In this framework, (3.17) describes a normal linear regression model (e.g. [36]). In the next section, we use a classical Bayesian approach to estimate parameters from observations having i.i.d. Gaussian residuals. It is worth mentioning the presented approach could also handle more complex and general forcing where, for example, the components of ${\mathcal{W}}_{k}(n)$ are not independent anymore (e.g. follow a joint Gaussian distribution) as can result from a Brownian forcing of the flow field, which produces correlated bead motion. We do not provide such analysis here, but this is a possible extension of the proposed approach.

In the following, as we found in the deterministic case that ${\mathfrak{F}}_{\mathrm{a}}(t)$ and ${\mathit{\Theta}}_{\mathrm{e}}$ are not simultaneously identifiable, we will mainly focus in the identification and estimation convergence of either one or the other.

##### (i) Identification of passive elastic parameters.

In this subsection, we consider (3.17) with ${\mathfrak{F}}^{\mathrm{\prime}}={\mathbf{{\rm Y}}}_{\mathrm{e}}{\mathit{\Theta}}_{\mathrm{e}}$ and we use the classical Bayesian approach for which the conditional probability $P({\mathit{\Theta}}_{\mathrm{e}}|\mathcal{V})$ for parameters ${\mathit{\Theta}}_{\mathrm{e}}$ given observations $\mathcal{V}$ reads

*N*independent observations are

*N*increases.

##### (ii) Identification of active forces.

In this subsection, we suppose that passive elastic properties are known and estimated. We consider (3.17) with ${\mathfrak{F}}^{\mathrm{\prime}}={\mathfrak{F}}_{\mathrm{a}}(t)+{\mathbf{{\rm Y}}}_{\mathrm{e}}(t){\hat{\mathit{\Theta}}}_{\mathrm{e}}$. In the following, we suppose that ${\mathfrak{F}}_{\mathrm{a}}(t)$ is a periodic function of *t* with period *T*. As we are now searching for a time-periodic parameter vector, it is not wise anymore to consider *any* discrete set of time ${t}_{n}$ to sample the problem. Thus, here we specifically consider that each time *t* is sampled *N* times along period *T*, at time ${t}_{n}=t+nT$ with $n\in [1,N]$. We consider the system to follow a closed periodic orbit without noise, so that $\mathcal{V}(t)$ and ${\mathfrak{F}}_{\mathrm{a}}(t)$ are two *T*-periodic functions. In practice, *T* can be found from analysing the time variations of $\mathcal{V}(t)$, for example, using the frequency associated with its Fourier transform highest peak. If passive properties are known, a very precise knowledge of the period *T* is not critical as it will only affect the evaluation of ${\mathfrak{F}}_{\mathrm{a}}(t)$ nearby value $t=T$. However, if one wishes to estimate beforehand these passive properties in the presence of active forces using (3.16), a precise estimate of *T* would be mandatory, which might depend on the number of period sampling of
$\mathcal{V}(t)$.

In the following, we thus keep with a continuous time representation of ${\mathfrak{F}}_{\mathrm{a}}(t)$ for notation purposes, but, obviously for numerical implementation, this continuous time can also be discretized later on. Defining ${\mathbb{K}}^{\mathrm{\prime}}=\mathbb{K}+{\mathbf{{\rm Y}}}_{\mathrm{e}}{\hat{\mathit{\Theta}}}_{\mathrm{e}}$, we now find the likelihood $\mathcal{L}({\mathfrak{F}}_{\mathrm{a}}\hspace{0.17em}|\hspace{0.17em}\mathcal{V})$ as

### 4. Validation tests

In this section, we illustrate through various examples the identification results obtained in the previous section. For conciseness, we focus our numerical tests to kinematic formulation. Nevertheless, in most cases we numerically evaluate both formulations without notable difference. Here, we consider a quiescent fluid ($\mathbb{K}=\mathbf{0}$). We provide for each case the numerical parameters used for readers to be able to reproduce our test cases. As our motivation is mainly a computational proof of the method, we did not used dimensional quantities, so that bead radius, viscosity are set to unity. We also give the relevant dimensionless numbers (sperm number, elasto-gravitational number) for the chosen test cases. Let us first consider the continuous description of an inextensible object which experiences bending torques and elastic forces. The Serret–Frenet frame triad $(\mathbf{t},\mathbf{n},\mathbf{b})$ needs the definition of tangent $\mathbf{t}(s)$

*s*of the fibre, the normal vector $\mathbf{n}(s)$

*i*th component of a bending rigidity elastic vector ${\mathit{\theta}}_{\mathrm{e}}$ associated with the bending rigidity of the link between bead

*i*and bead $i+1$ already defined in (2.12). The associated bending torque for each bead

*i*is given by

Finally, it is important to mention that for these numerical illustrations, we use a Rotne–Prager mobility matrix (cf. [29] for more details) and a third-order Adam–Bashforth scheme for time integration.

#### (a) Elastic relaxation of passive fibres without noise

In this example, we consider elastic relaxation towards equilibrium of a single filament discretized into ${N}_{b}-1$ rigid rods of length $2a$ without external force. The fibre is initialized with a given radius of curvature and restores its equilibrium straight shape along a transient controlled by elastic stresses and hydrodynamic interactions (figure 1). Table 1 provides the parameters used for this numerical computation where we have used a variable bending stiffness along the fibre.

parameter | value | parameter | value |
---|---|---|---|

time step | ${10}^{-1}$ | filament length | $L=40$ |

initial curvature radius | $R{u}_{\mathrm{ini}}=L/2$ | intrinsic curvature | $R{u}_{\mathrm{eq}}={10}^{8}\hspace{0.17em}L$ |

max bending stiffness | ${\mathit{\Theta}}_{\mathrm{em}}=4.71\times {10}^{3}$ | sperm number | $2\mu {a}^{4}/{\mathit{\Theta}}_{\mathrm{em}}=2/{\mathit{\Theta}}_{\mathrm{em}}$ |

Our aim is to estimate the bending rigidity ${\theta}_{\mathrm{e}}({s}_{i})$ on each bead *i*, i.e. the vector ${\mathit{\Theta}}_{\mathrm{e}}$, by using the identification formula (3.5). We choose ${\theta}_{\mathrm{e}}$ to vary along the line *s* to show that our method works in the general case of bodies with heterogeneous properties. We found in §3(i) that the vector ${\mathit{\Theta}}_{\mathrm{e}}$ was identifiable as long as the matrix $\mathfrak{M}\mathbf{{\rm Y}}$ has full column rank. This approach provides an exact (up to computer accuracy) estimation of ${\mathit{\Theta}}_{\mathrm{e}}$. We use positions and velocities given by a simulation of the model choosing an initial configuration and relaxing the fibre to its equilibrium position. Computing the ${L}_{2}$ relative error between the estimated parameter ${\mathit{\Theta}}_{\mathrm{e}}$ and the exact one ${\mathit{\Theta}}_{\mathrm{e}}^{\mathrm{E}}$ at each time *t*, $\parallel {\mathit{\Theta}}_{\mathrm{e}}^{\mathrm{E}}-{\mathit{\Theta}}_{\mathrm{e}}{\parallel}_{{L}_{2}}/\parallel {\mathit{\Theta}}_{\mathrm{e}}^{\mathrm{E}}{\parallel}_{{L}_{2}}$, we found a deviation of about ${10}^{-14}$, close to computer accuracy. It is interesting to mention that, in this passive case, choosing non-uniform elastic properties along the fibre, as, for example, done in [21] for bull sperm flagellum, does not reduce the precision of the estimation. This was not the case in [21] when both active and passive parameter estimation were concomitantly determined. This might be due to the fact that, as shown in §3, simultaneous active and passive identification is not unique.

#### (b) Settling passive fibres without noise

Here we estimate the bending stiffness of a flexible filament settling under gravity force ${\mathbf{F}}_{\perp}={F}_{\perp}{\mathbf{e}}_{\perp}$ acting perpendicularly to its major axis. The fibre is initialized with a straight shape orthogonal to gravity. The filament reaches a steady shape and settles at a constant speed resulting both from the balance of elastic stresses, gravitational acceleration and hydrodynamic interactions [29, 39]. Adding the following generalized external force ${\mathfrak{F}}_{\perp}={({\mathbf{F}}_{\perp}\hspace{0.17em}\mathbf{0}\hspace{0.17em}\dots \hspace{0.17em}{\mathbf{F}}_{\perp}\hspace{0.17em}\mathbf{0})}^{\mathrm{T}}$, to the right-hand side of (2.13), so that, without active force, the relation between generalized velocity and forces simply reads

parameter | value | parameter | value |
---|---|---|---|

time step | ${10}^{-4}$ | filament length | $L=40$ |

intrinsic curvature | $R{u}_{\mathrm{eq}}={10}^{7}L$ | initial curvature | $R{u}_{\mathrm{ini}}=R{u}_{\mathrm{eq}}$ |

equilibrium bending angle | ${\theta}_{{b}_{\mathrm{eq}}}=0$ | gravity force | ${F}_{\perp}=10$ |

elasto-gravitational number | $B={F}_{\perp}{L}^{3}/{\mathit{\Theta}}_{\mathrm{e}}=3\times {10}^{3}$ | bending stiffness | $\parallel {\mathit{\Theta}}_{\mathrm{e}}\parallel =213.3$ |

#### (c) Actuated passive filaments without noise

In order to test if the proposed formulation also permits identification for actuated filaments, we propagate a bending from one end along the passive elastic filament, similarly to [40–43]. To simulate the experiment described in [42, 43] (attached filament actuated by its base), we use actuated prescribed three-dimensional motion of the first two beads, so-called tethered base elements [29]. This actuation requires the addition of three vectorial kinematic constraints to the no-slip conditions, adding ${\mathbf{J}}^{\mathrm{act}}$ and ${\mathbb{M}}^{\mathrm{act}}$ to $\mathbf{J}$ and $\mathbb{B}$ in the constraints equation, as detailed in [29]. Table 3 provides the parameters used for this numerical computation where we consider a homogeneous bending stiffness over the entire fibre (figure 3).

parameter | value | parameter | value |
---|---|---|---|

time step | ${10}^{-3}$ | filament length | $L=20$ |

intrinsic curvature radius | $R{u}_{\mathrm{eq}}={10}^{8}L$ | initial curvature | $R{u}_{\mathrm{ini}}=R{u}_{\mathrm{eq}}$ |

bending stiffness | $\parallel {\mathit{\Theta}}_{\mathrm{e}}\parallel ={10}^{4}$ | number of periods | 5 |

frequency | $f=\pi $ | sperm number, ${S}_{p}$ | $6\pi f\hspace{0.17em}{L}^{4}/{\mathit{\Theta}}_{\mathrm{e}}=94.74$ |

Here, we just write the matrix formulation of the model and the parameter estimation. We note $\mathfrak{J}=[{\mathbf{J}}^{\mathrm{act}}\hspace{0.17em}\mathbf{J}]$ a vertical concatenation and $\mathbb{B}={\mathbb{B}}^{\mathrm{act}}$ in (2.10) reads

#### (d) Active filaments without noise

Locomotion of the nematode *C. elegans* is used here as in [44]. Here we are just interested in modelling the motion of the nematode using the framework of the BM. To do so, we use an oscillating driving torque ${\mathbf{t}}^{a}(s,t)$ to mimic the internal muscular contractions and a curvature model given by [44]. The active torque applied on bead *i* results from the difference in active bending moments across neighbouring links

parameter | value | parameter | value |
---|---|---|---|

number of spheres | 20 | filament length | $L=40$ |

initial radius of curvature | $R{u}_{\mathrm{ini}}=L/2$ | intrinsic curvature | $R{u}_{\mathrm{eq}}={10}^{8}\hspace{0.17em}L$ |

number of periods | 2 | frequency | $f=2$ |

wavenumber | $k=3\pi /2L$ | amplitude | ${\theta}_{0}^{a}=8.25/L$ |

sperm number | ${S}_{p}=L{(f\mu /{\mathit{\Theta}}_{\mathrm{e}})}^{1/4}=3$ | time step | ${10}^{-3}$ |

Table 5 provides the results for the active amplitude wave at a given time step. The results are obtained from the direct inversion of relation (3.9), while solving ${\mathfrak{F}}_{\mathrm{a}}=\mathfrak{M}{\mathrm{{\rm Y}}}_{a}{\mathit{\Theta}}_{a}$, at a given time-step (at all time-step). Again, computer accuracy is achieved if one knows exactly the elastic properties.

link number i |
exact value ${\mathit{\Theta}}^{{a}_{\mathrm{e}}}\times {10}^{2}$ | ${\mathit{\Theta}}^{a}\times {10}^{2}$ (${\mathit{\Theta}}_{\mathrm{e}}$ known) | relative error $(\times {10}^{-14})$ |
---|---|---|---|

1 | −1.318155924588519 | −1.318155924588518 | 0.043123440710212 |

2 | −0.445553256609255 | −0.445553256609253 | 0.350843360582397 |

3 | 0.524174207564081 | 0.524174207564081 | 0.027110938520310 |

4 | 1.379638534110982 | 1.379638534110980 | 0.144205863415543 |

5 | 1.740923696038537 | 1.740923696038537 | 0.032651298267784 |

6 | 1.447190937009256 | 1.447190937009254 | 0.098196128456760 |

7 | 0.874900978326801 | 0.874900978326801 | 0.064971259912771 |

8 | 0.315226288359876 | 0.315226288359876 | 0.078892518390509 |

#### (e) Elastic relaxation with noise

In this section, we test the estimator prediction in the stochastic case. We explore the elastic relaxation already investigated in §4a in the presence of noise. The main purpose of this section is to test the theoretical predictions (3.23) and test the convergence toward this prediction given by the deviation law (3.26). The deviation (3.26) provides a normal convergence with a rather complex covariance matrix which non-trivially depends on the number of independent observations *N*. Thus we compute the statistical quadratic error, when for any noisy configuration we estimate ${\mathit{\Theta}}_{\mathrm{e}}(i)$ for $i=1,N$

*N*. The relative error over

*N*independent observations collected during the fibre relaxation have been averaged over 10 different statistical samples in order to compare results to the expected theoretical error (3.26). Both are in very good agreement, as shown in figure 4. In the presented case of large noise amplitude $\sigma =1$, one can see that almost a thousand estimations are necessary to reach a 1% accuracy in the prediction. Obviously, the error is proportional to the noise amplitude, and thus can be reduced for lower noise amplitude.

This computation illustrates and provides consistent results with the prediction obtained in previous section for noisy measurements.

### 5. Conclusion and perspectives

We have achieved internal parameter identification of active or passive inextensible fibres experiencing elasto-hydrodynamic coupling with Stokes flows. The fibre dynamics is described by an assembly of beads. Using Stokes flow linearity and explicit computation of contact forces associated with kinematic constraints, we derived a linear relationship between kinematics and internal parameters. The relation involves an explicit combination of kinematic constraints (through the Jacobian matrix and/or the actuation vector), hydrodynamic coupling (through Mobility matrix), elastic behaviour (through internal constitutive matrix) and external flow. This linear relation has enabled us to provide identifiability conditions. We find that either active or passive parameter identification is possible but not both simultaneously. Nevertheless, in the case of a periodic active forcing, we propose two possible approaches to circumvent indeterminacy and provide successively both passive and active properties.

We also found that in each case (passive or active), the presence of noise does not lead to a ill-posed inversion, but a well-conditioned direct inversion problem. We derive theoretical convergent estimators of the parameters. We tested this convergence numerically for various configurations and found consistent results with theoretical predictions. In each case of passive, actuated or active fibre we provide a direct estimate for the internal parameters without need of an initial guess, at the very moderate cost of a small linear system inversion. For noisy data, we provide error estimations of the internal parameter versus evaluation number which are fully consistent with numerical test-cases.

These results have been derived in the case of linear dependence of constitutive laws with parameters, but in the general case of arbitrary deformations (not necessarily small) and possibly nonlinear dependence of constitutive laws with shape (e.g. elastic bending depends on local curvature which is a nonlinear function of shape). The hypothesis of constitutive laws linearity with parameters might be seen as a limitation of the proposed results. As a perspective, let us briefly touch on how part of the proposed approach might be useful to tackle nonlinear parameter estimation. Obviously, in the general case of nonlinear dependence of constitutive laws with parameters, most of the analytic explicit estimators will not apply anymore. Nevertheless, the BM approach might still be useful, if one is able to derive a discretization of these continuous constitutive laws. The main point is that part of the nonlinearity associated with the fluid–structure interaction issues can still be tackled by the contact force model. Indeed, the linear relation between general velocities and forces (2.9) is still valid for forces having nonlinear parameter dependence. From then, the inversion problems that we derived might be useful in various iterative schemes associated with nonlinear inversion problems (e.g. Newton methods). Hence, our derivation might also been seen as a first useful step to more complex parameter estimations.

To provide more comments on the physical and biological relevance of our result it is important to stress that we found that two independent evaluations of active and passive property might provide a complete inversion of the internal parameters of flagellum or cilia. Such independent evaluations are interesting because they can both be used under experimental conditions, provided a precise quantitative observation of the fibre kinematics. First, non-invasive, ‘on the fly’ observations of passive fibre deformations should permit to evaluate internal mechanical properties.

This approach should also appear fruitful to provide internal parameter estimation of active micro-swimmers in the future, if, by any mean, one is able to temporarily switch-off the active swimming forces (e.g. using ion-pumps inhibitors). Passive properties could first be evaluated while later on the active swimming force components could be estimated from switching-on active effects.

### Data accessibility

All parameters associated with the simulations are given in the tables.

### Author's contributions

F.P. has provided the theoretical background and writing of §§2 and 3 as well as appendix A. E.I.T. has performed the programming, testing and validation of §4 results, and also contributed to §3 and appendix A. B.D. and E.C. have helped in the programming of section 4 and careful reading, corrections, comments and modifications of the manuscript.

### Competing interests

We declare we have no competing interests.

### Funding

This work has been supported by the Agence Nationale pour la Recherche (ANR) under grant no. MOTIMO (ANR-11-MONU-009-01).

## Acknowledgments

The authors thanks Jérôme Fehrenbach, Alexandre Steyer and Florencio Balboa-Usabiaga for their enlightening discussions and comments.

## Appendix A

**(a) The internal constitutive matrix ${\mathrm{{\rm Y}}}_{\mathrm{e}}$ in bending elastic relaxation**

In the case of elastic relaxation, the bending moment reads

*i*, again, ${\mathbf{b}}_{i}$ is the local bi-normal of the local Serret–Frenet base at the centre of bead

*i*, and $\kappa ({s}_{i})$ is the local curvature of the bead. ${\theta}_{\mathrm{e}}({s}_{i})$ is the

*i*th component of bending rigidity elastic vector ${\mathit{\theta}}_{\mathrm{e}}$. In most cases studied in this paper the rigidity ${\theta}_{\mathrm{e}}({s}_{i})\equiv {\theta}_{\mathrm{e}}$ is chosen constant for all beads, but we keep a general formulation for which it could spatially differ. $\kappa ({s}_{i})$ is given by (4.3), consistently with [45], so that the intermediate vector $\mathit{\xi}({s}_{i})$ is useful to consider

*i*with $i=2,\dots ,{N}_{b}-1$. The bending torque derived from the moment of each bead reads

Although the generalized elastic force, when the collection of the force ${\mathbf{f}}_{i}$ was zero for all beads, is expressed as $\mathfrak{F}=({\mathfrak{f}}_{1},\dots ,{\mathfrak{f}}_{{N}_{b}})$ with ${\mathfrak{f}}_{i}=(0,{\mathbf{t}}_{i})$ we have to rearrange the previous matrix ${\mathrm{{\rm Y}}}_{\mathrm{e}}$ to obtain the ${\mathrm{{\rm Y}}}_{\mathrm{e}}$ matrix elastic parameter. Finally, we get

**(b) The internal constitutive matrix ${\mathrm{{\rm Y}}}_{\mathrm{a}}$ in swimming nematode**

Now, in the case of swimming nematodes, the active bending moment components are written so an active, time-periodic curvature wave is prescribed on each bead

The active bending torque is given by

### References

- 1
Park HM, Choi YJ . 2009 Investigation of structural parameters of dilute polymer solutions using velocity measurements.**J. Non-Newton. Fluid Mech.**, 29–34. (doi:10.1016/j.jnnfm.2009.08.001) Crossref, Web of Science, Google Scholar**164** - 2
Mader MM, Ez-Zahraouy H, Misbah C, Podgorski T . 2007 On coupling between the orientation and the shape of a vesicle under a shear flow.**Eur. Phys. J.**, 275–280. (doi:10.1140/epje/e2007-00029-6) Web of Science, Google Scholar**22** - 3
Biben T, Farutin A, Misbah C . 2011 Three-dimensional vesicles under shear flow: numerical study of dynamics and phase diagram.**Phys. Rev.E**, 031921. (doi:10.1103/PhysRevE.83.031921) Crossref, Web of Science, Google Scholar**767** - 4
Coupier G, Farutin A, Minetti C, Podgorski T, Misbah C . 2012 Shape diagram of vesicles in poiseuille flow.**Phys. Rev. Lett.**, 178106. (doi:10.1103/PhysRevLett.108.178106) Crossref, PubMed, Web of Science, Google Scholar**108** - 5
Risso F, Collé-Paillot F, Zagzoule M . 2006 Experimental investigation of a bioartificial capsule flowing in a narrow tube.**J. Fluid Mech.**, 149–173. (doi:10.1017/S0022112005007652) Crossref, Web of Science, Google Scholar**547** - 6
de Loubens C, Deschamps J, Edwards-Levy F, Leonetti M . 2016 Tank-treading of microcapsules in shear flow.**J. Fluid Mech.**, 750–767. (doi:10.1017/jfm.2015.758) Crossref, Web of Science, Google Scholar**789** - 7
Carin M, Barthès-Biesel D, Edward-Lévy F, Postel C, Andrei DC . 2003 Compression of biocompatible liquid-filled hsa-alginate capsules: determination of the membrane mechanical properties.**Biotechnol. Bioengng.**, 207–212. (doi:10.1002/bit.10559) Crossref, PubMed, Web of Science, Google Scholar**82** - 8
Lefebvre Y, Leclerc E, Barthès-Biesel D, Walter J, Edwards-Lévy F . 2003 Flow of artificial microcapsules in microfluidic channels: a method for determining the elastic properties of the membrane.**Phys. Fluids**, 123102. (doi:10.1063/1.3054128) Crossref, Web of Science, Google Scholar**20** - 9
de Loubens C, Deschamps J, Boedec G, Leonetti M . 2015 Stretching of capsules in an elongation flow a route to constitutive law.**J. Fluid Mech.**, R3. (doi:10.1017/jfm.2015.69) Crossref, Web of Science, Google Scholar**767** - 10
Pozrikidis C . 1992**Boundary integral and singularity methods for linearized viscous flow**. Cambridge, UK: Cambridge University Press. Crossref, Google Scholar - 11
Jayaraman G, Ramachandran S, Ghose S, Laskar A, Saad Bhamla M, Sunil Kumar PB, Adhikari R . 2012 Autonomous motility of active filaments due to spontaneous flow-symmetry breaking.**Phys. Rev. Lett.**, 158302. (doi:10.1103/PhysRevLett.109.158302) Crossref, PubMed, Web of Science, Google Scholar**109** - 12
Chelakkot R, Gopinath A, Mahadevan L, Hagan MF . 2013 Flagellar dynamics of a connected chain of active, polar, brownian particles.**J. R. Soc. Interface**, 20130884. (doi:10.1098/rsif.2013.0884) Link, Web of Science, Google Scholar**11** - 13
Laskar A, Singh R, Ghose S, Jayaraman G, Kumar PBS, Adhikari R . 2013 Hydrodynamic instabilities provide a generic route to spontaneous biomimetic oscillations in chemomechanically active filaments.**Sci. Rep.**, 1964. (doi:10.1038/srep01964) Crossref, PubMed, Web of Science, Google Scholar**3** - 14
Ghosh A, Gov NS . 2014 Dynamics of active semiflexible polymers.**Biophys. J.**, 1065–1073. (doi:10.1016/j.bpj.2014.07.034) Crossref, PubMed, Web of Science, Google Scholar**107** - 15
Eisenstecken T, Gompper G, Winkler RG . 2016 Conformational properties of active semiflexible polymers.**Polymers**, 304. (doi:10.3390/polym8080304) Crossref, Web of Science, Google Scholar**18** - 16
Laskar A, Adhikari R . 2016 Brownian microhydrodynamics of active filaments.**Soft Matter**, 9073–9085. (doi:10.1039/C5SM02021B) Crossref, Web of Science, Google Scholar**11** - 17
Winkler RG . 2016 Dynamics of flexible active Brownian dumbbells in the absence and the presence of shear flow.**Soft Matter**, 3737–3749. (doi:10.1039/C5SM02965A) Crossref, PubMed, Web of Science, Google Scholar**12** - 18
Lindstrom SB, Uesaka T . 2007 Simulation of the motion of flexible fibers in viscous fluid flow.**J. Fluid Mech.**, 113307. (doi:10.1063/1.2778937) Google Scholar**19** - 19
Gadêlha H, Gaffney EA, Goriely A . 2013 The counterbend phenomenon in flagellar axonemes and cross-linked filament bundles.**Proc. Natl Acad. Sci. USA**, 12 180–12 185. (doi:10.1073/pnas.1302113110) Crossref, Web of Science, Google Scholar**110** - 20
Dominic DW, Pelle W, Brokaw CJ, Lesich CB, Lindemann KA . 2009 Mechanical properties of the passive sea urchin sperm flagellum.**Cell Motil. Cytoskeleton**, 721–735. (doi:10.1002/cm.20401) Crossref, PubMed, Web of Science, Google Scholar**66** - 21
Riedel-Kruse IH, Hilfinger A, Howard J, Jülicher F . 2007 How molecular motors shape the flagellar beat.**HFSP J.**, 192–208. (doi:10.2976/1.2773861) Crossref, PubMed, Google Scholar**1** - 22
Swan JW, Brady JF, Moore RS . 2011 Modeling hydrodynamic self-propulsion with stokesian dynamics. Or teaching stokesian dynamics to swim.**Phys. Fluid**, 071901. (doi:10.1063/1.3594790) Crossref, Web of Science, Google Scholar**23** - 23
Bouffanais R, Sun J, Yue DKP . 2013 Physical limits on cellular directional mechanosensing.**Phys. Rev. E**, 052716. (doi:10.1103/PhysRevE.87.052716) Crossref, Web of Science, Google Scholar**87** - 24
Gauger E, Stark H . 2006 Numerical study of a microscopic artificial swimmer.**Phys. Rev. E**, 021907. (doi:10.1103/PhysRevE.74.021907) Crossref, Web of Science, Google Scholar**74** - 25
Manghi M, Schlagberger X, Netz RR . 2006 Propulsion with a rotating elastic nanorod.**Phys. Rev. Lett.**, 068101. (doi:10.1103/PhysRevLett.96.068101) Crossref, PubMed, Web of Science, Google Scholar**96** - 26
Skjetne P, Ross RF, Klingenberg DJ . 1997 Simulation of single fiber dynamics.**J. Chem. Phys.**, 2108. (doi:10.1063/1.474561) Crossref, Web of Science, Google Scholar**107** - 27
Yamamoto S, Matsuoka T . 1993 A method for dynamic simulation of rigid and flexible fibers in a flow field.**J. Chem. Phys.**, 644. (doi:10.1063/1.464607) Crossref, Web of Science, Google Scholar**98** - 28
Yamamoto S, Matsuoka T . 1995 Dynamic simulation of fiber suspensions in shear flow.**J. Chem. Phys.**, 2254–2260. (doi:10.1063/1.468746) Crossref, Web of Science, Google Scholar**102** - 29
Delmotte B, Climent E, Plouraboué F . 2015 A general formulation of bead models applied to flexible fibers and active filaments at low Reynolds number.**J. Comput. Phys.**, 14–37. (doi:10.1016/j.jcp.2015.01.026) Crossref, Web of Science, Google Scholar**286** - 30
Wajnryb E, Mizerski KA, Zuk PJ, Szymczak P . 2013 Generalization of the Rotne-Prager-Yamakawa mobility and shear disturbance tensors.**J. Fluid Mech.**, R3. (doi:10.1017/jfm.2013.402) Crossref, Web of Science, Google Scholar**731** - 31
Brady JF, Bossis G . 1988 Stokesian dynamics.**Annu. Rev. Fluid Mech.**, 111–157. (doi:10.1146/annurev.fl.20.010188.000551) Crossref, Web of Science, Google Scholar**20** - 32
Maxey MR, Patel BK . 2001 Localized force representations for particles sedimenting in stokes flow.**Int. J. Multiph. Flow**, 1603–1626. (doi:10.1016/S0301-9322(01)00014-3) Crossref, Web of Science, Google Scholar**27** - 33
Yeo K, Maxey MR . 2010 Simulation of concentrated suspensions using the force-coupling method.**J. Comput. Phys.**, 2401–2421. (doi:10.1016/j.jcp.2009.11.041) Crossref, Web of Science, Google Scholar**229** - 34
Roper M, Dreyfus R, Baudry J, Fermigier M, Bibette J, Stone HA . 2006 On the dynamics of magnetically driven elastic filaments.**J. Fluid Mech.**, 167. (doi:10.1017/S0022112006009049) Crossref, Web of Science, Google Scholar**554** - 35
Balboa Usabiaga F, Kallemov B, Delmotte B, Bhalla A, Griffith BE, Donev A . 2016 Hydrodynamics of suspensions of passive and active rigid particles: a rigid multiblob approach. (http://arxiv.org/abs/1602.02170) Google Scholar - 36
- 37
Bishop TC, Cortez R, Zhmudsky OO . 2004 Investigation of bend and shear waves in a geometrically exact elastic rod model.**J. Comput. Phys.**, 642–665. (doi:10.1016/j.jcp.2003.08.028) Crossref, Web of Science, Google Scholar**193** - 38
Landau LD, Lifshitz EM . 1975**Theory of elasticity**, 3rd edn. Oxford, UK: Pergamon Press. Google Scholar - 39
Lagomarsino MC, Capuani F, Lowe CP . 2003 A simulation study of the dynamics of a driven filament in an Aristotelian fluid.**J. Theor. Biol.**, 215–224. (doi:10.1016/S0022-5193(03)00159-0) Crossref, PubMed, Web of Science, Google Scholar**224** - 40
Wiggins CH, Riveline D, Ott A, Goldstein RE . 1998 Trapping and wiggling: elastohydrodynamics of driven microfilaments.**Biophys. J.**, 1043–1060. (doi:10.1016/S0006-3495(98)74029-9) Crossref, PubMed, Web of Science, Google Scholar**74** - 41
Yu TS, Lauga E, Hosoi AE . 2006 Experimental investigations of elastic tail propulsion at low Reynolds number.**Phys. Fluids**, 091701. (doi:10.1063/1.2349585) Crossref, Web of Science, Google Scholar**18** - 42
Coq N, du Roure O, Marthelot J, Bartolo D, Fermigier M . 2008 Rotational dynamics of a soft filament: wrapping transition and propulsive forces.**Phys. Fluids**, 051703. (doi:10.1063/1.2909603) Crossref, Web of Science, Google Scholar**20** - 43
Coq N, DuRoure O, Fermigier M, Bartolo D . 2011 The shape of an elastic filament in a two-dimensional corner flow.**J. Phys. Condens. Matter**, 063602. (doi:10.1063/1.3601446) Google Scholar**23** - 44
Majmudar T, Keaveny EE, Zhang J, Shelley MJ . 2012 Experiments and theory of undulatory locomotion in a simple structured medium.**J. R. Soc. Interface**, 1809–1823. (doi:10.1098/rsif.2011.0856) Link, Web of Science, Google Scholar**9** - 45
LJ Fauci, CS Peskin . 1988 A computational model of aquatic animal locomotion.**J. Comput. Phys.**, 85–108. (doi:10.1016/0021-9991(88)90158-1) Google Scholar**108**