Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences

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 [79] 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 [1117], flagellum or cilia [1821] 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 Nb 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 [1114, 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 ci 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 ci between two beads i and i+1, the contact point velocity Vci has to match exactly with each bead rigid body velocity so that

Vci=vi+ati×ωi=vi+1ati×ωi+1,2.1
where, a is the bead radius, ti is the unit vector connecting their centre of mass and vi,ω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; ri and pi denote the position and orientation of each bead i{1,Nb} as in [29]. The translation velocities vi are collected into the vector V[v1,v2,vNb], and the angular velocities ωi into Ω[ω1,ω2,ωNb]. Similarly, F[f1,f2,fNb] collects the forces fi and T[t1,t2,,tNb] the torques ti, 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 V[(v1ω1),(v2ω2)(vNbωNb)] and forces F[(f1t1),(f2t2)(fNbtNb)] are also used in the following for the sake of compactness. We consider Nc linear non-holonomic kinematic constraints associated with generalized velocities acting on the Nb beads of the form

JV+B=0,2.2
where J is the matrix associated with the constraints. As these constraints are linear here, J also represents the constraints' Jacobian (justifying J symbol), and might stand for a generic relation for linearized constraints. J is a Nc×6Nb sparse matrix and B is a vector of Nc components. It is worth mentioning that the rolling contact condition of the gear model [29] is a special case of (2.2) for which B=0 while J embeds all constraints (2.1). The equation (2.2) is linear because J and B might either depend on time t, positions ri and orientations pi but not on V. Vector 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 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
Fc=λTJ,2.3
where λ 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 Fh and non-hydrodynamic contributions F, is zero:

Fh+F=0.2.4
As a result, the generalized velocities V are linearly related to the generalized external forces F through the mobility matrix M:
V=MF+V,2.5
where V[(v1ω1),(v2ω2)(vNbωNb)] corresponds to the generalized ambient external velocity v and vorticity ω evaluated at the centre of mass of each bead. Note that if the particles are immersed in a linear shear flow, the response to the ambient rate of strain can be added without loss of generality [30]. This problem is nonlinear on the instantaneous positions of all particles of the system. Although being an explicit and linear relation, (2.5) is non-trivial as the exact mobility matrix results from the solution of a many-body problem. Approximated mobility matrices have nevertheless been proposed to circumvent this difficulty. The simplest model (free drain) assume a diagonal mobility matrix neglecting hydrodynamic interactions between neighboured spheres. More elaborated models include far field interactions using explicit analytic expressions known as Rotne–Prager–Yamakawa tensor [30]. Finally, other formulations compute numerically, at each time step, the complete many-body problem either using Stokesian dynamics (SD) [22, 31] or the force coupling method (FCM) [32, 33]. In the following, we do not focus on how the mobility matrix should be computed nor discuss the numerical issues associated with its calculation. We suppose this is a known quantity, and most of the theoretical results are obtained independently from its precise evaluation. When needed in the following sections, we will use the classic explicit formulation of Rotne–Prager–Yamakawa tensor to illustrate parameter identifications methodology. Non-hydrodynamic forces F are the sum of elastic forces Fe, inner forces associated with active mechanisms Fa and contact forces Fc. Let us note F=Fe+Fa, so that F=F+Fc. Relation (2.5) stands for the elasto-hydrodynamic coupling of coupled objects, and, provides a linear relationship between active forces, elastic forces and velocities. Nevertheless, as internal contact forces Fc are unknown and functions of all other quantities, the linearity is not obvious at this stage. The inextensibility condition provides the supplementary condition to set in order to get internal forces governing equations (e.g. [34]). In the framework of BM, instead of the inextensibility condition, we propose prescribing kinematic constraints [29] which then impose the inextensibility condition. The main point is then to find the internal forces needed for kinematic constraints to hold at every time. Using (2.5) and (2.2), we are able to find the Lagrangian multipliers for internal forces
MFc=VMFV,2.6
so that applying J to (2.6) provides the linear system for Lagrange multipliers estimation
JMJTλ=BJ(MF+V),2.7
while we have used the prescribed constraints (2.2) on the first term of the right-hand side of (2.7). As J is full column rank, and M is invertible, JMJT is also invertible (it is also simple to show that JMJT is symmetric). In the case of non-actuated, inextensible fibres, the actuation vector B vanishes, B=0. Nevertheless, for actuated fibres driven from one end there is a non-zero contribution associated with actuated beads which drive the motion of the others as exemplified in [29] as well as in §4c.

Hence, the contact forces can be found explicitly and depend linearly on other forces F and ambient flow

Fc=(JMJT)1{BJ(MF+V)}J.2.8
Using the explicit contact force expression (2.8) in the mobility matrix expression (2.5) leads to the following generalized velocities
V=MF+K,2.9
while using notations
M=M(I6NbJT(JMJT)1JM)2.10
and
K=MJT(JMJT)1B(MJT(JMJT)1JI6Nb)V,2.11
where I6Nb is the 6Nb×6Nb identity matrix. This expression provides an explicit linear formulation of the generalized velocities with external, internal and elastic forces, for a flexible object having a complex shape, as the internal contact forces have been explicitly expressed as linear functions of kinematic and external forces themselves in (2.8). Note that, M(V(t)) is configuration dependent, i.e. depends on bead position and orientation at each time step, through the constraints Jacobian J and mobility matrix M dependence. Furthermore, K(V(t),V) is also configuration dependent and external flow dependent. Using symbol M in (2.10) the definition is intentional as a reminder that this matrix acts as an 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 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 B=0 in (2.2), and furthermore with no external ambient flow, so that V=0, and thus K=0.

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

Fe=ΥeΘe,2.12
where the internal constitutive matrix Υe will be exemplified in §4 as well as explicitly given in appendix A. To avoid unnecessary complications in the following, we restrict our discussion to bending and will not cover the effect of torsion. We would nevertheless like to stress that there is no particular limitation in adding torsion in the forthcoming discussions. This matrix provides the structure of the (possibly complex) constitutive linear behaviour of the fibre so as to describe its twist and bending properties. For the GBM, (2.9) proves that kinematics depends linearly on these parameters.
VK=M(Fa+ΥeΘe).2.13

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 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 V(t) through the effective mobility matrix M(V(t)). They also depend linearly on static unknown parameter Θe, and dynamic unknown active force Fa(t). Let us also suppose that Υe in (2.12) is also known at each time step from linear constitutive properties of the deformable object. The inverse of matrix M is an effective resistance matrix. In the following, two viewpoints are possible to address identifiablility. One possible approach is to consider kinematic observable 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)

