A monolithic optimal control method for displacement tracking of Cosserat rod with application to reconstruction of C. elegans locomotion

This article considers an inverse problem for a Cosserat rod where we are given only the position of the centreline of the rod and must solve for external forces and torques as well as the orientation of the cross sections of the centreline. We formulate the inverse problem as an optimal control problem using the position of the centreline as an objective function with the external force and torque as control variables, with meaningful regularisation of the orientations. A monolithic, implicit numerical scheme is proposed in the sense that primal and adjoint equations are solved in a fully-coupled manner and all the nonlinear coefficients of the governing partial differential equations are updated to the current state variables. The forward formulation, determining rod configuration from external forces and torques, is first validated by a numerical benchmark; the solvability and stability of the inverse problem are then tested using data from forward simulations. The proposed optimal control method is motivated by reconstruction of the orientations of a rod’s cross sections, with its centreline being captured through imaging protocols. As a case study, we take the locomotion of the nematode, Caenorhabditis elegans. In this study we take laboratory data for its centreline and infer its cross-section orientation (muscle locations) with the control force and torque being interpreted as the reaction force, activated by C. elegans’ muscles, from the surrounding fluids. This method thus combines the mathematical modelling and laboratory data to study the locomotion of C. elegans, which gives us insights into the potential anatomical orientation of the worm beyond what can be observed through the laboratory data. The paper is completed with several additional remarks explaining the theoretical and technical details of the model.


Introduction
The beam and rod theories have been developed to model a typical three dimensional solid structure which is much longer in one dimension than the other two dimensions. The classical Euler-Bernoulli beam theory considers the extension and compression of a rod or beam and allows loads of stretching, compressing or bending [6,24], which is generally suitable for modelling a thin rod with small deformation. The Timoshenko-Ehrenfest beam theory was developed to take into account shear deformation induced by rotational bending effects, making it suitable for modelling thick beams with larger deformation [19,23]. The Kirchhoff-Love and Cosserat rod theories were developed to model rods with finite deformation, with the former allowing bending and twisting while ignoring stretching, compressing and shearing deformation [21,56]; and the latter allowing all types of loads and deformation [2,18,29,74,75]. This paper is to for-mulate an optimal control problem based upon the special Cosserat theory of rods. Cosserat rod theory is geometrically exact for modelling bending and torsion as well as extension and shear, which is considered as a geometrically nonlinear generalisation of the Timoshenko-Ehrenfest beam (while Kirchhoff-Love is a geometrically nonlineaer generalisation of Euler-Bernoulli beam) [2,74], and has been adopted to model the locomotion of Caenorhabditis elegans (C. elegans) in recent years [31,68].
In the context of optimal control of a rod or beam, the solution existence of optimal control of the longitudinal vibration of a viscoelastic rod by either a contact force or distributed force is discussed in [73], and the mean mechanical energy minimised by a boundary force is studied, using the methods of the calculus of variations [30], maximum principle [70] and Ritz method [51]; minimisation of the mean square deviation of the Timoshenko beam is investigated by controlling a distributed force [71] or by the angular acceleration [83], and singularity of its solution is discussed in [69]; optimal control of transverse vibration of Euler-Bernoulli beam is introduced in [79]. We will consider displacement tracking of the Cosserat rod in this paper, which, to the best of our knowledge, has not been studied before. Due to their long, thin nature, accurately capturing the rotation vectors of a rod is a challenging task whereas there are several well studied approaches for reconstructing the centreline [8,27,32]. To investigate the applicability of our method we will consider a case study on the nematode, C. elegans: We shall apply this optimal control formulation to the reconstruction of the locomotion of C. elegans based upon laboratory data [72].
Locomotion -the ability of an organism to move from one place to another -is achieved by animals through a variety of methods [9,16,35]. C. elegans is a transpar-ent nematode of about 1 mm long [11,76], whose planar undulatory locomotion has been widely studied [17,48] by laboratory experiments [8,25,47,54] or mathematical modelling [31,63,68]. In its natural habitat, this nematode moves in three dimensional environments. However, such locomotion has only recently begun to be recorded [52,72]. Of particular interest is a recent modelling study of a roll manoeuvre modelled as a torsional turn [10]. One key challenge with interpretation of 3D video footage of locomotion is that the images (and hence the centreline reconstructions) lack information linking the body shape to the local, anatomically meaningful frame of the body (its left, right, ventral and dorsal directions (see Figure 4), and information about internal torsion or twist along the body).
In this paper, we propose a method to combine laboratory data of the motion of C. elegans' centreline and mathematical modelling to reconstruct the whole picture of C. elegans locomotion: how does the worm wriggle and wiggle locally through its body (how does its anatomical frame evolve)?
The centreline data of C. elegans can be constructed using videos from three different perspectives [72]. However, it is challenging to construct the local frames (information about the internal torsion or twist along the body), which is the motivation to develop the proposed method in this paper. We point out that the proposed optimal control formulation is general and not limited to C. elegans but for simplicity we restrict the example formulation to neglect inertial terms [17] (consistent with C. elegans being a low Reynolds number swimmer).
The contributions of this paper are highlighted as follows: a monolithic optimal control method is developed based upon the special Cosserat theory of rods; an incompressibility condition is derived and integrated into the forward as well as the control problem; the formulation is implicit and the primal and adjoint equations are solved in a fully-coupled manner; this new optimal control method is applied to a challenging inverse problem: reconstruction of C. elegans locomotion based on its centreline from laboratory data; implementation in the open-source software package FreeFEM++, which is available on the public Github site.
The paper is organised as follows. The governing partial differential equations of the Cosserat rod are introduced in Section 2 with a focus on expressing these control equations in a closed component form. The optimisation problem with the corresponding primal and adjoint equations are derived in Section 3, followed by a monolithic optimal control formulation in Section 4. Numerical experiments are carried out in Section 5 to validate both the forward and the optimal control formulations, and the proposed optimal control method is applied to the reconstruction of C. elegans locomotion in Section 6. Finally, conclusions are drawn and future work are discussed in Section 7. First, two coordinate systems or frames, as well as their relations, are introduced in order to describe the geometry of the Cosserat rod. Then, the mechanics of Cosserat rod is described by the conservation of linear momentum and angular momentum, and a set of constitutive equations is introduced to close the system. Finally, the governing equations are expressed in terms of six unknown variables: three components of the position vector (x, y, z) and three components of the rotation vector (α, β, γ ). This formulation is based on the one presented in [12]: we derive the angular velocity and generalised curvature using a new method in Section 2.3 and rewrite all the control equations in a matrixvector format; in addition, we consider dilation of the cross section of the rod by differentiation of the reference and current arc lengths and derivation of the incompressibility condition in Section 2.4.