y(t)=M1(VK)3.1
and compare predicted forces with the one deduced from the kinematics. Both viewpoints are considered in the following. More precisely, we show how a simple quadratic error (the L2 norm of the difference between observable and model prediction in continuous time) minimization permits one to obtain an explicit direct linear system to either evaluate the internal properties in a given flow in the absence of active component, or, unknown active forces, when elastic parameter of the object are known (in the absence of any flow).

(i) Identification of passive elastic parameters.

Let us start with the first case, suppose y(t) is known from the observation V(t) and solving (3.1). Let us suppose there is no active force Fa(t)=0, but an external flow so that K0 in (2.9) and (3.1). Then, at each time step, one is able to define the following functional

F(Θe)=(y(t)Υe(t)Θe)2,3.2
the minimization of which is related to the explicit computation of its Jacobian, leading to the following linear system
Υe(t)TΥe(t)Θe=Υe(t)Ty(t).3.3
This least-squares problem has a unique solution if Υe(t) has full column rank, which is shown in appendix A.

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

F(Θe)=(V(t)KMΥe(t)Θe)2,3.4
leads to
((MΥe)T(MΥe))(t)Θe=(MΥe)T(t)(VK),3.5
which has a unique solution if MΥe(t) has full column rank. Both the force formulation (3.3) and the kinematic formulation (3.5) are interesting to consider for the evaluation of the elastic properties of passives fibres.

In both cases, a direct 6Nb×6Nb 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 M1 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

F(Fa(t))=(y(t)Υe(t)ΘeFa(t))23.6
leads to the direct result
Fa(t)=y(t)Υe(t)Θe.3.7
Similarly, in kinematic formulation, the functional reads
F(Fa(t))=(V(t)KM(Υe(t)Θe+Fa(t)))2,3.8
leading to
Fa(t)=(MTM)1(MT(V(t)K))Υe(t)Θe,3.9
which necessitates the inversion of MTM at each time step to find the active force Fa(t).

(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

F(Fa(t),Θe)=t(y(t)Fa(t)Υe(t)Θe)2,3.10
where we have implicitly discretized the time 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[0,T]. In the following, we denote each field evaluated at discrete time step ti=iT/N with i[1,N] with index i, e.g. Fa(ti)Fa(i). The unknown parameters are thus the collection vector Fa(1),Fa(2),Fa(N) and Θe. As the functional (3.10) is quadratic, its Jacobian JFFa(t),ΘeF provides the linear system to solve for, in order to find the unknown parameters, formally
JF[ΘeFa]=[ΣiΥe(i)Ty(i)y(t)].3.11
A formal identifiability condition follows from (3.11): detJF0. Computing the determinant of JF, we explicitly find that
detJF=iΥeT(i)Υe(i)+iΥeT(i)Υe(i)=0.3.12
Thus, we find that kinematic observations are not sufficient to provide enough information to distinguish both time-dependent active force and static elastic parameters. They are not identifiable simultaneously. This is true even if one is able to process an arbitrary number of measurements, with the hope that increasing the number of observations could compensate for parameter ignorance. Such compensation just does not hold here. Let us pause a little to examine this issue in the simplest case. Suppose we have only one observation at one time step, and suppose that we only want to identify both active force Fa and elastic one Fe. As the functional (3.10) depends only on the sum of those, there is already an infinite combination of Fa and Fe giving the same Fa+Fe for living the functional invariant. More precisely, in two dimensions there is an infinite continuous, one-parameter family of the possible sum of two vectors leading to a given vector: the ones built from any cord in the circle having the sum as its diameter. In higher dimensions d a continuous d1 parameter family of possible sum can be found. In other words, this system is always undetermined. Similarly for the kinematic formulation
F(Fa(t),Θe)=t(V(t)KM(Υe(t)Θe+Fa(t)))2,3.13
we find that the determinant of the Jacobian JF is also zero. It is easy to understand that whichever, kinematic or forces formulation, both active and passive forces cannot be simultaneously identified because, both functional (3.10) and (3.13) depend on the sum of active and passive contributions, none of which can be identified from the action of both on the kinematics. This simple statement is of general scope, and results from the linear underlying relation between kinematics and both active and elastic forces described in (2.13).

(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) ti within one period, one could try to evaluate at one given phase, e.g. ti, over several periods, e.g. ti+T, ti+2T, etc. For simplicity, let us hereby simply consider two distinct periods with two distinct evaluations, one at ti, and the other one at ti+T. As, Fa(ti)=Fa(ti+T), it is easy to realize that

M1(ti+T)(V(ti+T)K)=Υe(ti+T)Θe+Fa(ti)andM1(ti)(V(ti)K)=Υe(ti)Θe+Fa(ti).}3.14
Let us suppose that (V(ti+T)K)(V(ti)K), because K(ti+T)K(ti) as, for example, the external flow velocity also varies with time and is not 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
Δy=M(ti+T)1(V(ti+T)K)M(ti)1(V(ti)K)3.15
and ΔΥe=Υe(ti+T)Υe(ti), the elastic parameters can be found from solving
Θe=(ΔΥeTΔΥe)1ΔΥeTΔy.3.16
Once the passive properties are known, one could also use (3.9) to evaluate the active force. It is interesting to mention that this estimation of passive properties necessitates the precise knowledge of the active forcing period 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