Global and local frames
In order to describe all the types of deformation of the Cosserat rod, a global coordinate system [e 1 , e 2 , e 3 ] is first introduced as shown in Figure 1, which is assumed to form a fixed right-hand orthogonal unit basis (also called the fixed frame), to define the centreline of the rod by a threedimensional curve: r(s, t) = x(s, t)e 1 + y(s, t)e 2 +z(s, t)e 3 , where s ∈ [a(t), b(t)] is the arc-length parameter of the curve and t denotes the time; a local coordinate system (the moving frame) [d 1 (s, t), d 2 (s, t), d 3 (s, t)] (orthogonal unit basis) is also introduced everywhere at the centreline to describe the motion of the rod's cross section, and it is assumed that d 3 (s, t) is always perpendicular (not necessarily coinciding with the tangential ∂ s r(s, t) of the centreline) to the cross section to facilitate the expressions of the moment of inertia and constitutive relations. This local coordinate system can be constructed by the following three successive rotations from the global coordinate system.

Remark 1
The arc length s is a current configuration, which is generally not a constant especially when considering the case of large extension or compression [31]. We introduce a reference or initial configurations = s 0 ∈ [a 0 , b 0 ] to compute the strain (equations (9) and (10)), and consider s = s(s, t) as a function ofs and time t. Let us also introduce the deformation scalar j(s, t) = ds(s, t)/ds for the convenience of notation in the following sections.
Step 1: Rotate the e 1 − e 3 plane clockwise around the e 2 axis by an angle γ , so that the e 3 axis sits on the d 2 − d 3 plane as shown in Figure 2 (left), i.e.: perform a rotation operation [e 1 , e 2 , e 3 ] R T y (γ ) with Step 2: Rotate the e 2 − e 3 plane clockwise around the e 1 axis by an angle β, so that the e 3 axis overlaps with the d 3 axis as shown in Figure 2 (middle), i.e.: perform another rotation operation [e 1 , e 2 , Step 3: Rotate the e 1 − e 2 plane clockwise around e 3 axis by an angle α, so that the [e 1 , e 2 , e 3 ] overlaps with [d 1 , d 2 , d 3 ] as shown in Figure 2 (right), i.e.: perform the final rotation operation [e 1 , e 2 , The overall rotation matrix can be expressed as: cos α cos γ − sin α sin β sin γ − sin α cos β cos α sin γ + sin α sin β cos γ sin α cos γ + cos α sin β sin γ cos α cos β sin α sin γ − cos α sin β cos γ − cos β sin γ sin β cos β cos γ where all the three angles α, β and γ are functions of the arc length s and time t: α = α(s, t), β = β(s, t) and γ = γ (s, t). Therefore, the local coordinate system can be obtained by The components of any vector v in these two coordinates system have the following relations: if v is expanded in the global frame as which implies In the rest of this article, we use the superscript 'g' to indicate the components of a vector expanded in the global frame and 'l' in the local frame.

Conservation laws
The governing equations of the Cosserat rod are based on the conservation of linear momentum and conservation of angular momentum as follows [2]: (6) where n and m are the internal force and torque respectively, f and l are the external force and torque densities (per unit reference length) respectively, ρ(s) and A(s, t) are the density and area of the cross section respectively, and h is the angular momentum (per unit reference length).
In equations (5) and (6), it is convenient to express all the vectors in the local frame except r(s, t). Therefore, we first express these vectors in the local frame, and then transform them to the global frame using (4), which will finally be substituted into (5) and (6) in order to obtain an equation system in its component form.
The angular momentum h = 3 i=1 h l i d i can be expressed as: with ω = 3 i=1 ω l i d i denoting the generalised angular velocity, and I(s) denoting the moment of inertia (per unit reference length). Let (ξ, η, ζ ) denote the coordinates in the local frame, then I(s) can be computed as follows, noticing that d 3 is perpendicular to the cross section: In order to close the equation system (5) and (6), constitutive relations of n and m have to be established in terms of the unknown variables. In the local frame, we adopt a linear relation between the internal force n(s, t) and strain (s, t), and linear relation between internal torque m and the curvature κ(s, t) [12,31] as follows. Let and a linear relation, in the local frame, between n(s, t) and (s, t) can be expressed as where E and G are the Young's and shear moduli respectively, k is a numerical factor depending on the shape of the cross section at s [34], and A(s, t) is the area of the rod's cross section.
We assume A(s, t) is a function of space s and time t, and an incompressibility assumption will be used to determine A(s, t) in Section 2.4. Using the transformation (4) between the local and global coordinates, (11) can be expressed as ⎛ ⎝ n g 1 Similarly, let Then, a linear relation, in the local frame, between m(s, t) and κ(s, t) can be expressed as where Eξ 2 dξ dη, Remark 2 For the main context of this paper, we consider a circular cross section (with the exception of a rectangular cross section that is used in numerical test 5.1 for validation against a published result), and constant density ρ, Young's modulus E and shear modulus G. In which case, In the spirit of expressing all the unknown variables in terms of (x, y, z) and (α, β, γ ), we further express the angular velocity ω and curvature κ in terms of rotation angles (α, β, γ ) in the following section.

Expressions of angular velocity and curvature in terms of the angles of rotation
For any fixed-length vector function, say, of t, with constant c. This suggests that ∂ t v is always perpendicular to v. If vector v rotates according to an angular velocity ω as shown in Figure 3, we have Since d 1 , d 2 and d 3 are all unit vectors, we can apply the above property (17) to these three vectors and have Following the same argument, we also have: Now, let Q T = q 1 , q 2 , q 3 with q T i = (q i1 , q i2 , q i3 ), i = 1, 2, 3, being the row vectors of Q, then from (2) we have and Using the fact that for any two vectors we can compute the cross product in (18): Finally using (18) and (23) to (25), the angular velocity ω, in the local frame, can be expressed as: noticing that for i = j (i, j = 1, 2, 3) A further calculation based on (1) and (26) expresses rotation angles in the local frame as follows: Using the same procedure, the curvature at a point along the centreline can be expressed in the local frame as: Substituting equation (13) into equation (5), we express the conservation of linear momentum in its component form as follows.
Transforming the local coordinates in (27), (28), (7) and (15) into global coordinates by (4), then substituting them into equation (6), we express the conservation of angular momentum equation in its component form as follows:

Remark 3 It can be seen from
This observation will be adopted in Section 5 for numerical implementation.

Incompressibility assumption
We assume the rod is incompressible and derive a condition for its cross section A(s, t) in this section. An incompressible material requires the total volume to be constant, i.e.: from which we get This equation can be solved by separation of variables as follows: Considering the initial condition j(s, 0) = 1 and A(s, 0) = A 0 , and noticing that A and j are both positive, the solution of (33) can be expressed as:

Finite element weak form
Equations (29) and (30) can be solved either on the reference configurations (total Lagrangian formulation) or the current configuration s (updated Lagrangian formulation). These two formulations can be transformed from one to another using equation (34), and we introduce these two formulations in this section. Let then the weak form of (29) and (30) on the current configu- with δx and δα denoting the test functions corresponding to x and α respectively. Let and Then, (36) can rewritten, in the reference configurations, as: Remark 4 It is convenient to solve a forward problem ons with updating the deformation scaler j based on a fixedpoint iteration for example, while it is convenient to solve a backward problem on s (see Section 3) because we already have the current mesh s and j can be computed directly.

The optimal control problem
In this section, we formulate an optimal control problem based on the Cosserat rod model described in the previous section. The motivation is to reconstruct the full rod configuration by computing (α, β, γ ) from observed data (x g , y g , z g ). This is an inverse problem which we formulate as a control problem. In case of low Reynolds number rods, we neglect the inertia terms and consider the following optimisation problem: reducing the discrepancy between the centreline x and an objective position given by the observed data x g = (x g , y g , z g ), by optimisation of the external force f and torque l in (29) and (30).
Problem 1 (piecewise-in-time control) Given the state variables x n−1 and α n−1 at the previous time t n−1 (n = 1, 2, . . .), and an objective position vector x g (t n ) of the worm's centreline at current time t n , subject to and In the above, Λ κ = diag λ κ 1 , λ κ 2 , λ κ 3 and Λ ω = diag λ ω 1 , λ ω 2 , λ ω 3 are two diagonal matrices. The first term in (41) is the real objective to be minimised, and we choose λ g = 10 9 ∼ 1/ b a |x g | 2 , so that the first term would not become infinitely small during the process of minimisation. All the other terms are regularisation terms with regularisation parameters λ f , λ l , λ d , Λ κ and Λ ω . Too large regularisation parameters could make it difficult to achieve the real objective, while too small ones may cause convergence issues for the numerical scheme. These regularisation terms have different mathematical and computational purposes (see Section 3.4 for biological meanings of each term): generally speaking, we want to add reasonable constraints so that the control problem is solvable and a unique solution can be obtained; the second (λ f −term) and the third term (λ l −term) are constraints of the control variables so that they would not go to infinity; the term λ d is used to disallow arbitrary rotation angles and ensure the problem is solvable; the regularisation parameter of the last two terms Λ κ −term and Λ ω −term in (41) are diagonal matrices, whose elements are parameters to control the three components of generalised curvature and angular velocity correspondingly.

Remark 5
Because we lack data to control the rotation angle α, we set a control of the variation of α, with respect to space (Λ κ −term) and time (Λ ω −term), in the last two terms in (41), which is necessary for the solution existence and for the convergence of our numerical method as formulated in Problem 2 in Section 4.
We introduce the Lagrange multipliers (or adjoint variables)x,α to eliminate the constraints of Problem 1 (dropping the subscript 'n' for simplicity).
The following boundary conditions are also included in the above functional L: and Remark 6 We find that the proposed optimal control formulation is solvable with either a Dirichlet boundary condition α(a) = α a (first component of rotation α) or the regularisation Λ ω −term in (41). We do not have the Dirichlet data for all frames unfortunately, but we notice from the last term in (41) that α 0 must be given. Therefore, we solve the first frame using Dirichlet boundary condition α(a) = 0 instead of Λ ω −term, and from the second frame we use the Λ ω −term.

Remark 7
The solvability of Problem 1 is generally a difficult question, and there is one case we are sure is unsolvable: suppose the rod undergoes a pure twist induced only by the third component of the torque l in (43), in which case there is no way to determine the frames only using the data from the centreline. In all other cases (l l 3 = 0), the deformation of the centreline is coupled with the frames and hopefully we can detect the frames through this coupling by solving Problem 1. We shall validate this idea in numerical test 5.2. Luckily, it is parsimonious to assume l l 3 = 0 when modelling C. elegans because its longitudinal body wall muscles which may not generate l l 3 torque.
After integration by parts, L can be further expressed as follows: The following Karush-Kuhn-Tucker (KKT) conditions are the first-order necessary conditions to minimise (47). with being the GO ateaux derivative with respect to variable p along the direction q [67]. If q is an arbitrary direction from p, it is usually expressed as q = δp (variation of p) [7], in which case it is convenient to abbreviate δL(p)[p; δp] as δL(p).
Let  ([a, b]) be the subspace of H 1 ([a, b]) whose functions satisfy the Dirichlet boundary condition in (45), in particular H 1 0 ([a, b]), the homogeneous Dirichlet boundary conditions. The above optimality conditions: (48) to (50), lead to the following partial differential equations (in weak forms).

Primal equation
The optimality condition (48) gives the primal equation in its weak form as follows . Find (x, α)

Adjoint equation
The optimality condition (49) gives the adjoint equation in its weak form as follows (neglecting the variation of the matrix Q A, B and vector q 3 ). Find x,α ∈ H 1 0 , such that Remark 8 We have implemented the adjoint equation with consideration of variation of q 3 , and we found that our optimal control algorithm in Section 4 struggled to converge. It is worth investigating the reason for this convergence issue, and testing the case in which the variation of all these terms is included in the future.

Optimality equation
The optimality condition (50) gives the relation between the control force and adjoint variable:

Relation to C. elegans locomotion problem
Our locomotion dataset, obtained from a 3D microscopic set up [72] contains many different trajectories of centreline The controls, the forces and the torques, can be interpreted as the reaction force from the surrounding fluids, which is initially activated by the worm's muscles. The muscles generating locomotion in C. elegans, are called body wall muscles (see Figure 4) because they are tethered to the 'wall' (or cuticle) of the animal, acting longitudinally to contract or relax the local side of the body. As C. elegans contains 95 body wall muscles that span the entire body length, we consider the action of the muscles continuously along the body. The directionality of the muscle contraction at every point along the animal is determined by α, β and γ which themselves are unknowns as part of the control problem. To represent this muscle action, we only allow d 3 close to the tangential direction of the body and also restrict the twisting movement of the worm as captured by the λ d −term.
The regularisations given by Λ κ −term and Λ ω −term are also biologically motivated. For example, it is known from the anatomy that the left and right muscle quadrants do not receive distinct neural connections along the body and tail (i.e. the posterior two thirds of the body which lie beyond the neck of the animal). This results in the majority of bending occurring in the dorsal-ventral directions and less in the left-right directions with the exception of the head and tail. We therefore adjust the magnitude of λ κ 1 and λ κ 2 to favour solutions that have more bending around d 1 than d 2 . The longitudinal muscles also restrict the twisting motion of the body which can be considered by setting a larger λ κ 3 . These additional adjustments of constraints may result in a frame that more closely matches the anatomically meaningful frame of the animal as it moves around in 3D. Λ ω −term models the internal friction, which can also be different in three local directions based on the worm's anatomy.
For the forward simulations of biological worms, the Neumann boundary condition is usually adopted at both ends of the rod. For the backward simulations, we can fully use the information from the data and adopt appropriate Dirichlet boundary conditions as considered above in (45).

Remark 9
We have the data x g of the worm's centreline for every time frame, which means the mesh (arc length) s(s, t n ) is known at the current time frame t n . Therefore, we can directly compute the deformation scaler j = ∂ss(s, t n ). This is generally not true for a forward problem, which usually requires j to be iteratively computed.

A monolithic optimal control formulation
Substituting the optimality condition (54), specifically its strong form f =x/λ f and f =α/λ l , into equation (52), we have a monolithic scheme to solve the optimisation Problem 1 as follows.

Problem 2
Given the state variables (x n−1 , α n−1 ) at the previous time t n−1 (n = 1, 2, . . .), and an objective position vector x g (t n ) at the current time t n , compute Remark 10 For the above fixed-point iteration, a relaxation parameter 0 ≤ w ≤ 1 is introduced to stabilise the algorithm: instead of directly updating x k , α k after solving (61), a weighted w x k , α k + (1 − w) x k−1 , α k−1 is adopted. We use w = 0.5 for all our simulations.

Numerical tests
We first validate the formulation (40) for simulation of a forward problem with a time discretisation scheme introduced in Appendix A, and then apply the optimal control formulation (61) to data from a forward simulation, in which case we have the ground truth rotations of the local frames and a quantitative comparison can be performed. Finally, we apply the optimal control method to data from laboratory experiments and infer the frames of rotation. All the numerical tests are implemented using open-source library FreeFem++ [42]. For code and results, see Data Availability in declarations section below.

Forward simulation of a cantilever beam
We consider a cantilever beam with a dynamic load and reproduce the result presented in [12]. The beam's length L = 1m (which is a constant for this test due to a small deformation), with density ρ = 2.73 × 10 3 kg/m 3 , Young's modulus E = 7.10×10 10 Pa and shear modulus G = 2.69×10 10 Pa. The cross section of the beam is a rectangle with width a = 0.06m and height h = 0.04m, and the numerical shear correction factor for this cross section is set to be k = 0.833. The moment of inertia is I 11 = ρab 3 /12,  (40), is expressed as: with the natural frequency of the system ω 0 = 207.0236s −1 . The beam is discretised by 100 segments and the total computational time T = 0.06 is divided into 1000 steps. The displacement and rotation at the end of the beam are plotted in

Optimal control using data from a forward simulation
In this example, we modify the previous test of the cantilever beam so that it has similar material properties to C. elegans and undergoes a large deformation. As discussed in Section 3, we also neglect the inertia terms in equation (40) to model C. elegans locomotion. Our motivation is to generate a dataset to validate the proposed optimal control method. The new beam now has an initial length L 0 = 10 −3 m and circular cross section with radius r 0 = L 0 /40m [5]. The numerical correction factor for a circular cross section is taken to be k = 4/3 [34]. We adopt values for its Young's modulus E = 1.1 × 10 5 Pa and shear modulus G = 5.0 × 10 4 Pa [5]. Before generating the dataset by applying a complicated external force and torque, we first apply a simple force F at the end of the beam, as shown in Figure 6, to validate the approach against the analytical solution based on the Timoshenko beam theory [34]: The rod is discretised by 100 segments (for which the mesh has converged), and we compute the deflection of the rod by solving the primal equation (52). It can be seen from Figure 7 that the result of the Cosserat model agrees very well with the prediction of the Timoshenko theory for up to 10% deflection of the rod. The first test of the proposed control algorithm is to use a dataset generated by a distributed force along the rod: with F max = 10 −4 . The beam undergoes a large deformation as shown in Figure 8 (left); meanwhile a curvature κ l 1 (see formula (28)) is also generated along the rod. Using the proposed control formulation in Section 4 and control parameters of λ f = 1, λ l = 10 −10 , λ d = 10 −6 , λ κ 3 = 10 −10 (λ κ 1 = λ κ 2 = 0) and Λ ω = 0, both the position and the curvature can be recovered accurately as shown in Figure  8. Note that all the other components of the position vector and generalised curvature are zero although they are not presented here. Before moving to other test cases, let us test the convergence of the objective x − x g , the output curvature κ l − κ l f (κ l f is from the forward simulation), the tangential direction ∂ s x − d 3 as well as the algorithm itself measured by the relative error of (x, α) between the current and previous fixed-point iterations, with regards to the control parameters λ f , λ l and λ d .
The main findings are summarised as follows: (1) If we only use f (notice that this does not mean λ l = 0; we have to remove l term in (61)) as the control, the proposed algorithm struggles to converge no matter how we play with the other parameters. Therefore, the regularization λ l − term in (41) does play an important stabilisation role, although l = 0 when we generate the dataset; (2) If we plot the above convergence measures as shown in Figure 9 (λ f = 1, λ d = 10 −6 and λ κ 3 = 10 −20 ) and vary λ l from magnitude 10 −20 to 10 3 , we find that these convergence curves are exactly the same. The only difference we notice is that the magnitude of the adjoint variable δα varies correspondingly from 10 −13 to 10 10 so that the control torque l = δα/λ l always has a magnitude of 10 −7 . Therefore, the algorithm (at least for this test) is not sensitive to the regularisation parameter λ l , although it is required to stabilise the algorithm as pointed out above; (3) The proposed algorithm can converge stably with the regularisation parameter λ d varying from 10 1 to 10 −10 , and the convergence of κ l − κ l f and ∂ s x − d 3 is plotted in Figure 10 with λ f = 1, λ l = 10 −10 and λ κ 3 = 10 −20 ; (4) The proposed algorithm converges for a range of parameters λ f from 10 −10 to 10 6 . Despite the steady convergence of the algorithm, the value of λ f must be sufficiently small for the objectives to be sufficiently reduced. To demonstrate this, convergence plots for the relevant quantities given two extreme values (10 −10 and 10 6 ) of λ f , are compared in Figure 11; (5) The purpose of the λ κ 3 parameter is to control twist along the rod. In this test, the algorithm can converge stably and the curvature error can be reduced sufficiently with λ κ 3 from 1 to 10 −15 .
We next consider a dataset generated by the following force and torque, which creates both curvature and torsion along the rod as shown in Figure 12.
Again we test all the regularisation parameters systematically, and our findings are summarised as follows: (1) λ l is necessary for stability, and it can be chosen from a magnitude of 10 −30 to 10 −2 (based on a test with λ f = 1, λ d = 10 −6 and λ κ 3 = 10 −20 ); (2) the recommended value for λ d ranges between 10 −10 and 10 3 , otherwise the objective cannot be sufficiently reduced when λ d > 10 3 , or the algorithm cannot converge when λ d < 10 −10 (based on a test with λ f = 1, λ l = 10 −6 and λ κ 3 = 10 −20 ); (3) the proposed algorithm can converge steadily for a range of λ f from 10 −30 to 10 2 , and the suggested values are λ f < 1, otherwise the objective cannot be reduced sufficiently (based on test of λ l = 10 −6 , λ d = 10 −6 and λ κ 3 = 10 −20 ); (4) The algorithm converges stably with λ κ 3 for magnitude of 10 −15 to 1, but the recommended value is > 10 −25 otherwise a too large κ l 3 could be generated. A convergence of relevant quantities with a specific parameter set is plotted in Figure 13. The comparison of the position and curvature between the forward and backward computations are displayed in Figure 12 and 14 respectively.  It can be seen that the positions between the forward and backward simulations match very well along the rod, and the curvatures also match well except the end where the Dirichlet boundary condition is applied.

Remark 11
We find that non-trivial forward data, involving large bend and torsion for example, is not easy to generate, because it is not straightforward to provide or design a force f or torque l so that the forward problem can converge easily. While once a dataset is given, the control problem is easy to converge -converging to the same position vector x and rotation α (as the forward simulation results) even with a different control force and torque. This can be understood and by noting that we do not expect that the force and torque (producing the same x and α) are unique. As an example, we show the control force and torque for the previous test in Figure 15, from which it can be seen that the magnitude is similar to the designed one in (65), but the distribution is different.

Reconstruction of C. elegans locomotion based on experimental data
We tested the proposed method on three examples of C. elegans locomotion in 3D volumes [72]. The data represent the body-midlines that were reconstructed from microscopyvideo footage of freely and spontaneously moving worms that were immersed in different fluids. We present one test case in this section and all the three tests (including the dataset, FreeFem++ code and simulation results) can be found from public GitHub repository: https://github.com/ yongxingwang. The first test case consists of a sequence of 1160 time frames with a sampling interval of 0.04s (as inertia is not considered in our model and the time step only appears in the regularisation term in (41)), and each reconstructed body centerline consists of 128 discrete three-dimensional spatial points. These examples are of interest due to the threedimensional postures and motion of the swimmer: in this clip, the worm exhibits large bend, large torsion and moves forward and backward in a three-dimensional space for about 46s. The physical body of the worm is modelled by a cylindrical Cosserat rod with a circular cross section of initial radius 2 × 10 −5 m, Young's modulus E = 1.1 × 10 5 Pa and shear modulus E/(1 + ν)/2Pa with ν = 0.4 [5,20]. Four typical postures of the worm are shown in Figure 16: the worm initially moves from the right to the left and starts a manoeuvre to reverse its motion at around the 400 th time frame; after another 350 steps, the worm suddenly bends to resemble a capital Ω (left-bottom in Figure 16) and moves to the right.
To construct the first frame, we set α a = 0 and Λ ω = 0; and from the second frame, we use a non-zero Λ ω (at least non-zero λ ω 3 for the sake of convergence) without any Dirichlet data. Three components of generalised curvature are plotted along the worm's body for all the time frames in a two-dimensional plane as shown in Figure 17, from which a bending (κ l 1 and κ l 2 ) wave can be seen propagating from the worm's head to the tail; the twisting (κ l 3 ) wave is not obvious but some twisting can still be observed. The propagated wave is consistent with the moving direction of the worm as shown in Figure 16 and analysed in the above. Starting from a converged parameter set as shown in Table  1 (for the results in Figures 16 and 17), with the converged objectives in (41) being shown in Figure 18, we vary these parameters, study the convergence of the algorithm and compare corresponding results in the following. Notice that this set of parameters has the minimal non-zero parameters to make sure the proposed algorithm can converge (please refer to Remark 5 and 6 and Section 3.4 for explanations).
Parameter λ f : with other parameters frozen, the proposed algorithm converges stably for λ f from 10 −20 to 10; the fixed-point iteration becomes slower for larger λ f . We find that all the objectives stay the same except the control f (see Figure 19).
Parameter λ l : we then keep λ f = 10 −20 and the proposed algorithm still converges stably with the magnitude of λ l from 10 −10 to 1. We observe that all the objectives are still almost the same except the control l as shown in Figure 19.
Parameter λ d : based on the above two tests, we realise that the total objective function (41) is dominated by the λ d −term.
With other parameters frozen, we find that the convergence range for λ d is approximately between 10 −2 and 10 3 . A comparison of converged objectives between Parameter-0 in Table 1 and its variation case with λ d = 10 −2 (reduced from λ d = 10 2 ) is plotted in Figure 20, from which it can be seen that (i) the real objective x − x g / x g is reduced by two orders of magnitude, with oscillations for some frames which is expected for such a small regularisation parameter; (ii) the  Table 1 Parameter-0: minimal non-zero parameters for the sake of convergence of the proposed algorithm 10 −3 10 −6 10 2 diag (0, 0, 0) diag 0, 0, 10 −20 λ d −term increases as its regularisation parameter decreases from λ d = 10 2 to λ d = 10 −2 . We notice that the magnitude of ∂ s x − d 3 / d 3 increases to 10 −1 for some frames, which means d 3 detaches from the tangential ∂ s x of the centreline 10%. For example, Figure 21 shows frame number 700 where the normal direction d 3 of the cross section detaches from the tangential direction ∂ s x of the centreline; (iii) none of the other objective terms show a significant change except the control f which varies according to x − x g / x g as expected.

Remark 12
The detachment of d 3 from ∂ s x is an important feature of Cosserat rods. Otherwise, the Cosserat rod approaches to the Kirchhoff rod (if the deformation scaler j = 1, which is true for our case study system of C. elegans: we observe that | j − 1| < 10 −3 always holds numerically), in which case it is assumed that d 3 = ∂ s x [31].

Remark 13
In keeping λ d = 10 −2 (which now cannot dominate the total objective function of (41)) and varying λ f and λ l , we again can observe a variation of other objective terms although we would not present all these tests here. However, too small regularisation parameters can cause stability issues as we have already seen when reducing λ d from 10 2 to 10 −2 although the algorithm still converged.
The parameter Λ κ provides a constraint of the curvature along the worm's body, which allows us to consider the anatomical muscle structure of C. elegans as explained in Section 3.4. Based on Parameter-0, we now simply choose Λ κ = diag 0, 0, 10 −20 to restrict the twist motion of the worm, and the generalised curvature is plotted in Figure 23. Comparing with Figure 17, we can see that not only does the magnitude of κ l 3 dramatically decrease, but a clear twisting wave also appears along the worm. In addition, there is also an influence on the first two components κ l 1 and κ l 2 : the wave propagation is clearer although the curvature magnitudes are almost the same as the case of Λ κ = 0. We also notice that the reversing manoeuvre (starting at around frame 700) becomes more distinct: κ l 3 is much larger near the worm's head at frame 700 than elsewhere or any other frames, as can be observed from Figure 22 and 23. Similarly, using non-zero λ κ 1 or λ κ 2 would allow us to favour bending in the dorsal-ventral directions which is consistent with the C. elegans neuromusculature.
Λ ω is used to model the internal friction of the worm which is necessary from the second time frame. For the sake of convergence of our algorithm, only λ ω 3 is required. If we increase Λ ω in Parameter-0 from diag 0, 0, 10 −20 to some value less than diag (1, 1, 1), the proposed algorithm converges stably and the angular velocity stays almost the same as shown in Figure 24. However, there is a big change in the curvature (see Figure 24 and Figure 25): wave propagation is no longer apparent because larger Λ ω tends to keep the rotation angles (along the body and overtime) the same in time, consequently the space derivative (curvature) along the worm's body does not change significantly in time.

Conclusion and discussion
This paper presents three contributions: the forward formulation of a Cosserat rod, the optimal control method and the reconstruction of C. elegans locomotion.
The forward formulation of Cosserat rod is developed from [12], in which the Cosserat rod is described by three components of the position vector (x, y, z) and three components of the rotation vector (α, β, γ ). We derive the angular velocity and generalised curvature using a new method in Section 2.3 and rewrite all the control equations in a matrix-    [12] and found that this formulation is robust and convenient for analysis of complex dynamic behaviour of slender rods.
A well-posed inverse problem may be solving for external forces and torques, given both the position vector and rotation vector. However, accessing both position and rotational information may not be practical. For example, for a biological worm moving freely in a fluid environment, it is difficult to measure the worm's local orientation (rotation vector) while its centreline (position vector in the global frame) can be reconstructed from video footage [72]. A complementary problem that may be tackled with an analogous approach considers a robotic worm exploring an unknown space. Given sufficient local body sensors, the body posture (including bending and twisting) would be reliably detected by such a robot. However, in the absence of external location data, the robot may lack positional information. Motivated by these biological and engineering problems, we consider an ill-posed inverse problem which solves for rotation vector, external forces and torques given only the position vector. We present a robust and efficient optimal control method to solve this inverse problem: the objective is to minimise the discrepancy between the position vector and a given centreline of the Cosserat rod, and the control variables are the external forces and torques, with regularisations of the rotation vector. The regularisation terms provide constraints of the rotation vector so that the inverse problem is solvable. We have tested the proposed optimal control formulation using data from forward simulations and shown that the rotation vector can be accurately computed with appropriate and controllable regularisation parameters.
The proposed optimal control is applied to reconstruction of C. elegans locomotion based upon its centreline data from laboratory recordings. The solvability of this challenging inverse problem relies on meaningful regularisation terms. The proposed approach allows us to add different terms conveniently to model C. elegans' neuromusculature, and our inverse model is demonstrated to be robust to a range of regularisation parameters. There are five parameters (nine if considering components of Λ κ and Λ ω ) as shown in Table 1, which also indicates the minimally-required non-zero parameters for the sake of convergence of our proposed method. λ f and λ l correspond to the control force f and torque l respectively, which can stop f and l becoming infinite and have the effect of stabilising the proposed method. λ f and λ l can be robustly chosen from a range of values based on our numerical experiments, without a significant influence on the main outputs such as the centreline x and curvature κ. λ d has a biological and anatomical grounding because it keeps the normal direction d 3 of the worm's cross section close to the tangential direction ∂ s x of its body's centreline. It also has a numerical effect of stabilising the proposed method and differentiating the Cosserat rod and Kirchhoff rod models: larger λ d tends to force d 3 to be the same as ∂ s x (approaching the Kirchhoff rod consequently).
λ ω 3 in Λ ω = diag λ ω 1 , λ ω 2 , λ ω 3 has a clear numerical purpose, because our method needs it to be non-zero for convergence. The other two components of Λ ω and all the three components of Λ κ = diag λ κ 1 , λ κ 2 , λ κ 3 can be zero. However, setting Λ ω and Λ κ to non-zero values allows us to model the muscular of C. elegans as pointed out in Section 3.4.
Several interesting topics have been stimulated by this study, which are briefly summarised as follows: The proposed optimal control formulation is based on a combination of laboratory data and modelling of the worm's muscle structure. If we can collect the data of at least one cross section's movement of C. elegans, we can then apply a Dirichlet boundary condition of α and use less regularisation terms as commented in Remark 6 (λ ω 3 can be zero then). Measuring the movement of cross sections of a hair-thin C. elegans in laboratory is technically difficult. However, setting a mark and following one cross section may be possible and would provide a possibility to validate the predictions of our inverse model based on the above assumptions.
Having computed the local frames (rotation vectors), we can then formulate these rotations into the objective function, and compute the external force f and torque l without regularisation λ d −, λ ω −, and λ κ − terms. These force and torque terms will provide us C. elegans' muscle force quantitatively, which will help us to model and understand its neuromuscular system.
One more interesting topic is modelling time evolution. In this paper, a friction term (λ ω − term) is introduced to link different time frames. An alternative approach is to introduce a viscoelastic constitutive model as adopted in [53], which is expected to more appropriate for modelling nematodes locomotion [5].
Another way to formulate the underlining problem is to apply a model for the force f, such as slender body theory [49], then only use l as a control variable. Hopefully, this would lead to a well-posed problem without additional regularisation terms.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Time discretisation
We introduce a time discretisation scheme for weak form (36) or (40). Let M and S be the mass and stiff matrix respectively after spacial discretisation, and F be force vector, the spatial discretisation of (36) or (40) leads to the following algebra partial differential system: If the time domain is disretised as t 0 = 0, t 1 , . . . with a time frame of Δt = t n − t n−1 , then (66) may be discretised as follows.