V(t)=M(t)F(t)+K(t)+W(t),3.17
where W(t) is a Wiener process vector for which each component is an identically independently distributed (i.i.d.) Wiener process. For each discrete time step, using same notation as in the previous section, we denote W(tn)W(n) and then each i.i.d. component k of vector W(n) follows Wk(n)N(0,σ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, V are the observables, bead positions and orientations are the variables and Fa(t) and/or Θ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 Wk(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 Fa(t) and Θ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 F=ΥeΘe and we use the classical Bayesian approach for which the conditional probability P(Θe|V) for parameters Θe given observations V reads

P(Θe|V)=(P(V|Θe)P(Θe)P(V))L(Θe|V)P(Θe),3.18
depending on prior distribution of parameters P(Θe). The likelihood L(Θe|V) can be found from (3.17) using the fact that each component of W(n) is i.i.d. and Gaussian, from a 6Nb components vector (as V), the likelihood of N independent observations are
L(Θe|V)=1N1(2π)3Nbσ6Nbe((VKMΥeΘe)T(VKMΥeΘe))/2σ2.3.19
In the following, we will estimate the maximum of (3.18) in the special case of uniform (constant) P(Θe) as no prior knowledge on parameter values is supposed. If needed, it is possible to include supplementary terms associated with proper prior distribution of parameters P(Θe), as an extension of the proposed approach. Hence, equivalently, the maximum of (3.18) is also the maximum of the logarithm of L(Θe|V), the so-called log-likelihood
lnL(Θe|V)=3NNbln2π3NNblnσ21N12σ2(VKMΥeΘe)T(VKMΥeΘe).3.20
It is worth noting that, up to constants, the log-likelihood displays a quadratic functional form similar to (3.4) or (3.13). Again, if needed the prior distribution of parameters P(Θe) would add a supplementary term to (3.20), e.g. also a quadratic one, for a Gaussian prior. The computation of the maximum of this log-likelihood essentially follows the same footsteps. The gradient of lnL in parameters directions being zero leads to the following linear system, the noise amplitude now being part of the unknown estimate
[lnLΘelnLσ2]=[1σ21N(MΥe)T(n)(VKMΥeΘe)(n)3NNbσ2+1N(VKMΥeΘe)T(VKMΥeΘe)2σ4]=[00].3.21
It leads to the following linear system inversion for the estimation of the parameter Θ^e
1N(MΥe)T(MΥe)(n)Θ^e=1N(MΥe)T(n)(VK)(n),3.22
so that
Θ^e=(MΥe)TMΥe1(MΥe)T(VK),3.23
with 1N/N is the time average of each matrix and vector. (3.23) is thus very similar to (3.5) derived in the case K=0, except for the sum over different time associated with some noise-averaging convergence of the stochastic framework. The classical result for the variance also holds from the second line of (3.21)
σ^2=16Nb(VKMΥeΘ^e)T(VKMΥΘ^e).3.24
In this case, the convergence towards the estimated parameters Θ^e and σ^ is provided by the Cramér–Rao lower bound associated with the eigenvalue of the expectation of the inverse of the likelihood Hessian. In the case of Gaussian noise, the result is explicit, as the Hessian expectation reads
E[2lnLΘe22lnLΘeσ22lnL[σ2]22lnLΘeσ2]=[1N(MΥ)TMΥσ2003NbNσ4]3.25
and thus the inverse of the Hessian expectation provides the covariance matrix associated with the convergence toward the estimators, i.e.
Θ^eΘeN(0,σ2NMΥ(MΥ)T1MΥ(MΥ)TMΥ1(MΥ)T),3.26
which shows that the covariance of the fluctuation to the estimator vector Θ^e is governed by a covariance matrix directly related to the matrix needed for computing Θ^e (3.23). More importantly, a central-limit theorem holds for the convergence toward the estimator Θ^e as the number of considered evaluation 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 F=Fa(t)+Υe(t)Θ^e. In the following, we suppose that Fa(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 tn to sample the problem. Thus, here we specifically consider that each time t is sampled N times along period T, at time tn=t+nT with n[1,N]. We consider the system to follow a closed periodic orbit without noise, so that V(t) and Fa(t) are two T-periodic functions. In practice, T can be found from analysing the time variations of 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 Fa(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 V(t).

In the following, we thus keep with a continuous time representation of Fa(t) for notation purposes, but, obviously for numerical implementation, this continuous time can also be discretized later on. Defining K=K+ΥeΘ^e, we now find the likelihood L(Fa|V) as

L(Fa|V)=1N1(2π)3Nbσ6Nbe((VKMFa)T(VKMFa))/2σ2.3.27
The footsteps of the previous section are almost identical leading to the linear system for the estimation of vector parameter F^a
1NMTM(n)F^a=1NMT(n)(VK)(n),3.28
so that
F^a(t)=MTM(t)1MT(VK)(t),3.29
where to be precise here, (t)=(1/N)n=1N(t+nT). Furthermore, following the same footsteps as in the previous section, it is easy to find the deviation distribution to F^a
F^aFaN(0,σ2NMMT1MMTM1MT).3.30

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 (K=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 (t,n,b) needs the definition of tangent t(s)

t(s)=r(s)s1r/s,4.1
where r(s) is the continuous vector position varying along the centreline arclength s of the fibre, the normal vector n(s)
n(s)=ts1t/s,4.2
so that the binormal vector is given by b(s)=t(s)×n(s). The local curvature κ(s) of the fibre is also given by
κ(s)=ts(s)n(s)=ts(s).4.3
The bending moment is given by the constitutive law of elasticity [37, 38]
m(s)=θe(s)t(s)×ts(s)=θe(s)κ(s)b(s),4.4
where θe(s) is the bending rigidity, which can vary along the fibre centreline, and, again, κ(s) the local curvature. Let us now describe the discrete implementation of the bending moment (4.4) (which is a bending moment with no intrinsic curvature or zero equilibrium curvature)
m(si)=θe(si)κ(si)bi.4.5
θe(si) is the ith component of a bending rigidity elastic vector θ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
tib=m(si+1)m(si1).4.6
Furthermore, using first and last bending moments along the fibre as
t1b=m(s2),t2b=m(s3),tNb1b=m(sNb2)andtNbb=m(sNb1)
permits explicitly building the linear response matrix Υe defined in (2.12). Appendix A provides an explicit formulation of how this matrix is built.

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 Nb1 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.

Figure 1.

Figure 1. Various snapshots of the elastic relaxation of fibre of N=20 beads towards equilibrium along time (grey levels). (xc,yc) corresponds to the centre of mass of the fibre and a the radius of the beads. In this case, without noise, measuring the generalized velocity V permits the exact evaluation of elastic properties along the fibre Θe using (3.5).

Table 1.Parameters of elastic relaxation of passive fibre simulation. Number of beads N=20.

parameter value parameter value
time step 101 filament length L=40
initial curvature radius Ruini=L/2 intrinsic curvature Rueq=108L
max bending stiffness Θem=4.71×103 sperm number 2μa4/Θem=2/Θem

Our aim is to estimate the bending rigidity θe(si) on each bead i, i.e. the vector Θe, by using the identification formula (3.5). We choose θ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 Θe was identifiable as long as the matrix MΥ has full column rank. This approach provides an exact (up to computer accuracy) estimation of Θ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 L2 relative error between the estimated parameter Θe and the exact one ΘeE at each time t, ΘeEΘeL2/ΘeEL2, we found a deviation of about 1014, 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 F=Fe 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 F=(F0F0)T, to the right-hand side of (2.13), so that, without active force, the relation between generalized velocity and forces simply reads

V=M(F+ΥΘe),4.7
so that one can use relation (3.5) from replacing V by the known quantity VMF (as M is known when the bead positions are given), and thus estimate Θe. Table 2 provides the parameters used for this numerical computation where we consider a homogeneous bending stiffness over the entire fibre (figure 2). The obtained estimation of Θe is correct, up to computer accuracy, when the fibre gets its equilibrium position at the final time. The L2 relative error between the estimated parameter Θe and the exact one ΘeE at the final time is 109. Finally, note that one needs to know the passive settling force F beforehand when performing this identification.
Figure 2.

Figure 2. Various snapshot of the bending relaxation of fibre composed of N=20 beads settling under gravity in an infinite fluid domain along time (grey levels). The model used in this example has been validated in [29]. (xc,yc) corresponds to the centre of mass of the fibre and a the radius of the beads. In this case without noise, the estimation of the general velocity V permits the evaluation of elastic properties using (3.5).

Table 2.Parameters for the settling passive fibre identification. We use N=20 beads for the evaluation.

parameter value parameter value
time step 104 filament length L=40
intrinsic curvature Rueq=107L initial curvature Ruini=Rueq
equilibrium bending angle θbeq=0 gravity force F=10
elasto-gravitational number B=FL3/Θe=3×103 bending stiffness Θe=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 [4043]. 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 Jact and Mact to J and 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).

Figure 3.

Figure 3. Illustration of the GBM for an actuated filament with N=10 beads tethered at the base. (a) Snapshots of the filament centreline at 10 times along the actuation cycle (grey levels): the fibre displays an hellicoidal tri-dimensional shape, the form of which depends on internal elastic properties and (b) discretization with the BM.

Table 3.Parameters of the actuated passive filament identification. Number of beads N=10.

parameter value parameter value
time step 103 filament length L=20
intrinsic curvature radius Rueq=108L initial curvature Ruini=Rueq
bending stiffness Θe=104 number of periods 5
frequency f=π sperm number, Sp 6πfL4/Θe=94.74

Here, we just write the matrix formulation of the model and the parameter estimation. We note J=[JactJ] a vertical concatenation and B=Bact in (2.10) reads

Mact=M(I6NbJT(JMJT)1JM),4.8
while (2.11) becomes
K=MJT(JMJT)1Bact4.9
so that (2.9) writes
V=MactΥΘe+K4.10
and the estimation linear system is
Θe=((MactΥ)T(MactΥ))1(MactΥ)T(VexpK).4.11
Again our direct estimate provides the correct estimate for the elastic bending properties up to computer accuracy. The L2 relative error between the estimated parameter Θe and the exact one ΘeE is about 1011. Again, note here that one needs to know the actuation Bact to evaluate K in (4.9) beforehand when performing this identification.

(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 ta(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

ti,ta=ma(si+1,t)ma(si1,t),4.12
with
ma(si,t)=θ0a(si)sin(ksi2πft)b(si)=θa(si)b(si),4.13
using notation θa=θ0asin(ksi2πft). The total bending results from the contributions of the active moment and the passive one. Firstly, we can write the relation between velocities and forces without external forces to identify these parameters according to consistent equations. Writing again the general velocity as
V=M(Fe+Fa)=MΥe(t)Θe+MΥa(t)Θa(t),
where we have introduced a known matrix Υa(t) which describes a ‘wave curvature’ model previously proposed in [44] (details are given in appendix A) and unknown active amplitudes Θa associated with the wave propagation along the fibre. Supposing that the elastic internal parameters Θe are known, let us now discuss the results found for Θa estimation. Table 4 provides the parameters used for this numerical computation.

Table 4.Parameters of the active filament identification. Number of beads N=20.

parameter value parameter value
number of spheres 20 filament length L=40
initial radius of curvature Ruini=L/2 intrinsic curvature Rueq=108L
number of periods 2 frequency f=2
wavenumber k=3π/2L amplitude θ0a=8.25/L
sperm number Sp=L(fμ/Θe)1/4=3 time step 103

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 Fa=MΥaΘa, at a given time-step (at all time-step). Again, computer accuracy is achieved if one knows exactly the elastic properties.

Table 5.Estimation of Θa(t) at a given time t in the swimming nematode case with 10 beads.

link number i exact value Θae×102 Θa×102 (Θe known) relative error (×1014)
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 Θe(i) for i=1,N

eN=ΘeEΘeL2NΘeEL2,4.14
for various numbers of independent observations 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 σ=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.
Figure 4.

Figure 4. Error (4.14) to prediction (3.23) for the elastic relaxation versus the number of observations. The dotted line represents the prediction given by the theoretical estimate (3.26), whereas black bullets are the results obtained with numerical stochastic average in the presence of O(1) noise amplitude σ=1. (Online version in colour.)

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 Υe in bending elastic relaxation

In the case of elastic relaxation, the bending moment reads

m(si)=θe(si)κ(si)bi,A 1
where si denotes the curvilinear position of bead i, again, bi is the local bi-normal of the local Serret–Frenet base at the centre of bead i, and κ(si) is the local curvature of the bead. θe(si) is the ith component of bending rigidity elastic vector θe. In most cases studied in this paper the rigidity θe(si)θe is chosen constant for all beads, but we keep a general formulation for which it could spatially differ. κ(si) is given by (4.3), consistently with [45], so that the intermediate vector ξ(si) is useful to consider
ξ(si)=κ(si)bi=12ati1×ti,A 2
for each bead, i with i=2,,Nb1. The bending torque derived from the moment of each bead reads
tib=m(si+1)m(si1)A 3
and
tib=κ(si+1)bi+1θe(si+1)κ(si1)bi1θe(si1).A 4
Then, using notation (A 2)
tib=ξ(si+1)θe(si+1)ξ(si1)θe(si1).A 5
To be explicit, the general bending torque γb for six beads reads
(t1bt2bt3bt4bt5bt6b)=(ξ(s2)0000ξ(s3)00ξ(s2)0ξ(s4)00ξ(s3)0ξ(s5)00ξ(s4)0000ξ(s5))(θe(s2)θe(s3)θe(s4)θe(s5)),A 6
where obviously in (A 6) all tib and ξ(si) are three-component vectors and 0=(000).

Although the generalized elastic force, when the collection of the force fi was zero for all beads, is expressed as F=(f1,,fNb) with fi=(0,ti) we have to rearrange the previous matrix Υe to obtain the Υe matrix elastic parameter. Finally, we get

Fe=(0t1b0t2b0t3b0t4b0t5b0t6b)=(0000ξ(s2)00000000ξ(s3)000000ξ(s2)0ξ(s4)000000ξ(s3)0ξ(s5)000000ξ(s4)00000000ξ(s5))Υe(θe(s2)θe(s3)θe(s4)θe(s5))Θe.A 7

(b) The internal constitutive matrix Υ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

ma(si,t)=θa(si)sin(ksi2πft)b(si).A 8
Now, using notation Θa(si)=θa(si)sin(ksi2πft), we simplify (A 8)
ma(si,t)=Θa(si,t)b(si).A 9

The active bending torque is given by

tia=ma(si+1)ma(si1)=Θa(si+1,t)b(si+1)Θa(si1,t)ba(si1).
For six beads the total force associated with the active torque t=ta is
Fa=(0t1a0t2a0t3a0t4a0t5a0t6b)=(0000b(s2)00000000b(s3)000000b(s2)0b(s4)000000b(s3)0b(s5)000000b(s4)00000000b(s5))Υa(Θa(s2,t)Θa(s3,t)Θa(s4,t)Θa(s5,t))Θa,A 10
where again 0 is a column vector with three zeros.

Footnotes

Published by the Royal Society. All rights reserved.

References