# A continuum active structure model for the interaction of cilia with a viscous fluid

Astrid Decoene, Sébastien Martin, Fabien Vergnet

► To cite this version:

Astrid Decoene, Sébastien Martin, Fabien Vergnet. A continuum active structure model for the interaction of cilia with a viscous fluid. Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 2023, 10.1002/zamm.202100534. hal-02493513

HAL Id: hal-02493513

https://hal.science/hal-02493513

Submitted on 27 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A continuum active structure model for the interaction of cilia with a viscous fluid

Astrid Decoene\*, Sébastien Martin\*, Fabien Vergnet\*

February 27, 2020

## Abstract

This paper presents a modeling, analysis and simulation of a fluid-structure interaction model with an active thin structure, reproducing the behaviour of cilia or flagella immersed in a viscous flow. In the context of linear or nonlinear elasticity, the model is based upon the definition of a suitable internal Piola-Kirchoff tensor mimicking the action of the internal dyneins that induce the motility of the structure. In the subsequent fluid-structure interaction problem, two difficulties arise: on the one hand the internal activity of the structure which leads to more restricted well-posedness conditions and, on the other hand, the coupling conditions between the fluid and the structure that require a specific numerical treatment. In the context of numerical simulation, a weak formulation of the time-discretized problem is derived in functional frameworks that include the coupling conditions but, for numerical purposes, an equivalent formulation using Lagrange multipliers is introduced in order to get rid of the constraints in the functional spaces: this new formulation allows for the use of standard (fluid and structure) solvers, up to an iterative procedure. Numerical simulations are presented, including the beating of one or two cilia in 2d, discussing the competition between the magnitude of the internal activity and the viscosity of the surrounding fluid.

**Keywords.** Stokes flow, elasticity, fluid-structure interaction, active structure, numerical simulation
MSC. 74F10, 76Z05, 76D07, 74B20, 65M60, 74H15

## 1 Introduction

Cilia and flagella are motile elongated structures, involved in swimming and/or transport mechanisms that arise in many living organisms. Flagella are usually used by micro-swimmers such as sperm-cells, bacteria or algae for motility purpose at low Reynolds number, while cilia are generally involved in the transport of proteins, nutrients or dust inside bigger organisms. At the origin of all these mechanisms are two essential ingredients. The first one is the capacity for cilia and flagella to modify their shapes by generating active internal deformations and stresses, even without external load. The second one is the strong reciprocal interaction between these structures and the surrounding fluid. The problem we are interested in is the capacity of such microorganisms to deform themselves by mean of internal biological motors and to interact with the surrounding fluid. In the present article, we present a model for the actuation of elongated cilia-like structures, which fits within the framework of continuum mechanics, and study the fluid-structure interaction problem they are involved in.

Eukaryotic cilia (or flagella) are elongated deformable structures with a typical diameter between 0.1 and 0.3 µm, whereas the length of cilia can vary from 5 µm (in the lung [1]), to 80 µm long (for the tail of spermatozoon in mice). Cilia are membrane-bounded structures composed of a microtubule cytoskeleton, called axoneme, consisting of a ring of nine doublets microtubules surrounding a central pair of microtubules. Outer doublets microtubules are linked to the central pair by radial proteins and to each other by nexin links, which strengthen the structure. Beating movements of cilia are induced by internal motors producing a bending in the whole structure, when two outer doublets microtubules slide with respect to one another. This sliding is produced by proteins, called dyneins, that are synthesized on one doublet microtubules and attach to the neighboring one. The radial connections resist the

<footnote>\* Université Paris-Sud, LMO (CNRS UMR 8628), Bâtiment 307, 91405 Orsay cedex, France.
astrid.decoene@math.u-psud.fr</footnote>

<footnote>† Université de Paris, MAP5 (CNRS UMR 8145) & FP2M (CNRS FR 2036), 45 rue des Saints-Pères, 75270 Paris cedex, France.
sebastien.martin@parisdescartes.fr</footnote>

<footnote>‡ École polytechnique, CMAP (CNRS UMR 7641), Route de Saclay, 91128 Palaiseau cedex, France.
fabien.vergnet@polytechnique.edu</footnote>
sliding and contribute to the bending of the structure, since the cilium is anchored at the bottom. This mechanism appears all along the length of a cilium and between all doublets microtubules, which contributes to the emergence of different sliding patterns and different shapes of deformation for the cilium. The precise nature of the spatial and temporal control mechanisms regulating the various ciliary beats is still unknown [2]. Further details on the internal structure and mechanisms of cilia can be found in [3].

The presence of ciliary propulsion in almost all living organisms, from bacteria to mammals, has encouraged numerous scientists to model and study this universal phenomenon. The first work in that sense goes back to 1951 and is due to Taylor [4], who initiated the mathematical study of microorganism propulsion. This work presents the swimming of a extensible sheet in a viscous fluid modeled by the Stokes equations. The author works in the rest frame of the sheet, whose deformations are modeled by the propagation of a wave of small amplitude. Thus, the unknown of the problem is the velocity of the fluid far from the sheet, which also represents the velocity of the sheet in the laboratory frame. Extensions to an infinite cylinder [5, 6] and to finite objects [7, 8] have subsequently been studied. When the amplitude of deformations is large, the resistive force theory (also known as local drag theory) developed in the pioneering work of Gray and Hancock [9], describes the cilia as several cylinders and uses the linear property of the Stokes equations to compute the flow induced in the fluid. A similar but more accurate method is the slender body theory, started by Hancock in [10] and then improved by Lightill in [11], which makes use of the long and thin geometry of cilia. It consists in modeling a cilium by a distribution of stokeslets and dipoles, which impose a force on the surrounding fluid. This method is also known as the sublayer method or the stokeslet method and has been extensively applied to the simulation of thin flagellar propulsion when the deformations of the structure are imposed (see for example [12, 13, 14]). More recently, the immersed boundary method have been used for the simulation of thin beating cilia in a Newtonian fluid in [15]. As in the stokeslet method, the idea is to impose a distribution of forces in the fluid. However, in that case the force does not come from the slender body theory but is used to impose the equality of the fluid and solid velocities on the fluid-structure interface. In this work, the velocity of the structure is imposed and its action on the fluid is studied. A different approach is also considered in [16], where cilia are three-dimensional structures whose deformations are reproduced by solving a one-dimensional transport equation. The fluid velocity on the fluid-structure interface is imposed using a penalization method, thus no retro-action from the fluid to cilia is taking into account. In [17] cilia are modeled by thin structures whose deformations are given by the same transport equation than in [16]. The equality of the fluid and solid velocities on the fluid-structure interface is imposed with the immersed boundary method and the action of the fluid on the structure is taken into account by changing the velocity of the structure, but cilia always follow the same beat pattern. In all the works that have been previously mentioned, the cyclic shape change of a cilium is imposed whereas its beatform should be an emergent property of a coupled system involving the internal mechanisms of the cilium, the elastic properties of the structure and the surrounding viscous fluid.

The first work that attempted to take into account the internal activity of cilia through local deformations is due to Machin [18]. The cilium is considered as an elastic filament immersed in a viscous fluid, whose action on the structure is given by the resistive force theory. Moreover the internal activity of the cilium is modeled by adding an active bending moment distributed all along the structure. With the active bending moment, wave-like displacements are observed whereas, with a passive elastic filament driven from its proximal end, the forms of the wave do not match the actual shapes of cilia. Thus Machin brought to light the importance of local contractility in the deformation of cilia; in [19] a modification of this model is proposed, by considering two active filaments with regular cross-connections whose contractility in activated when the passive bending reach a critical value. In [20] a similar three-dimensional model is proposed with a more realistic geometry of the internal structure of a cilium. Subsequently, several class of models have been proposed such that curvature-controlled models [21, 22, 23] and self-oscillatory models [24, 25]. A comprehensive review of these different models is presented in [26].

Another approach is presented in [27] and further studied in [28, 29, 30], where a discrete description of the axoneme is proposed in two space dimension. The cilium is composed of elastic filaments connected by a finite number of springs that represent nexin and dynein links. Then, the deformation of the structure is produced by the connection scenario of dynein links which depends on the geometry of the structure. At the fluid-structure interface the continuity of the velocity is considered and is treated numerically with the immersed boundary method. A similar model is considered in [31] and [32] in three dimension, where a precise description of the “9+2” structure of the cilium is proposed. In both works, the emerging beating patterns are realistic.

Unlike all previous works on cilia and self-propelled microorganisms moving in a viscous fluid, we aim to model the behavior of active biological structures in the framework of continuum mechanics, without using a detailed description of their internal structure. The reasons for this study are twofold. First, since the chemical, biological and even mechanical processes for the internal activity of eukaryotic cilia are not yet completely understood, we do not intend to model the nexin and dynein links at the nanometric scale. Instead, we rather take into account
the activity in a more phenomenological manner by mean of an internal stress. The context of two and three dimensional elasticity is particularly suitable to reproduce realistic deformations of cilia and flagella. Second, the framework of continuum mechanics enables to fully consider the fluid-structure interaction, which is one of the most important ingredient of the system and which is often neglected in other studies. The model that we develop in this paper is suitable for both the mathematical study and the numerical simulation of the fluid-structure interaction with active structures and a viscous fluid.

The paper is organized as follows. In Section 2 we present the equations of active elasticity that have been introduced in biomechanics for the study of biological tissues, but never used or mathematically studied for active microorganisms at low Reynolds number. Examples of activity scenarios are illustrated as well. In Section 3, we couple the elasticity equations to the Stokes equations and study the fluid-structure interaction problem. For the numerical simulation of active structures beating in a viscous fluid we introduce a saddle-point formulation of the problem, where the condition of equality of the fluid and structure velocities on the fluid-structure interface is treated by a Lagrangian multiplier. In particular, this enables the use of standard finite element methods and solvers. In Section 4, we present the numerical resolution process and some numerical results for one and two cilia with prescribed internal activity. The influences of the viscosity of the surrounding fluid, the phase shift between internal activities and the distance between cilia are investigated.

# 2 A continuum active structure model

The purpose of the present section is to develop a macroscopic model for the internal activity of cilia-like structures. We start by giving a brief review of these models.

## 2.1 The active-stress method

For the study of biological tissues, two popular macroscopic approaches are used to model the muscles activity, namely the active-stress and active-strain methods (see [33] for a review). The former consists in adding an active component to the passive stress tensor usually derived from the strain energy law, while the later adopts a multiplicative decomposition of the tensor gradient of deformation in which the activation acts as a pre-strain. Both techniques have been extensively used in myocardium, arteries and even face muscles studies ([34, 35, 36, 37]), but, to our knowledge, not in the context of microswimmers.

When the structural organization of the active components in the tissue is known, but rather complicated to model individually (as it is the case for cilia), the active-stress method appears to be a more suitable approach. Indeed, in this case, the internal stress can be approximated by averaging the geometric arrangement of the active elements at the micro scale, in order to exhibit a macroscopic fiber-like structure. Then, we suppose that the active behavior of the tissue is only due to elastic deformations in the direction of these fiber-like structures, which are called active fibers. More precisely, if $e_a$ denotes a unit vector field in the direction of active fibers within the tissue, which depends on the material position and the time, the active stress tensor, which is denoted by $\Sigma^*$, writes

$$ \Sigma^* = \Sigma_a e_a \otimes e_a, $$

where $\Sigma_a$ is a scalar function, which also depends on the time and the material position, and $\otimes$ denotes the tensor product. Thus, in the general case, the internal activity is given by a scalar function $\Sigma_a$, that we call the activity scenario, and by a unit vector field $e_a$ which, at each material point, points in the direction of active fibers. For more information on models for contractile organs, see [38, Chapter 2] and references therein.

## 2.2 The problem of active elasticity

Let $\Omega_s$ be a Lipschitz open connected bounded subset of $\mathbb{R}^n$, with $n \in \{2, 3\}$. Its boundary, $\partial\Omega_s$, is divided in two parts denoted $\Gamma$ and $\Gamma_s$ such that the boundaries satisfy $\partial\Omega_s = \Gamma \cup \Gamma_s$ and $\Gamma \cap \Gamma_s = \emptyset$. Moreover we denote by $n_s$ the exterior unit normal vector to $\Omega_s$. We suppose that $\Omega_s$ is filled with an elastic active medium, subjected to a time dependent body force, denoted by $f_s$. The internal activity of the structure $\Sigma^*$ is described using the active-stress method and is supposed to depend only on the time and the material position. Then, the quasi-static problem of active elasticity with homogeneous Dirichlet and Neumann boundary conditions, is to find the displacement of the structure $d_s$, solution for all time $t \ge 0$, of the following set of equations:

$$ \begin{aligned} -\mathrm{div}((I + \nabla d_s(t))(\Sigma_s(d_s(t)) - \Sigma^*(t))) &= f_s(t), && \text{in } \Omega_s, \\ (I + \nabla d_s(t))(\Sigma_s(d_s(t)) - \Sigma^*(t))n_s &= 0, && \text{on } \Gamma, \\ d_s(t) &= 0, && \text{on } \Gamma_s. \end{aligned} \quad (1) $$
The matrix I is the identity matrix of $\mathbb{R}^n$ and $\Sigma_s(d_s(t))$ is the so-called second Piola-Kirchhoff stress tensor at time $t$, which describes the passive elastic behavior of the structure. For simplicity, we will always assume that the elastic medium follows the Saint Venant-Kirchhoff law, i.e. that the second Piola-Kirchhoff stress tensor writes

$$ 
\begin{aligned}
\Sigma_s(d_s(t)) &= 2\mu_s E(d_s(t)) + \lambda_s \operatorname{tr}(E(d_s(t)))I, \\
E(d_s(t)) &= \frac{1}{2}(\nabla d_s(t) + \nabla d_s(t)^T + \nabla d_s(t)^T \nabla d_s(t)),
\end{aligned}
\qquad (2) $$

where $\mu_s > 0$ and $\lambda_s > 0$ are Lamé's parameters and $E(d_s(t))$ is known as the Green-Lagrange strain tensor at time $t$. The elasticity parameters are usually given by mean of Young's modulus $E_s$, which represents the stiffness of the medium, and Poisson's ratio $\nu_s$, which represents its compressibility:

$$ \mu_s = \frac{E_s}{2(1 + \nu_s)}, \quad \lambda_s = \frac{E_s \nu_s}{(1 + \nu_s)(1 - 2\nu_s)}. $$

Problem (1) differs from the classical elasticity equations by the presence of the stress tensor $\Sigma^*$, which acts in two different ways on the structure. First, it modifies the resulting forces that act on the structure by adding a body force which writes $\operatorname{div}(\Sigma^*(t))$ in $\Omega_s$ and a surface force which writes $\Sigma^*(t)n_s$ on $\Gamma$. Second, it modifies the elasticity operator by adding the term $\operatorname{div}(\nabla d_s(t)\Sigma^*(t))$ in the left-hand side of the continuity equation.

In particular, if the body force $f_s$ is null, a displacement $d^*$ which satisfies

$$ 
\begin{aligned}
\Sigma_s(d^*(t, x)) &= \Sigma^*(t, x), & \forall t \ge 0, & \forall x \in \Omega_s, \\
d^*(t, x) &= 0, & \forall t \ge 0, & \forall x \in \Gamma_s,
\end{aligned}
 $$

is a solution of problem (1). This means that the internal activity acts as a constraint on the second Piola-Kirchhoff stress. At the infinitesimal scale, the second Piola-Kirchhoff stress describes the forces in the reference configuration that each particle of the elastic medium applies on its neighbors by unit area in the reference configuration. Thus, the active stress tensor $\Sigma^*$ can be seen as a constraint on the internal forces that neighboring particles exert on each other.

## 2.3 Application to cilia-like structures

In the present study, a cilium-like structure is supposed to be an elastic active medium whose passive component satisfies the Saint Venant-Kirchhoff law and whose reference configuration, $\Omega_s$, is a straight vertical cylinder of finite length. Moreover, the structure is supposed to be anchored at its bottom boundary, denoted $\Gamma_s$, and we denote by $\Gamma$ the remaining of the boundary. In order to model the internal activity of the cilium-like structure, we apply the active-stress method, based on the knowledge of the biological structure of cilia. As we explained, the bending mechanics inside a cilium comes from the activation of several molecules which ends up in the sliding of the microtubules, the elongated rod-like structures located at the periphery of the cilium. At a more macroscopic scale, these dynamics can be seen as local elastic deformations in the direction of the microtubules which induce, because the cilium is anchored at the bottom, a bending deformation.

As a consequence, we suppose that a cilium is embedded with vertical active fibers. Then the unit vector field $e_a$ is constant, in time and in space, and the active stress tensor $\Sigma^*$ is given by

$$ \Sigma^*(t, x) = \Sigma_a(t, x)e_a \otimes e_a, \quad t \ge 0, \quad x \in \Omega_s, \qquad (3) $$

where the activity scenario $\Sigma_a$ is a scalar field which only depends on the time and the material position in the reference configuration $\Omega_s$. In particular, if the activity scenario is constant in time and in space, this induces an elongation or a contraction of the whole structure in the direction given by $e_a$, depending on the sign of $\Sigma_a$. If $\Sigma_a$ is positive the solid stretches, whereas if $\Sigma_a$ is negative it shrinks. More generally, at a given time $t$ and at a given point $x$, the sign of $\Sigma_a(t, x)$ indicates whether the structure locally expanses or contracts. Thus, in the case of cilia-like bodies, we propose a model for active structures that only depends on the choice of an activity scenario $\Sigma_a$. In particular, this enables to easily reproduce biomimetic self-induced deformations of elongated elastic structures.

## 2.4 Examples of internal activity

In this subsection, we aim to mimic the characteristic flapping deformations of cilia and flagella. Since the structure is anchored at the bottom, the local expansion or shrinking of the medium should induce the bending of the whole structure, if the activity scenario is well-chosen. Let $L_c > 0$ be the length of the cilium-like structure, $r_c > 0$ be its
radius and $x_c$ be its mean position on the abscissa axis. Then, in two space dimensions, $\Omega_s$ is the rectangle defined by

$$ \Omega_s = \{(x_1, x_2) \in \mathbb{R}^2; x_c - r_c \le x_1 \le x_c + r_c, 0 \le x_2 \le L_c\}. \quad (4) $$

Moreover, we suppose that the structure is anchored at $x_2 = 0$. We recall that in the case of a cilium-like structure, all active fibers are oriented in the direction of the vector $e_a = (0, 1)$.

**Bending** The first scenario that we study is the case of the periodic (in time) bending of a two-dimensional structure. We consider a scenario which only depends on the time and on the first coordinate $x_1$ and which is proportional to the difference $x_c - x_1$. It writes

$$ \Sigma_a(t, (x_1, x_2)) = \frac{C_a}{L_c r_c} \sin(2\pi f_a t)(x_c - x_1), \quad \forall t \ge 0, \forall (x_1, x_2) \in \Omega_s, \quad (5) $$

where $f_a$ is the beating frequency and $C_a > 0$ is the intensity of the internal activity. Actually, if the sign of $\sin(2\pi f_a t)$ is positive, the structure locally stretches in the half-domain defined by

$$ \{(x_1, x_2); x_c - r_c \le x_1 < x_c, 0 \le x_2 \le L_c\} $$

and locally shrinks in the half-domain defined by

$$ \{(2x_c - x_1, x_2); x_c < x_1 \le x_c + r_c, 0 \le x_2 \le L_c\}. $$

Thus, because the cilium is anchored at the bottom, it ends in the bending of the whole structure to the right. On the contrary, if the sign of $\sin(2\pi f_a t)$ is negative, the cilium bends to the left. To illustrate this bending behavior, we numerically solve problem (1) with $\Sigma_a$ defined by (5) and $e_a = (0, 1)$. To that aim we use the finite element method with $P_1$ elements and a Newton solver to handle the nonlinear elasticity operator. The resulting deformations of the structure are presented in Fig. 1, with the set of parameters be given in Table 1.

<table><caption>Table 1: Set of parameters for the bending scenario of activity.</caption><thead><tr><th>L<sub>c</sub> (μm)</th><th>r<sub>c</sub> (μm)</th><th>x<sub>c</sub> (μm)</th><th>C<sub>a</sub> (pN·μm<sup>2</sup>)</th><th>f<sub>a</sub> (Hz)</th><th>E<sub>s</sub></th><th>ν<sub>s</sub> (pN·μm<sup>-2</sup>)</th></tr></thead><tbody><tr><td>6.5</td><td>0.2</td><td>0</td><td>1.3</td><td>10</td><td>10<sup>6</sup></td><td>0.49</td></tr></tbody></table>

The elasticity parameters $\mu_s$ and $\lambda_s$ are given by Young's modulus $E_s$ and Poisson's ratio $\nu_s$. The chosen value of Poisson's ratio means that the structure is nearly incompressible. For Young's modulus, this value comes from experimental studies [39]. The activity scenario is plotted on the mesh of the structure. For a time $t$ between 0s and 0.025s, we observe in Fig. 1 that the cilium bends to the right, since the sign of $\sin(2\pi f_a t)$ is positive. For greater values of $t$, the structure returns to its reference configuration while $\sin(2\pi f_a t)$ decreases and starts to bend to the left when the sign of $\sin(2\pi f_a t)$ becomes negative. At $t = 0.1$s the structure is back in its reference configuration and is about to bend to the right one more time.

**Flapping** Now we study the case of a two-dimensional flapping structure, by making the activity scenario depend also on the second coordinate $x_2$. We thus consider an activity scenario representing a wave propagation along the length of the cilium, defined by

$$ \Sigma_a(t, (x_1, x_2)) = \frac{C_a}{L_c r_c} \sin(2\pi f_a (x_2 - v_a t))(x_c - x_1), \quad \forall t \ge 0, \forall (x_1, x_2) \in \Omega_s, \quad (6) $$

where $f_a$ denotes the frequency of the wave, $C_a$ still denotes the intensity of the internal activity and $v_a$ is the velocity of propagation of the wave along the length of the cilium. In this situation, elongations and contractions of the medium are not constant along the length of the cilium and it results in undulations of the whole structure. As before, we numerically solve problem (1) with $\Sigma_a$ defined by (6) and $e_a = (0, 1)$, using the finite element method and Newton's algorithm. The results are shown in Fig. 2 with the set of parameters given in Table 2. From Fig. 2 a)

<table><caption>Table 2: Set of parameters for the flapping scenario of activity.</caption><thead><tr><th>L<sub>c</sub> (μm)</th><th>r<sub>c</sub> (μm)</th><th>x<sub>c</sub> (μm)</th><th>C<sub>a</sub> (pN·μm<sup>2</sup>)</th><th>f<sub>a</sub> (μm<sup>-1</sup>)</th><th>E<sub>s</sub> (pN·μm<sup>-2</sup>)</th><th>ν<sub>s</sub></th><th>v<sub>a</sub> (μm·s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>6.5</td><td>0.2</td><td>0</td><td>1.2</td><td>1</td><td>10<sup>6</sup></td><td>0.49</td><td>1</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1: Bending of an elongated elastic structure with the internal activity given by (5) at different times: a) $t = 0$s, b) $t = 0.012$s and c) $t = 0.025$s. The activity scenario $\Sigma_a$ is plotted, with the set of parameters be given in Table 1.</figcaption></figure>

to d) we observe that, while the wave propagates along the length of the cilium, the structure deforms itself with an oscillatory motion. For example, in Fig. 2 c), the sign of the activity scenario changes along the direction of the active fibers and is zero at approximate mid-point of the length of the cilium. This results in a double curved deformation of the structure since two opposite contraction behaviors are induced by the internal activity: at the bottom of the cilium the medium stretches in the left and contracts in the right, whereas at the top it contracts in the left and stretches in the right.

**A non symmetric scenario** In the study of the locomotion of microorganisms, it is well-known that in order to efficiently swim or propel the surrounding fluid, the movement of a cilium or a flagellum has to be non symmetric in time. This is the statement of Purcell's scallop Theorem [40] and this is due to reversibility properties of viscous fluids at low Reynold's number (which is the case we are interested in). Thus, to be able to model non symmetric internal activities is of primary importance for the study of active structures in a viscous fluid. In particular, we propose an activity scenario to mimic the deformations of cilia:

$$ \Sigma_a(t, (x_1, x_2)) = \frac{C_a}{L_c r_c} (\sigma_{a,1}(t, x_2) + \sigma_{a,2}(t, x_2)) (x_c - x_1), \quad \forall t \ge 0, \forall (x_1, x_2) \in \Omega_s, \quad (7) $$

where $\sigma_{a,1}(t, x_2)$ and $\sigma_{a,2}(t, x_2)$ are defined, for all $t \ge 0$ and for all $x_2$ in $[0, L_c]$ by

$$ \sigma_{a,1}(t, x_2) = (x_2 - L_c)^2 \sin\left(2\pi f_a \left(t - \frac{T_a}{4}\right)\right), $$

and

$$ \sigma_a(t, x_2) = \begin{cases} x_2 \left(x_2 - \frac{L_c}{2}\right) \cos\left(2\pi f_a \left(t - \frac{T_a}{4}\right)\right), & \text{if } t - \left\lfloor \frac{t}{T_a} \right\rfloor T_a \ge \frac{1}{2}, \\ 0, & \text{otherwise.} \end{cases} $$

As before, $f_a$ and $C_a$ still denote the beating frequency and the intensity of the internal activity. The parameters $T_a$ is the period of the beating, i.e. $T_a = \frac{1}{f_a}$. The results shown in Fig. 3 correspond to the set of parameters given in Table 3. They have been obtained by solving problem (1) with $\Sigma_a$ defined by (7) and using the finite element

<table><caption>Table 3: Set of parameters for the non symmetric scenario of activity.</caption><thead><tr><th>L<sub>c</sub> (μm)</th><th>r<sub>c</sub> (μm)</th><th>x<sub>c</sub> (μm)</th><th>C<sub>a</sub> (pN·μm)</th><th>f<sub>a</sub> (Hz)</th><th>E<sub>s</sub> (pN·μm<sup>-2</sup>)</th><th>ν<sub>s</sub></th></tr></thead><tbody><tr><td>6.5</td><td>0.2</td><td>0</td><td>3</td><td>10</td><td>10<sup>6</sup></td><td>0.49</td></tr></tbody></table>

method combined with a Newton's algorithm. We can observe that the activity scenario (7) is divided in two phases. During the first phase (Fig. 3 a) to Fig. 3 d)), taking place during the first half period of the beat, the function $\sigma_{a,2}$ is null and the cilium starts from a deformed configuration represented in Fig. 3 a). Then a bending
<figure><img src="image_3.png" /><figcaption>Figure 2: Deformation of an elongated structure subjected to the internal activity defined by (6) at different times: a) $t = 0$s, b) $t = 0.012$s, c) $t = 0.025$s and d) $t = 0.037$s. The activity scenario $\Sigma_a$ is represented with the set of parameters be given in Table 2.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Elongated structure subjected to the internal activity defined by (7) inducing cilia-like deformations at different time: a) $t = 0$s, b) $t = 0.015$s, c) $t = 0.035$s, d) $t = 0.050$s, e) $t = 0.050$s, f) $t = 0.065$s, g) $t = 0.065$s, h) $t = 0.085$s and e) $t = 0.1$s. The activity scenario $\Sigma_a$ is represented with the set of parameters be given in Table 3.</figcaption></figure>

deformation occurs to the right with the intensity being the greater close to the bottom of the cilium because of the term $(x_2 - L_c)^2$. At the end of the first phase the cilium is in the configuration represented in Fig. 3. During the second phase (Fig. 3 h) to Fig. 3 e)), the function $\sigma_{a,2}$ becomes non null and the deformation of the cilium is induced by the sum of both contributions of $\sigma_{a,1}$ and $\sigma_{a,2}$. In this phase, the bending occurs to the left. In particular the term $x_2(x_2 - \frac{L_c}{2})$ induces a double curved deformation that we can observe in Fig. 3 f).

**Remark.** The activity scenario (7) is used in Section 4 to perform numerical simulations of active structures beating in a viscous fluid.

**A scenario in three space dimensions** The framework developed in the present section is independent of the space dimension and can easily be adapted to model three-dimensional active structures. Let $L_c$, $r_c$ and $x_c$ still denote the length, the radius and the mean position on the $x_1$-axis of the cilium-like structure. Moreover, we denote by $y_c$ the mean position of the structure on the $x_2$-axis. In three space dimensions, $\Omega_s$ is the cylinder in $\mathbb{R}^3$ defined by

$$ \Omega_s = \{(x_1, x_2, x_3) \in \mathbb{R}^3; (x_c - x_1)^2 + (y_c - x_2)^2 \le r_c^2, 0 \le x_3 \le L_c\}. $$

Then, we consider the following activity scenario defined for all time $t \ge 0$ and for all $x = (x_1, x_2, x_3)$ in $\Omega_s$ by

$$ \Sigma_a(t, x) = \frac{C_a}{L_c \pi r_c^2} \sin(2\pi f_a (x_3 - v_a t)) (\cos(3\pi f_a t)(x_1 - x_c) + \sin(3\pi f_a t)(x_2 - y_c)), \quad (8) $$

where, as before, $f_a$ is the beating frequency, $C_a$ is the intensity of the activity and $v_a$ is the velocity of the wave that propagates along the length of the cilia. The deformations produced by the previous activity scenario are presented in Fig. 4, with the set of parameters be given in Table 4. It results in a twirling movement of the structure.

<table><caption>Table 4: Set of parameters for the twirling scenario of activity in three space dimensions.</caption><thead><tr><th>L<sub>c</sub> (μm)</th><th>r<sub>c</sub> (μm)</th><th>(x<sub>c</sub>, y<sub>c</sub>) (μm)</th><th>C<sub>a</sub> (pN · μm<sup>3</sup>)</th><th>f<sub>a</sub> (μm<sup>-1</sup>)</th><th>E<sub>s</sub> (pN · μm<sup>-2</sup>)</th><th>ν<sub>s</sub></th><th>v<sub>a</sub> (μm · s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>6.5</td><td>0.2</td><td>(0,0)</td><td>0.01</td><td>0.5</td><td>10<sup>6</sup></td><td>0.49</td><td>1</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 4: Elongated three-dimensional structure subjected to the internal activity defined by (8) inducing twirling at different time: a) $t = 0$s, b) $t = 0.052$s, c) $t = 0.082$s, d) $t = 0.107$s and e) $t = 0.128$s. The activity scenario $\Sigma_a$ in represented with the set of parameters be given in Table 4.</figcaption></figure>

# 3 A coupled fluid-structure interaction problem with cilia-like structures

In this section, we introduce a coupled fluid-structure interaction problem involving active structures, modeled by the active elasticity equations, and a Newtonian viscous incompressible homogeneous fluid, modeled by the Stokes equations. The choice of the Stokes equations is justified by the motivation of our study. Indeed, we recall that we aim at modeling the interaction between a biological fluid (essentially water) and microscopic ciliated organisms. Because of their microscopic size and their slow velocity, the Reynolds number associated to this system is of the order of $10^{-6} - 10^{-4}$ according to the microorganisms under consideration. Thus, the inertial effects of the system can be neglected and we consider non-inertial equations for both the fluid and the structures.

## 3.1 Description of the model

The problem presented in this section is intended to be used for numerical simulations. For that purpose we consider specific hypotheses on the domains to be meshed. Let $\Omega$ be a $n$-dimensional box ($n = 2$ or 3), i.e. the domain defined by $\Omega = \prod_{i=1}^n [0, L_i]$, where the $(L_i)_{i=1,...,n}$ are the dimensions of the box in every directions. The domain $\Omega$ is divided in two subsets, $\Omega = \Omega_s \cup \Omega_f$ such that $\Omega_f \cap \Omega_s = \emptyset$, where $\Omega_s$ is the solid domain, defined as the union of a finite number of cylindrical domains which represent cilia-like structures in their reference configuration, and $\Omega_f$ is the fluid domain in the reference configuration. The fluid-structure interface between $\Omega_s$ and $\Omega_f$ is denoted by $\Gamma$ and we define the remaining boundaries of $\Omega_s$ by $\Gamma_s = \partial\Omega_s \setminus \Gamma$, where homogeneous Dirichlet boundary conditions are applied. The boundary of $\Omega_f$ is divided in four parts, where different boundary conditions are applied. The first one is the fluid-structure interface $\Gamma$. Then, on the bottom boundary $\Gamma_{f,D} = \{x = (x_1, \dots, x_n) \in \Omega; x_n = 0\}$ we consider homogeneous Dirichlet boundary conditions, while on the top boundary $\Gamma_{f,S} = \{x = (x_1, \dots, x_n) \in \Omega; x_n = L_n\}$ we consider slip boundary conditions. Finally, on the lateral boundaries $\Gamma_{f,N} = \partial\Omega_f \setminus (\Gamma \cup \Gamma_{f,D} \cup \Gamma_{f,S})$, homogeneous Neumann boundary conditions are applied. Moreover, we denote by $\Gamma_f$ the union of $\Gamma_{f,D}$, $\Gamma_{f,S}$ and $\Gamma_{f,N}$:

$$ \Gamma_f = \Gamma_{f,D} \cup \Gamma_{f,S} \cup \Gamma_{f,N}. \qquad (9) $$

The domain $\Omega_s$ is filled with an active elastic medium which satisfies the Saint Venant-Kirchhoff law and whose activity is taken into account with the continuum mechanics framework developed in Section 2. Given a time dependent body force $f_s$ and a time dependent internal activity $\Sigma^*$, the problem of elasticity we are considering is to find the displacement $d_s$ of the solid medium such that, for all $t \ge 0$,

$$ \begin{aligned} -\mathrm{div}((I + \nabla d_s(t))(\Sigma_s(d_s(t)) - \Sigma^*(t))) &= f_s(t), & \text{in } \Omega_s, \\ d_s(t) &= 0, & \text{on } \Gamma_s, \end{aligned} \qquad (10) $$
where $\Sigma_s(d_s(t))$ is defined by (2).

The solid problem is always solved in the reference configuration of the structure $\Omega_s$, i.e. in Lagrangian coordinates, whereas the fluid problem is usually set up in Eulerian coordinates, i.e. in the current (deformed) configuration. This deformed configuration at time $t \ge 0$ depends on the displacement of the solid medium on the interface $\Gamma$ by means of a transformation denoted by $\Phi(d_s(t))$, which satisfies

$$ \Phi(d_s(t)) = \mathcal{I} + d_s(t), \quad \text{on } \Gamma. $$

The mapping $\mathcal{I}$ denotes the identity mapping in $\mathbb{R}^n$. Moreover, let $\gamma_\Gamma$ be the trace operator from $\Omega_s$ onto $\Gamma$ and $\mathcal{R}$ be a lifting operator from $\Gamma$ into $\Omega_f$ (in spaces made precise later on), then the transformation $\Phi(d_s(t))$ is defined in the whole domain $\Omega_f$ by

$$ \Phi(d_s(t)) = \mathcal{I} + \mathcal{R}(\gamma_\Gamma(d_s(t))), \quad \text{in } \Omega_f, $$

such that $\mathcal{R}(\gamma_\Gamma(d_s(t)))$ is null on the exterior boundary $\Gamma_f$, defined by (9). Thus, the transformation $\Phi(d_s(t))$ maps the reference fluid domain $\Omega_f$ to the deformed fluid domain at time $t$, $\Phi(d_s(t))(\Omega_f)$. In practice, the construction of the lifting operator $\mathcal{R}$ is done by solving an elliptic problem (typically a Laplace equation or the equations of linearized elasticity) on $\Omega_f$ with homogeneous Dirichlet boundary conditions on $\Gamma_f$ and applying $\gamma_\Gamma(d_s(t))$ as a Dirichlet boundary condition on the interface $\Gamma$.

At time $t$, the domain $\Phi(d_s(t))(\Omega_f)$ is filled with a Newtonian viscous incompressible homogeneous fluid, whose viscosity, $\mu_f$, is positive. Given a time dependent body force $f_f$, the velocity of the fluid $u_f$ and the pressure of the fluid $p_f$ satisfy, for all $t \ge 0$, the Stokes equations with mixed Dirichlet, Neumann and slip boundary conditions:

$$ \begin{aligned} -\operatorname{div}(\sigma_f(u_f(t), p_f(t))) &= f_f(t), & \text{in } \Phi(d_s(t))(\Omega_f), \\ \operatorname{div}(u_f(t)) &= 0, & \text{in } \Phi(d_s(t))(\Omega_f), \\ u_f(t) &= 0, & \text{on } \Gamma_{f,D}, \\ u_f(t) \cdot n_f(t) &= 0, & \text{on } \Gamma_{f,S}, \\ (\sigma_f(u_f(t), p_f(t))n_f(t)) \cdot \tau_f(t) &= 0, & \text{on } \Gamma_{f,S}, \\ \sigma_f(u_f(t), p_f(t))n_f(t) &= 0, & \text{on } \Gamma_{f,N}. \end{aligned} \quad (11) $$

The tensor $\sigma_f(u_f(t), p_f(t))$ is the fluid stress tensor defined by

$$ \sigma_f(u_f(t), p_f(t)) = 2\mu_f D(u_f(t)) - p_f(t)I, $$

where $D(u_f(t)) = \frac{1}{2}(\nabla u_f(t) + \nabla u_f(t)^T)$ is the symmetric gradient of $u_f(t)$ and $I$ denotes the identity matrix of $\mathbb{R}^n$.

Equations (10) and (11) are completed by the usual fluid-structure coupling conditions on the interface $\Gamma$, namely the equality of the velocities and the continuity of the normal components of stress tensors. Because these conditions are written in the reference configuration, we introduce the velocity of the fluid written in the fluid reference configuration, $w_f$, and the pressure of the fluid written in the fluid reference configuration, $q_f$, defined by

$$ w_f(t, \cdot) = u_f(t, \Phi(d_s(t, \cdot))) \quad \text{and} \quad q_f(t, \cdot) = p_f(t, \Phi(d_s(t, \cdot))), \quad \text{in } \Omega_f. \quad (12) $$

Moreover, the fluid stress tensor at time $t$ written in the fluid reference configuration, denoted by $\Pi_f(w_f(t), q_f(t))$, is defined by

$$ \Pi_f(w_f(t), q_f(t)) = \mu_f (\nabla w_f(t) F(d_s(t)) + (\nabla(\Phi(d_s(t))))^{-T} \nabla w_f(t)^T G(d_s(t))) - q_f(t) G(d_s(t)), $$

where $F(d_s(t))$ and $G(d_s(t))$ are the following matrices:

$$ \begin{aligned} F(d_s(t)) &= (\nabla(\Phi(d_s(t))))^{-1}\mathrm{cov}(\nabla(\Phi(d_s(t)))), \\ G(d_s(t)) &= \mathrm{cov}(\nabla(\Phi(d_s(t)))). \end{aligned} \quad (13) $$

Then, for all $t \ge 0$, the coupling conditions on $\Gamma$ write

$$ \begin{aligned} \frac{\partial d_s(t)}{\partial t} &= w_f(t), & \text{on } \Gamma, \\ (I + \nabla d_s(t))(\Sigma_s(d_s(t)) - \Sigma^*(t))n_s &= \Pi_f(w_f(t), q_f(t))n_s, & \text{on } \Gamma. \end{aligned} \quad (14) $$

The condition on the continuity of the velocity through the interface $\Gamma$, considered in (14), requires an initial condition concerning the initial displacement of the structure on $\Gamma$. Then, we suppose that all cilia-like structures are in their reference configuration at time $t = 0$, i.e. that $d_s(0) = 0$ on $\Gamma$.
**Remark.** In problem (11), the current fluid domain $\Phi(d_s(t))(\Omega_f)$ is an unknown of the problem, since it depends on the displacement of the structure. In the present article, we discretize in time equations (10), (11) and (14), in order to make the configuration of the current fluid domain at a given time step entirely determined by the displacement of the structure at the previous time step. This is the purpose of the next subsection.

**Remark.** In definition (13), the matrix $F(d_s(t))$, that appears in the expression of the fluid stress tensor written in the fluid reference configuration, is well-defined if the mapping $\Phi(d_s(t))$ is, for example, a $C^1$-diffeomorphism. This kind of regularity is not intended to be proved in the present article but it is studied in [41]. In the mean time, we can justify our approach by recalling that our purpose here is to compute numerical simulations based on the model constructed within the present section and that, in a context of numerical simulations, the deformation at a given time step is easily invertible. For this reason, we assume in the following that the deformations are regular enough to ensure the definition of the problem.

## 3.2 The discrete-in-time fluid-structure interaction problem

Let us introduce a discretization of $\mathbb{R}^+$ for the time variable $t$. Let $\delta t > 0$ be a constant time step. We construct a sequence $(t_k)_{k \in \mathbb{N}}$ such that

$$ 
\begin{aligned}
t_0 &= 0, \\
t_{k+1} &= t_k + \delta t, \forall k \ge 0,
\end{aligned}
$$

and, for all $k \ge 0$, we define the time-discretizations of the displacement of the structure $d_s^k$, the velocity of the fluid $u_f^k$ and the pressure of the fluid $p_f^k$ by

$$ d_s^k \approx d_s(t_k), \quad u_f^k \approx u_f(t_k), \quad \text{and} \quad p_f^k \approx p_f(t_k). $$

Similarly, the time-discretization of the body forces $f_f^k$ and $f_s^k$ and the activity $\Sigma_k^*$ are defined, for all $k \ge 0$ by

$$ f_f^k = f_f(t_k), \quad f_s^k = f_s(t_k) \quad \text{and} \quad \Sigma_k^* = \Sigma^*(t_k) $$

Moreover, the first coupling condition in (14) is discretized using the implicit Euler scheme:

$$ 
\begin{aligned}
d_s^0 &= 0, \\
d_s^{k+1} &= d_s^k + \delta t w_f^{k+1}, \forall k \ge 0.
\end{aligned}
$$

For $k > 0$, the deformed fluid domain at time $t_k$, $\Phi(d_s(t_k))(\Omega_f)$, is denoted by $\Omega_f^k$ and depends on the displacement of the structure at the previous time step by

$$ \Omega_f^k = \Phi(d_s^{k-1})(\Omega_f). $$

Similarly, the fluid-structure interface at time $t_k$ writes

$$ \Gamma^k = \Phi(d_s^{k-1})(\Gamma), $$

and the remaining fluid boundaries in the deformed configuration do not change, i.e.

$$ 
\begin{aligned}
\Gamma_{f,D}^k &= \Gamma_{f,D}, \\
\Gamma_{f,S}^k &= \Gamma_{f,S}, \\
\Gamma_{f,N}^k &= \Gamma_{f,N}.
\end{aligned}
$$

For $k = 0$, we suppose by convention that $d_s^{0-1} = 0$ on $\Gamma$, such that $\Omega_f^0 = \Omega_f$ and $\Gamma^0 = \Gamma$.

Then the discrete-in-time fluid-structure interaction problem that we consider is to find, for all $k \ge 0$, the displacement of the structure $d_s^k$, the velocity of the fluid $u_f^k$ and the pressure of the fluid $p_f^k$ which satisfy

$$ 
\begin{aligned}
-\mathrm{div}((I + \nabla d_s^k)(\Sigma_s(d_s^k) - \Sigma_k^*)) &= f_s^k, && \text{in } \Omega_s, \\
d_s^k &= 0, && \text{on } \Gamma_s,
\end{aligned}
\quad (15) $$

$$ 
\begin{aligned}
-\mathrm{div}(\sigma_f(u_f^k, p_f^k)) &= f_f^k, && \text{in } \Omega_f^k, \\
\mathrm{div}(u_f^k) &= 0, && \text{in } \Omega_f^k, \\
u_f^k &= 0, && \text{on } \Gamma_{f,D}, \\
u_f^k \cdot n_f^k &= 0, && \text{on } \Gamma_{f,S}, \\
(\sigma_f(u_f^k, p_f^k)n_f^k) \cdot \tau_f^k &= 0, && \text{on } \Gamma_{f,S}, \\
\sigma_f(u_f^k, p_f^k)n_f^k &= 0, && \text{on } \Gamma_{f,N},
\end{aligned}
\quad (16) $$
$$ 
\begin{aligned}
d_s^k &= d_s^{k-1} + \delta t w_f^k, & \text{on } \Gamma, \\
(I + \nabla d_s^k)(\Sigma_s(d_s^k) - \Sigma_k^*) n_s &= \Pi_f(w_f^k, q_f^k) n_s, & \text{on } \Gamma,
\end{aligned} 
\quad (17) $$

where

$$ 
\begin{aligned}
w_f^k(\cdot) &= u_f^k(\Phi(d_s^{k-1})(\cdot)), & \text{in } \Omega_f, \\
q_f^k(\cdot) &= p_f^k(\Phi(d_s^{k-1})(\cdot)), & \text{in } \Omega_f.
\end{aligned}
 $$

In Section 4 we present the numerical method used for the simulation of active structures in a viscous fluid. The method is based on a saddle-point formulation for the discrete-in-time coupled system of equations (15), (16) and (17), that is well-suited for the numerical simulation with standard finite element techniques and solvers, and that will be introduced later on. To that aim, we first study the weak formulation of problem (15), (16) and (17) in Subsection 3.3, then, the constraint on the continuity of velocities is expressed with a Lagrange multiplier and the saddle-point problem is introduced in Subsection 3.4. For the sake of simplicity, the well-posedness of both problems is studied in the linearized case, i.e. the case where the active elasticity equations are linearized around the equilibrium. To that aim, we introduce the linearized active elasticity problem, which consists in finding the displacement of the structure $d_s^k$ which satisfies

$$ 
\begin{aligned}
-\mathrm{div}(\sigma_s(d_s^k) - \nabla d_s^k \Sigma_k^*) &= f_s^k + \mathrm{div}(\Sigma^*), & \text{in } \Omega_s, \\
d_s^k &= 0, & \text{on } \Gamma_s,
\end{aligned} 
\quad (18) $$

where $\sigma_s(d_s^k)$ is the linearized stress tensor of the structure around the equilibrium at time $t_k$ which writes

$$ \sigma_s(d_s^k) = 2\mu_f D(d_s^k) + \lambda_s \mathrm{div}(d_s^k) I. \quad (19) $$

Coupling conditions (17) are also linearized near the equilibrium and become

$$ 
\begin{aligned}
d_s^k &= d_s^{k-1} + \delta t w_f^k, & \text{on } \Gamma, \\
(\sigma_s(d_s^k) - \nabla d_s^k \Sigma_k^*) n_s &= \Pi_f(w_f^k, q_f^k) n_s + \Sigma^* n_s, & \text{on } \Gamma.
\end{aligned} 
\quad (20) $$

Thus, the linearized coupled fluid-structure problem we are interested in and that will be studied in the next subsection is the system of equations (16), (18) and (20).

## 3.3 Existence and uniqueness of a weak solution to the linearized fluid-structure interaction problem

In this subsection, we define and study the weak formulation of the linearized fluid-structure problem (16), (18) and (20). First, let us introduce the following function spaces, for all $k \ge 0$:

$$ 
\begin{aligned}
V_f^k &= \{ v \in (H^1(\Omega_f^k))^n; \gamma_{\Gamma_{f,D}}(v) = 0, \gamma_{\Gamma_{f,S}}(v) \cdot n_f^k = 0 \}, \\
V_s &= \{ v \in (H^1(\Omega_s))^n; \gamma_{\Gamma_s}(v) = 0 \}. \\
W_u &= \{ (v_f, v_s) \in V_f^k \times V_s; \gamma_{\Gamma}(v_f \circ \Phi(d_s^{k-1})) = \gamma_{\Gamma}(v_s) \}, \\
W_d &= \{ (v_f, d_s) \in V_f^k \times V_s; \delta t \gamma_{\Gamma}(v_f \circ \Phi(d_s^{k-1})) + d_s^{k-1} = \gamma_{\Gamma}(d_s) \}.
\end{aligned} 
\quad (21) $$

Then, the weak formulation of problem (16), (18) and (20) is given by the following lemma.

**Lemma 1.** The weak formulation of problem (16), (18) and (20) is defined by

$$ 
\left\{
\begin{array}{l}
\displaystyle \text{find } (u_f^k, d_s^k) \in W_d \text{ and } p_f^k \in L^2(\Omega_f^k) \text{ such that} \\
\displaystyle \int_{\Omega_f^k} \sigma_f(u_f^k, p_f^k) : \nabla v_f + \int_{\Omega_s} (\sigma_s(d_s^k) - \nabla d_s^k \Sigma_k^*) : \nabla v_s \\
\displaystyle = \int_{\Omega_f^k} f_f^k \cdot v_f + \int_{\Omega_s} f_s^k \cdot v_s - \int_{\Omega_s} \Sigma_k^* : \nabla v_s, \quad \forall (v_f, v_s) \in W_u, \\
\displaystyle \int_{\Omega_f^k} q_f \mathrm{div}(u_f^k) = 0, \quad \forall q_f \in L^2(\Omega_f^k).
\end{array}
\right. 
\quad (22) $$

*Proof.* Let $(v_f, v_s)$ be a regular function in $W_u$ and suppose that $u_f^k$, $d_s^k$ and $p_f^k$ are sufficiently regular. We can formally multiply the first equation in (16) by $v_f$ and the first equation in (18) by $v_s$ and integrate respectively over $\Omega_f^k$ and $\Omega_s$. After an integration by part, we obtain

$$ \int_{\Omega_f^k} \sigma_f(u_f^k, p_f^k) : \nabla v_f - \int_{\Gamma^k \cup \Gamma_{f,S}} (\sigma_f(u_f^k, p_f^k) n_f^k) \cdot v_f = \int_{\Omega_f^k} f_f^k \cdot v_f, \quad (23) $$
and

$$ \begin{aligned} & \int_{\Omega_s} (\sigma_s(d_s^k) - \nabla d_s^k \Sigma_k^*) : \nabla v_s - \int_{\Gamma} ((\sigma_s(d_s^k) - \nabla d_s^k \Sigma_k^*) n_s) \cdot v_s \\ & = \int_{\Omega_s} f_s^k \cdot v_s - \int_{\Omega_s} \Sigma_k^* : \nabla v_s + \int_{\Gamma} (\Sigma_k^* n_s) \cdot v_s. \end{aligned} \tag{24} $$

Moreover, decomposing the normal stress $\sigma_f(u_f^k, p_f^k)n_f^k$ on $\Gamma_{f,S}$ in its normal and tangential components and using the slip boundary conditions on $\Gamma_{f,S}$ for the fluid problem, it follows that

$$ \int_{\Gamma_{f,S}} (\sigma_f(u_f^k, p_f^k)n_f^k) \cdot v_f = \int_{\Gamma_{f,S}} (((\sigma_f(u_f^k, p_f^k)n_f^k) \cdot n_f^k) + ((\sigma_f(u_f^k, p_f^k)n_f^k) \cdot \tau_f^k) \cdot v_f = 0. $$

Then, after a change of variables and making use of the second coupling condition in (20), we have

$$ \begin{aligned} & \int_{\Gamma^k} (\sigma_f(u_f^k, p_f^k)n_f^k) \cdot v_f = \int_{\Gamma} (\Pi_f(w_f^k, q_f^k)n_f) \cdot (v_f \circ \Phi(d_s^{k-1})), \\ & = \int_{\Gamma} ((\sigma_s(d_s^k) - \nabla d_s^k \Sigma_k^*)n_s) \cdot v_s - \int_{\Gamma} (\Sigma_k^* n_s) \cdot v_s. \end{aligned} $$

Now, summing equation (23) and equation (24) it comes

$$ \int_{\Omega_f^k} \sigma_f(u_f^k, p_f^k) : \nabla v_f + \int_{\Omega_s} (\sigma_s(d_s^k) - \nabla d_s^k \Sigma_k^*) : \nabla v_s = \int_{\Omega_f^k} f_f^k \cdot v_f + \int_{\Omega_s} f_s^k \cdot v_s - \int_{\Omega_s} \Sigma_k^* : \nabla v_s. $$

Similarly, let $q_f$ be in $L^2(\Omega_f^k)$. Formally, we multiply the second equation in (16) by $q_f$ and integrate over $\Omega_f^k$. We obtain

$$ \int_{\Omega_f^k} q_f \mathrm{div}(u_f^k) = 0. $$

□

To prove that problem (22) is well-posed, we perform a change in variables in the displacement $d_s^k$ in order to work on a velocity-velocity formulation of the fluid-structure problem (instead of a velocity-displacement formulation). Let us introduce the discrete-in-time velocity of the structure at time $t_k$, $u_s^k$, defined by

$$ u_s^k = \frac{1}{\delta t} (d_s^k - d_s^{k-1}). \tag{25} $$

Thus, the previous weak problem is equivalent to the problem where $d_s^k$ has been replaced by $\delta t u_f^k + d_s^{k-1}$:

$$ \left\{ \begin{array}{l} \text{find } (u_f^k, u_s^k) \in W_u \text{ and } p_f^k \in L^2(\Omega_f^k) \text{ such that} \\ a^k((u_f^k, u_s^k), (v_f, v_s)) - (B(v_f, v_s), p_f^k)_{L^2(\Omega_f^k)} = L^k(v_f, v_s), \quad (v_f, v_s) \in W_u, \\ (B(u_f^k, u_s^k), q)_{L^2(\Omega_f^k)} = 0, \quad \forall q \in L^2(\Omega_f^k), \end{array} \right. \tag{26} $$

where $(\cdot, \cdot)_{L^2(\Omega_f^k)}$ denotes the scalar product in $L^2(\Omega_f^k)$ and $a^k, L^k$ and $B$ are defined by:

$$ - \text{ for all } (u_f, u_s), (v_f, v_s) \in V_f^k \times V_s, $$

$$ a^k((u_f, u_s), (v_f, v_s)) = 2\mu_f \int_{\Omega_f^k} D(u_f) : D(v_f) + \delta_t \int_{\Omega_s} \sigma_s(u_s) : \nabla v_s - \delta_t \int_{\Omega_s} (\nabla u_s \Sigma_k^*) : \nabla v_s, $$

$$ L^k(v_f, v_s) = \int_{\Omega_f^k} f_f^k \cdot v_f + \int_{\Omega_s} f_s^k \cdot v_s - \int_{\Omega_s} \Sigma_k^* : \nabla v_s - \int_{\Omega_s} (\sigma_s(d_s^{k-1}) - \nabla d_s^{k-1} \Sigma_k^*) : \nabla v_s, $$

$$ - \text{ for all } (v_f, v_s) \in W_u, $$

$$ B(v_f, v_s) = \mathrm{div}(v_f). $$
**Theorem 3.1.** Let $k \ge 0$ and suppose that the force $f_f^k$ belongs to $(L^2(\Omega_f^k))^n$, the force $f_s^k$ belongs to $(L^2(\Omega_s))^n$, the displacement $d_s^{k-1}$ belongs to $V_s$ and the activity tensor $\Sigma_k^*$ belongs to $(L^\infty(\Omega_s))^{n \times n}$. There exists a constant $C(n, \Omega_s)$ which only depends on the dimension $n$ and on the domain $\Omega_s$ such that, if $\Sigma_k^*$ satisfies

$$ \|\Sigma^*\|_{L^\infty(\Omega_s)} < C(\Omega_s)\mu_s, \qquad (27) $$

then there exists a unique solution to problem (26).

*Proof.* Let us show that $a$ is a continuous bilinear coercive form on $(V_f^k \times V_s)^2$, that $L$ is a continuous linear form on $V_f^k \times V_s$ and that $B$ is a continuous linear surjective operator from $W_u$ to $L^2(\Omega_f^k)$. Because $\Sigma^*$ belongs to $(L^\infty(\Omega_s))^{n \times n}$, it is clear that $a$ is a continuous bilinear form on $(V_f^k \times V_s)^2$ and that $L$ is a continuous linear form on $V_f^k \times V_s$. In addition, the operator $B$ is also continuous and linear from $W_u$ to $L^2(\Omega_f^k)$. Now we show that $a^k$ is coercive under condition (27). Using the $L^\infty$-regularity of $\Sigma^*$ and Korn's inequality we have, for all $(u_f, u_s)$ in $V_f^k \times V_s$,

$$ a^k((u_f, u_s), (u_f, u_s)) = 2\mu_f \|D(u_f)\|_{L^2(\Omega_f^k)}^2 + 2\mu_s \delta t \|D(u_s)\|_{L^2(\Omega_s)}^2 + \lambda_s \delta t \|\mathrm{div}(u_s)\|_{L^2(\Omega_s)}^2 - \delta t \int_{\Omega_s} (\nabla u_s \Sigma_k^*) : \nabla u_s. $$

Yet, the last integral writes

$$ \begin{aligned} \int_{\Omega_s} (\nabla u_s \Sigma_k^*) : \nabla u_s &= \sum_{i,j=1}^n \int_{\Omega_s} (\nabla u_s \Sigma_k^*)_{ij} (\nabla u_s)_{ij}, \\ &= \sum_{i,j=1}^n \int_{\Omega_s} \left( \sum_{k=1}^n (\nabla u_s)_{ik} (\Sigma_k^*)_{kj} \right) (\nabla u_s)_{ij}, \end{aligned} $$

and using the $L^\infty$-regularity of $\Sigma^*$ it follows that

$$ \begin{aligned} \int_{\Omega_s} (\nabla u_s \Sigma_k^*) : \nabla u_s &\le \|\Sigma_k^*\|_{L^\infty(\Omega_s)} \sum_{i,j=1}^n \int_{\Omega_s} \left(\sum_{k=1}^n (\nabla u_s)_{ik}\right) (\nabla u_s)_{ij}, \\ &\le \|\Sigma_k^*\|_{L^\infty(\Omega_s)} \sum_{i=1}^n \int_{\Omega_s} \left(\sum_{j=1}^n (\nabla u_s)_{ij}\right)^2. \end{aligned} $$

Moreover, using Young's inequality, it comes

$$ \begin{aligned} \int_{\Omega_s} (\nabla u_s \Sigma_k^*) : \nabla u_s &\le \|\Sigma_k^*\|_{L^\infty(\Omega_s)} \sum_{i=1}^n \int_{\Omega_s} n \sum_{j=1}^n (\nabla u_s)_{ij}^2, \\ &\le n \|\Sigma_k^*\|_{L^\infty(\Omega_s)} \|\nabla u_s\|_{L^2(\Omega_s)}^2. \end{aligned} $$

Now using Poincaré and Korn inequalities, we have

$$ \begin{aligned} a^k((u_f, u_s), (u_f, u_s)) &\ge 2\mu_f C_K(\Omega_f^k) \|u_f\|_{H^1(\Omega_f^k)}^2 + 2\mu_s \delta t C_K(\Omega_s) \|u_s\|_{H^1(\Omega_s)}^2 - n\delta t C_P(\Omega_s) \|\Sigma^*\|_{L^\infty(\Omega_s)} \|u_s\|_{H^1(\Omega_s)}^2, \\ &\ge 2\mu_f C_K(\Omega_f^k) \|u_f\|_{H^1(\Omega_f^k)}^2 + (2\mu_s C_K(\Omega_s) - nC_P(\Omega_s) \|\Sigma^*\|_{L^\infty(\Omega_s)}) \|u_s\|_{H^1(\Omega_s)}^2, \end{aligned} $$

where $C_K(\Omega_f^k)$, $C_K(\Omega_s)$ and $C_P(\Omega_s)$ are positive constants which only depend on the domains $\Omega_s^k$ and $\Omega_s$. So, denoting

$$ C(\Omega_s) = \frac{2C_K(\Omega_s)}{nC_P(\Omega_s)} $$

and under condition (27), the continuous bilinear form $a^k$ is coercive on $V_f^k \times V_s$.

To conclude, it remains to show that the operator $B$ is surjective from $W_u$ to $L^2(\Omega_f^k)$. Indeed, let $q$ be in $L^2(\Omega_f^k)$. We extend function $q$ in the whole space $L_0^2(\Omega)$ considering the extension operator $E_p$, defined by

$$ E_p q = \begin{cases} q & \text{in } \Omega_f^k, \\ \frac{-1}{|\Omega \setminus \Omega_f^k|} \int_{\Omega_f^k} q & \text{in } \Omega \setminus \Omega_f^k, \end{cases} $$
where $|\Omega \setminus \Omega_f^k|$ is the volume of the subset $\Omega \setminus \Omega_f^k$. Since $E_p q$ belongs to $L_0^2(\Omega)$, Bogovskii's result, see [42], ensures the existence of a function $\tilde{u}$ in $(H_0^1(\Omega))^n$ such that $\text{div}(\tilde{u}) = E_p q$. Then, we define $v_f = \tilde{u}_{|\Omega_f^k}$ and $v_s = \tilde{u}_{|\Omega_s^k}$ and it follows that the couple $(v_f, v_s \circ \Phi(d_s^{k-1}))$ belongs to $W_u$ and satisfies

$$B(v_f, v_s \circ \Phi(d_s^{k-1})) = \text{div}(v_f) = q.$$

This proves the surjectivity of the operator $B$. As a conclusion, according to [43], problem (26) admits a unique solution $(u_f^k, u_s^k, p_f^k)$.

**Remark.** As problems (22) and (26) are equivalent, Theorem 3.1 also applies to problem (22).

## 3.4 A saddle-point formulation for the fluid-active structure interaction problem

Even though problems (22) and (26) are well-posed, they are rather complicated to solve in the context of numerical simulations with standard finite element techniques. The main reason is that the fluid equations (16) and the solid equations (18) are written in two different configurations, with transmission conditions on the fluid-structure interface. For the direct simulation using the finite element method, this means that one is supposed to construct a basis of finite element functions that should approximate the whole space $W_u \times L^2(\Omega_f^k)$, which does not enable to use standard finite element solvers.

Another strategy is to use an iterative method and solve both problems separately. It has the advantage to make use of existing solvers for both problems, but it requires a particular method to treat the coupling conditions on the fluid-structure interface. In the present article, in the context of direct numerical simulations, we consider a fitted-mesh method based on a Lagrangian multiplier to impose the continuity of the velocity through the fluid-structure interface. Actually this provides an easy to implement method. For that purpose, the present subsection is dedicated to the introduction of a different formulation of problem (26), where the constraint of equality of the fluid and solid velocities on $\Gamma$, which appears in the function space $W_u$, is treated by duality and enforced with a Lagrangian multiplier.

Let us introduce the constraint operator $K$, defined by

$$K: \begin{array}{rcl} V_f^k \times V_s & \rightarrow & L^2(\Omega_f^k) \times \Upsilon \\ (v_f, v_s) & \mapsto & (\text{div}(v_f), \gamma_\Gamma(v_f \circ \Phi(d_s^{k-1}))) - \gamma_\Gamma(v_s), \end{array}$$

where the space $\Upsilon$, defined by

$$\Upsilon = (H_{00}^{1/2}(\Gamma))^n,$$

denotes the image of $V_s$ (resp. $V_f$) by the trace operator on $\Gamma$. In other words, $\Upsilon$ is the space of functions in $(H^{1/2}(\Gamma))^n$ whose extension by zero on $\Gamma_s$ (resp. $\Gamma_f$) belongs to $(H^{1/2}(\partial\Omega_s))^n$ (resp. $(H^{1/2}(\partial\Omega_f))^n$). Then, we define and study the well-posedness of the following (non-constrained) saddle-point problem:

$$\left\{ \begin{array}{l} \text{find } (u_f^k, u_s^k) \in V_f^k \times V_s, p_f^k \in L^2(\Omega_f^k) \text{ and } \lambda^k \text{ in } \Upsilon \text{ such that} \\ \qquad a^k((u_f^k, u_s^k), (v_f, v_s)) - (p_f^k, \lambda^k), K(v_f, v_s) \Big|_{L^2(\Omega_f^k) \times \Upsilon} = L^k(v_f, v_s), \quad \forall (v_f, v_s) \in V_f^k \times V_s, \\ \qquad (q_f, \mu), K(u_f^k, u_s^k) \Big|_{L^2(\Omega_f^k) \times \Upsilon} = 0, \quad \forall (q_f, \mu) \in L^2(\Omega_f^k) \times \Upsilon. \end{array} \right. \qquad (28)$$

**Theorem 3.2.** Let $k \ge 0$ and suppose that the force $f_f^k$ belongs to $(L^2(\Omega_f^k))^n$, the force $f_s^k$ belongs to $(L^2(\Omega_s))^n$, the displacement $d_s^{k-1}$ belongs to $V_s$ and the active tensor $\Sigma_k^*$ belongs to $(L^\infty(\Omega_s))^{n \times n}$. If $\Sigma_k^*$ satisfies condition (27), then there exists a unique solution to problem (28).

*Proof.* As before, the well-posedness of problem (28) is proved using standard results on saddle-point problems, see [43]. It has already been argued in the proof of Theorem 3.1 that $a^k$ is a continuous coercive bilinear form on $(V_f^k \times V_s)^2$ and that $L^k$ is a continuous linear form on $V_f^k \times V_s$. Moreover, the operator $K$ is clearly linear and continuous from $V_f^k \times V_s$ to $L^2(\Omega_f^k) \times \Upsilon$ due to the continuity of the divergence operator and the continuity of the trace operators from $V_f$ to $\Gamma$ and from $V_s$ to $\Gamma$.

Let us show that $K$ is surjective from $V_f^k \times V_s$ to $L^2(\Omega_f^k) \times \Upsilon$. Given a couple $(q, \mu)$ in $L^2(\Omega_f^k) \times \Upsilon$, we aim to find a couple $(v_f, v_s)$ in $V_f^k \times V_s$ such that $K(v_f, v_s) = (q, \mu)$. First, suppose that $v_s$ is known. Then, we construct $\hat{v}_f$ in $(H_{\Gamma_f}^1(\Omega_f))^n$ such that the trace of $\hat{v}_f$ on $\Gamma$ satisfies

$$\gamma_{\Gamma}(\hat{v}_f) = \mu + \gamma_{\Gamma}(v_s). \qquad (29)$$
This is due to the existence of a continuous linear lifting operator from $\Upsilon$ to $(H_{\Gamma_f}^1(\Omega_f))^n$, since $\Omega_f$ is a Lipschitz domain (see [44, app. B]), and because $\mu + \gamma_\Gamma(v_s)$ belongs to $\Upsilon$. Now, suppose that

$$ \int_{\Omega_f^k} q - \operatorname{div}(\hat{v}_f \circ \Phi^{-1}(d_s^{k-1})) = 0. \qquad (30) $$

Then, Bogovskii's result in [42] ensures that there exists a function $\tilde{v}_f$ in the space $(H_0^1(\Omega_f^k))^n$ such that

$$ \operatorname{div}(\tilde{v}_f) = q - \operatorname{div}(\hat{v}_f \circ \Phi^{-1}(d_s^{k-1})). $$

Defining $v_f = \hat{v}_f \circ \Phi^{-1}(d_s^{k-1}) + \tilde{v}_f$, then $v_f$ belongs to $(H_{\Gamma_f}^1(\Omega_f^k))^n$ (a subspace of $V_f^k$) and one has $\operatorname{div}(v_f) = q$. It remains to construct $v_s$ such that condition (30) holds true. In fact, using the Piola identity and equation (29), condition (30) becomes

$$ \int_{\Omega_f^k} q = \int_{\Omega_f^k} \operatorname{div}(\hat{v}_f \circ \Phi^{-1}(d_s^k)) = \int_{\Omega_f} \operatorname{div}(G(d_s^{k-1})^T \hat{v}_f) = \int_{\Gamma} (G(d_s^{k-1})^T \hat{v}_f) n_f = \int_{\Gamma} (G(d_s^{k-1})^T (\mu + v_s)) n_f. $$

Thus, we need to construct $v_s$ in $V_s$ such that

$$ \int_{\Gamma} (G(d_s^{k-1})^T v_s) n_f = \int_{\Omega_f^k} q - \int_{\Gamma} (G(d_s^{k-1})^T \mu) n_f. $$

Taking any function $\tilde{v}_s$ in $V_s$ such that

$$ \int_{\Gamma} (G(d_s^{k-1})^T \tilde{v}_s) n_f \neq 0, $$

we define

$$ v_s = \frac{\int_{\Omega_f^k} q - \int_{\Gamma} (G(d_s^{k-1})^T \mu) n_f}{\int_{\Gamma} (G(d_s^{k-1})^T \tilde{v}_s) n_f} = \tilde{v}_s. $$

Then $v_s$ belongs to $V_s$ and satisfies condition (30). As a consequence, the couple $(v_f, v_s)$ belongs to $V_f^k \times V_s$ and satisfies $K(v_f, v_s) = (q, \mu)$. Finally, the operator $K$ is surjective from $V_f^k \times V_s$ to $L^2(\Omega_f^k) \times \Upsilon$ and problem (28) admits a unique solution $(u_f^k, u_s^k, p_f^k, \lambda^k)$. $\square$

The main difference between problem (26) and problem (28) is that the function spaces involved in problem (28) are free of constraints, whereas the function space $W_u$ involved in problem (26) is not. In particular, recalling that $u_s^k$ is defined by (25), problem (28) is equivalent to the following problem:

$$ \left\{ \begin{array}{ll} \text{find } u_f^k \in V_f^k, p_f^k \in L^2(\Omega_f^k), d_s^k \in V_s \text{ and } \lambda^k \text{ in } \Upsilon \text{ such that} & \\ \left. \begin{array}{l} 2\mu_f \int_{\Omega_f^k} D(u_f^k) : D(v_f) - \int_{\Omega_f^k} p_f^k \operatorname{div}(v_f) = \int_{\Omega_f^k} f_f^k \cdot v_f + (\lambda^k, v_f \circ \Phi(d_s^{k-1}))_{\Upsilon}, \quad \forall v_f \in V_f^k, \\ \int_{\Omega_f^k} q_f \operatorname{div}(u_f^k) = 0, \quad \forall q_f \in L^2(\Omega_f^k), \\ \int_{\Omega_s} \sigma_s(d_s^k) : \nabla v_s - \int_{\Omega_s} (\nabla d_s \Sigma_k^*) : \nabla v_s = \int_{\Omega_s} f_s^k \cdot v_s - (\lambda^k, v_s)_{\Upsilon}, \quad \forall v_s \in V_s, \\ \left( \mu, \gamma_\Gamma(u_f^k \circ \Phi(d_s^{k-1})) - \frac{1}{\delta t}(\gamma_\Gamma(d_s^k) - \gamma_\Gamma(d_s^{k-1})) \right)_{\Upsilon} = 0, \quad \forall \mu \in \Upsilon. \end{array} \right. & (31) \end{array} \right. $$

# 4 Numerical simulations of active structures in a viscous fluid

For the direct simulation of fluid-structure problems, the saddle-point formulation (31) of the problem is particularly interesting since the resolution of the fluid and structure problems may take advantage of the introduction of the Lagrange multipliers $\lambda^k$ by using an iterative method, such as Uzawa's algorithm. This is the purpose of the next section.
## 4.1 Description of the method

Coming back to our initial problem, where the structure satisfies the (nonlinear) equations of elasticity, we readily adapt, for all $k \ge 0$, by analogy with problem (31), the saddle-point formulation of equations (15), (16) and (17) by

$$ 
\left\{
\begin{array}{l}
\text{find } u_f^k \in V_f^k, p_f^k \in L^2(\Omega_f^k), d_s^k \in V_s \text{ and } \lambda^k \text{ in } \Upsilon \text{ such that} \\
2\mu_f \int_{\Omega_f^k} D(u_f^k) : D(v_f) - \int_{\Omega_f^k} p_f^k \operatorname{div}(v_f) = \int_{\Omega_f^k} f_f^k \cdot v_f + (\lambda^k, v_f \circ \Phi(d_s^{k-1}))_{\Upsilon}, \quad \forall v_f \in V_f^k, \\
\qquad \int_{\Omega_f^k} q_f \operatorname{div}(u_f^k) = 0, \quad \forall q_f \in L^2(\Omega_f^k), \\
\int_{\Omega_s} (I + \nabla d_s^k)(\Sigma_s(d_s^k) - \Sigma_k^*) : \nabla v_s = \int_{\Omega_s} f_s^k \cdot v_s - (\lambda^k, v_s)_{\Upsilon}, \quad \forall v_s \in V_s, \\
\left( \mu, \gamma_{\Gamma}(u_f^k \circ \Phi(d_s^{k-1})) - \frac{1}{\delta t}(\gamma_{\Gamma}(d_s^k) - \gamma_{\Gamma}(d_s^{k-1})) \right)_{\Upsilon} = 0, \quad \forall \mu \in \Upsilon.
\end{array}
\right.
\quad (32) $$

Problem (32) is solved at each time $t_k$ using Uzawa's algorithm. It consists in constructing a sequence $(\lambda^{k,j})_{j \in \mathbb{N}}$ in the space $\Upsilon$ which converges, under assumption, to the Lagrange multiplier solution of problem (32). At each iteration $j$ of Uzawa's algorithm, the Lagrange multiplier $\lambda^{k,j}$ is known. Thus, the resolution of problem (32) reduces to the resolution of a Stokes problem and an elasticity problem, where $\lambda^{k,j}$ is seen as a Neumann boundary condition on the fluid-structure interface. In practice, solutions of both problems are approximated with the finite element method on conformal meshes and standard finite element solvers are used, since both variational problems are classical. In two space dimensions, the Stokes problem is solved in mixed formulation with a direct solver and using Mini elements. The nonlinear equations of elasticity are solved with a Newton solver and using $P_1$ Lagrange elements. Such a choice leads to compatible elements on the fluid-structure interface, since the meshes are conformal. Fig. 5 provides an example of discretization mesh.

<figure><img src="image_6.png" /><figcaption>Figure 5: Examples of mesh in the fluid-structure problem, around the head of the cilium. a) in the fluid domain the mesh is deformed after several time steps. b) a remeshing technique in the fluid domain is used.</figcaption></figure>

Uzawa's algorithm is summarized in Algorithm 1. At time $t_k$, given a parameter $\rho > 0$ and an initial guess for the Lagrange multiplier $\lambda^{k,0}$, we construct for all $j \ge 0$ a new function $\lambda^{k,j+1}$ by $i$) solving the Stokes problem appearing in problem (32) with $\lambda^k = \lambda^{k,j}$ as Neumann boundary condition on the fluid-structure interface, ii) solving the elasticity problem with $\lambda^k = \lambda^{k,j}$ as Neumann boundary condition on $\Gamma$ and iii) updating the Lagrange multiplier using the fourth equation in (32):

$$ \lambda^{k,j+1} = \lambda^{k,j} + \rho(\delta t \gamma_{\Gamma}(u_f^{k,j} \circ \Phi(d_s^{k-1})) - \gamma_{\Gamma}(d_s^{k,j}) + \gamma_{\Gamma}(d_s^{k-1})), \text{ on } \Gamma. $$

In practice, the function $\lambda^{k,j}$ is approximated by a $P_1$ function in the whole domain $\Omega_s$, but only its values on the boundary $\Gamma$ are used in the resolution of both the Stokes and elasticity problems. This algorithm is known to converge in the sense that, if the parameter $\rho$ is chosen small enough, the sequences $(u_f^{k,j})_{j \in \mathbb{N}}$, $(p_f^{k,j})_{j \in \mathbb{N}}$, $(d_s^{k,j})_{j \in \mathbb{N}}$ and $(\lambda^{k,j})_{j \in \mathbb{N}}$ converge to $u_f^k, p_f^k, d_s^k$ and $\lambda^k$, for all $k \ge 0$.
**Algorithm 1** Resolution of problem (32) at time $t_k$ (Uzawa's algorithm).

Choose a parameter $\rho$ and an initial guess $\lambda^{k,0}$ for the Lagrange multiplier.
$j = 0$.

**while** convergence criteria are not satisfied **do**

Compute the solution of the fluid problem ($u_f^{k,j}, p_f^{k,j}$), with $\lambda^k = \lambda^{k,j}$.

Compute the solution of the structure problem $d_s^{k,j}$ with $\lambda^k = \lambda^{k,j}$.

Update the Lagrange multiplier:

$$ \lambda^{k,j+1} = \lambda^{k,j} + \rho(\delta t \gamma_{\Gamma}(u_f^{k,j} \circ \Phi(d_s^{k-1})) - \gamma_{\Gamma}(d_s^{k,j}) + \gamma_{\Gamma}(d_s^{k-1})), \text{ on } \Gamma. $$

Update the number of iterations $j = j + 1$.

**end while**

Once Uzawa's algorithm has converged, we recover the velocity of the fluid at time $t_k$, $u_f^k$, and the pressure of the fluid at time $t_k$, $p_f^k$, both defined in the deformed domain $\Omega_f^k$. We also obtain the displacement of the structure at time $t_k$, $d_s^k$, expressed in the reference solid configuration $\Omega_s$. The Lagrange multiplier at time $t_k$, $\lambda^k$, is also obtained in $\Omega_s$ and is used at the next time step as an initialization for Uzawa's algorithm. Besides, to go to the next time step, we move the fluid and structure meshes using the displacement of the structure. The solid mesh at time $t_{k+1}$ is directly obtained by moving the mesh representing the domain $\Omega_s$ with the displacement of the structure at time $t_k$. For the fluid domain, we can not use directly the fluid displacement since fluid recirculations may occur and it could lead to poor quality meshes. Then, we construct the deformation $\Phi(d_s^k)$ which maps the reference fluid domain $\Omega_f$ to the deformed fluid domain at time $t_{k+1}$, $\Omega_f^{k+1}$. For that purpose, we solve a problem of linearized elasticity in the domain $\Omega_f^k$ with Dirichlet boundary conditions on the deformed fluid-structure interface $\Gamma^k$ given by the displacement of the structure $d_s^k$:

$$ \left\{ \begin{array}{ll} \text{find } d_f^k : \Omega_f^k \to \mathbb{R}^n \text{ such that} & \\ -\text{div}(\sigma_s(d_f^k)) = 0, & \text{in } \Omega_f^k, \\ d_f^k = d_s^k \circ \Phi^{-1}(d_s^{k-1}), & \text{on } \Gamma^k, \\ d_f^k = 0, & \text{on } \Gamma_f^k, \end{array} \right. \qquad (33) $$

where $\sigma_s(d_f^k)$ is the linearized elasticity stress tensor defined by (19). Thus, the deformation of the fluid domain is given by

$$ \Phi(d_s^k) = \Phi(d_s^{k-1}) + d_f^k, \text{ in } \Omega_f. $$

The displacement of the fluid domain, $d_f^k$, obtained through this process is smoother than the real displacement of the fluid, since the operator of linearized elasticity extends the displacement of the structure from the interface $\Gamma^k$ to the whole domain $\Omega_f^k$ with a diffusion process. However, this regularization process is not enough when the displacement of the structure is large, in which case the fluid domain has to be remeshed. More particularly, the interior of the deformed fluid domain is remeshed but we never modify the boundary, since we want the fluid mesh and the solid mesh to be conformal at the interface $\Gamma$. The algorithm is summarized in Algorithm 2.

The results shown in the remaining of the present section have been obtained using the finite element softwares FEniCS, see [45], and FreeFem++, see [46]. The remeshing of the fluid domain is done with the Mmg software, see [47].

## 4.2 Influence of the fluid viscosity

We start our numerical investigations with the study of the influence of the fluid viscosity on the fluid-structure system. In particular, we examine its effects on the deformations of the structure and on the displacements of the fluid. We consider the case of one cilium beating in a viscous fluid of viscosity $\mu_f$ with the activity scenario be given by (7). The computational domain is a two-dimensional box of dimensions $L_1 = 20\mu m$ and $L_2 = 10\mu m$. Initially, the cilium is represented as a vertical thin cylinder of length $L_c = 6.5\mu m$ and radius $r_c = 0.2\mu m$, anchored at the bottom. For the elasticity parameters of the structure we take $E_s = 10^6 \text{pN} \cdot \text{m}^{-2}$ and $\nu_s = 0.49$.

In Fig. 6, we show the result of the simulation, at different times, of the fluid-structure problem with a fluid viscosity $\mu_f = 0.01 \text{pN} \cdot \text{m}^{-2} \cdot \text{s}$ (ten times more viscous than water). Even though the internal activity of the cilium is imposed, the emerging beating pattern results from both the elasticity properties of the solid and the strong coupling with the surrounding fluid. The fluid velocity is represented as streamlines and glyphs while, in the
**Algorithm 2** Algorithm for the resolution of the fluid-structure problem.

Choose a maximal number of iterations $N$.

Choose an initial guess $\lambda^{-1}$ for the Lagrange multiplier at $t = 0$.

Set $k = 0$.

while $k < N$ do

Initialize Uzawa's algorithm with $\lambda^{k,0} = \lambda^{k-1}$.

Solve fluid-structure problem (32) with Uzawa's algorithm (1) and obtain

$$ (u_f^k, p_f^k, d_s^k, \lambda^k). $$

Find the displacement of the fluid domain $d_f^k$ by solving problem (33).

Move the solid and fluid meshes.

Remesh the fluid domain if necessary.

Update the number of iterations $k = k + 1$.

end while

<figure><img src="image_7.png" /><figcaption>Figure 6: One cilium beating in a viscous fluid of viscosity $\mu_f = 0.01\text{pN} \cdot \text{\textmu m} \cdot \text{s}$ at different times: a) $t = 0.091\text{s}$, b) $t = 0.110\text{s}$, c) $t = 0.120\text{s}$ and d) $t = 0.166\text{s}$. The activity scenario of the structure is given by (7).</figcaption></figure>
deformed solid domain, we plot at each point of the mesh the Frobenius norm of the Green-Lagrange strain tensor, defined in the reference configuration by

$$ \sqrt{E(d_s(x)) : E(d_s(x))}, \quad \forall x \in \Omega_s. $$

The Green-Lagrange strain tensor $E(d_s)$ measures the deformations of the solid material. It can also be seen as the variations of the deformations of the structure compared to rigid deformations. In the present context, the Frobenius norm of the Green-Lagrange tensor gives us a general scalar criterion to observe the action of the fluid viscosity on the deformations of the structure.

In Fig. 7, we compare the solutions of the fluid-structure problem for different values of the fluid viscosity at different times, when the cilium beats with the same activity scenario. As the viscosity of the fluid increases, we remark that the displacements of the cilium are less important because the surrounding fluid acts like a damper on the structure. In particular, this result enlights the particularity of our model, which considers active structures with a finite internal energy and whose displacements result from a strong interaction between the surrounding fluid and the (elastic and internal) properties of the cilium. In comparison, models considered in [14] and [17] for the cilia does not take the effects of the fluid into account, thus the beating of the cilia does not change when the viscosity of the fluid is modified. To quantify this phenomenon, we study the evolution of the mean of the Frobenius

<figure><img src="image_8.png" /><figcaption>Figure 7: One cilium beating in a viscous fluid at time $t = 0.128s$ (left) and at time $t = 0.171s$ (right) for different values of the fluid viscosity: $\mu_f = 0.01pN \cdot \mu m^{-1} \cdot s$ (a, b), $\mu_f = 0.02pN \cdot \mu m^{-1} \cdot s$ (c, d) and $\mu_f = 0.04pN \cdot \mu m^{-1} \cdot s$ (e, f). The activity scenario of the structure is given by (7).</figcaption></figure>

norm of the gradient of the displacement in $\Omega_s$, defined at time $t_k$ by

$$ \frac{1}{|\Omega_s|} \int_{\Omega_s} \sqrt{\nabla d_s^k : \nabla d_s^k} $$
Actually, this quantity is a measure of the deformation of the structure and gives information on the damping effects of the fluid. In Fig. 8, we plot the solid deformations as a function of the time and observe the increase of the damping effect of the fluid on the structure with the increase of the viscosity. Indeed, the more important the

<figure><img src="image_9.png" /><figcaption>Figure 8: Mean deformation of the solid as a function of the time for different values of the fluid viscosity (in pN·µm⁻²·s).</figcaption></figure>

viscosity is, the smaller are the deformations of the structure. To go further, we represent in Fig. 9 a) the time average of the deformation with respect to the viscosity of the fluid, defined for a number of iterations $N > 0$ by

$$ \frac{1}{t_N |\Omega_s|} \sum_{k=1}^{N} (t_k - t_{k-1}) \int_{\Omega_s} \sqrt{\nabla d_s^k : \nabla d_s^k} $$

In particular, the deformations of the cilium decrease by half when the viscosity increases from 0.01 to 0.05pN·µm⁻²·s. Moreover, it suggests that if the fluid is too viscous, the cilium will not be able to deform at all. In Fig. 8, we also remark that the increase in the fluid viscosity induces a time shift in the beating of the cilium. This time shift is represented in Fig. 9 b), with respect to the viscosity of the fluid.

Another important criterion for the study of this fluid-structure system, is the capacity of the cilium to propel the surrounding fluid. To study this phenomenon we plot in Fig. 10 the evolution of the average horizontal fluid velocity in $\Omega_f^k$ with respect to the time, defined at time $t_k$ by

$$ \int_{\Omega_f^k} u_f^k \cdot e_1, $$

where $e_1$ is the first vector of the canonical basis of $\mathbb{R}^2$. The increase in the fluid viscosity induces a diminution of the instantaneous velocity of the fluid and we also observe the same time shift as in Fig. 9 b). Moreover, the order of $100\mu m \cdot s^{-1}$ that we find for the fluid velocity, is consistent with the typical velocity of micro-swimmer at low-Reynolds number. To quantify the transport of fluid by the cilium, we represent in Fig. 11 b) the mean horizontal fluid velocity in time with respect to the viscosity, defined for a number of iterations $N > 0$ by

$$ \frac{1}{t_N} \sum_{k=1}^{N} (t_k - t_{k-1}) \int_{\Omega_f^k} u_f^k \cdot e_1. \qquad (34) $$

With the chosen activity scenario for the cilium, the mean horizontal fluid velocity is negative, which means that the fluid goes to the left on average. Furthermore, the velocity tends to zero as the viscosity increases, which is
<figure><img src="image_10.png" /><figcaption>Figure 9: Time average of the solid deformation as a function of the viscosity of the fluid.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Mean horizontal velocity of the fluid in function of the time for different values of the fluid viscosity (in pN · µm⁻² · s).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: a) Time average of the horizontal fluid velocity in function of the viscosity of the fluid and b) horizontal velocity of the fluid as a function of the height for different values of the viscosity.</figcaption></figure>

coherent with what is found for the deformation of the structure: when the viscosity increases the damping effects of the fluid tend to immobilize the whole system.

In particular, this behavior illustrates well the decrease of the mucociliary clearance efficiency in the lungs when the viscosity of the mucus increases, which is for example the case in diseases such as cystic fibrosis. In Fig. 11 b), we plot the mean horizontal velocity of the fluid as a function of the height for different values of the viscosity. Recalling that the length of the cilium is 6.5µm, it appears that the velocity of the fluid grows linearly with the height close to the cilium and is almost constant far from the bottom of the domain. When the viscosity of the fluid increases, the damping effects make the fluid go slower. However, it is interesting to remark that the velocity of the fluid at the top of the domain is more important for a viscosity of $\mu_f = 0.02\text{pN} \cdot \mu\text{m}^{-2} \cdot \text{s}$ than for $\mu_f = 0.01\text{pN} \cdot \mu\text{m}^{-2} \cdot \text{s}$, which may suggest the existence of an optimal value of the viscosity for the transport of the fluid far from the bottom. These velocity profiles are very similar to those found in [14] and in [48], only considering the bottom layer of their bi-fluid models for the mucociliary clearance process.

## 4.3 Interaction of two cilia in a viscous fluid

A cilium rarely beats alone in a fluid but is often in interaction with other cilia. Thus, the study of the interaction of two cilia beating in a viscous fluid is of particular interest if we aim to study a more complex system. The numerical simulations that we present in this subsection correspond to two cilia of length $L_c = 6.5\mu\text{m}$ and radius $r_c = 0.2\mu\text{m}$ anchored at positions $(x_{c1}, 0)$ and $(x_{c2}, 0)$. The whole domain of simulation is the same as the one described in the previous subsection, i.e. a two-dimensional box of dimensions $L_1 = 20\mu\text{m}$ and $L_2 = 10\mu\text{m}$. The activity scenario for both cilia is the one defined by (7), but we add a phase shift between cilia whose effects on the fluid-structure system is studied. In addition to the phase shift, the effects of the distance between cilia on the fluid-structure system is also investigated. In Fig. 12, we represent the two cilia beating in the fluid at two different times. In this simulation, the cilia are anchored at positions $x_{c1} = 7.5\mu\text{m}$ and $x_{c2} = 12.5\mu\text{m}$. The elasticity parameters of the structure are $E_s = 10^6\text{pN} \cdot \mu\text{m}^{-2}$ and $\nu_s = 0.49$ and the viscosity of the fluid is $\mu_f = 0.01\text{pN} \cdot \mu\text{m}^{-2} \cdot \text{s}$. The beating frequency of both cilia is $f_a = 10\text{Hz}$ and the intensity of the activity is $C_a = 3\text{pN} \cdot \mu\text{m}$. Moreover, we consider a phase shift $\delta\phi = 0.075\text{s}$ between the activity scenarios of the cilia, i.e. the second cilia (at position $x_{c2}$) starts to beat at $t = 0.075\text{s}$. The emerging beating pattern for both cilia is more complicated than what we observed in the previous subsection, since each cilium is now subjected to the fluid flow generated not only by its own activity but also by the one of the other cilium.

Actually, this behavior can be seen in Fig. 13, where we plot the trajectories of the cilia summits during one beating period, for two different values of the phase shift: $\delta\phi = 0.025\text{s}$ in Fig. 13 a) and $\delta\phi = 0.075\text{s}$ in Fig. 13 b). In each simulation the two cilia have completely different beating patterns, which are also different from one simulation to the other. Thus, both cilia act on each other by means of hydrodynamic interactions, which depend on several parameters: the viscosity of the fluid, the elastic parameters of the structures, the parameters of the activity scenario
<figure><img src="image_13.png" /><figcaption>Figure 12: Two cilia beating in a viscous fluid with $\mu_f = 0.01\text{pN} \cdot \text{\textmu m}^{-2} \cdot \text{s}$ at a) $t = 0.788\text{s}$ and b) $t = 0.83\text{s}$. The cilia are anchored at $x_{c1} = 7.5\text{\textmu m}$ and $x_{c2} = 12.5\text{\textmu m}$. The phase shift between the two cilia is $\delta\phi = 0.075\text{s}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Trajectories of the two cilia summits for different values of the phase shift between their activity scenario: a) $\delta\phi = 0.025$ and b) $\delta\phi = 0.075$.</figcaption></figure>
and the geometrical parameters.

In this subsection, we assume that the two cilia are identical, except for the presence of a phase shift between their activity scenario, and that the viscosity of the fluid is fixed at $\mu_f = 0.01\text{pN} \cdot \mu\text{m}^{-2} \cdot \text{s}$. Thus, it only remains two parameters of interest to study, namely the phase shift between the activity scenarios of the cilia, denoted by $\delta\phi$, and the distance between their anchorage points, denoted by $\delta x_c$.

We start with the study of the influence of the phase shift. For that matter, we solve the fluid-structure problem during six beating periods and consider different values for the phase shift from $\delta\phi = 0\text{s}$ to $\delta\phi = 0.1\text{s}$. In these simulations, all other parameters are fixed. In particular, the distance between the two cilia is $\delta x_c = 5\mu\text{m}$ and the activity scenario is given by (7). In Fig. 14 a), we plot the mean horizontal fluid velocity as a function of the height for different values of the phase shift. Depending on the phase shift, it appears that the behavior of the system is totally different. Actually, for $\delta\phi = 0.025\text{s}$ the mean velocity of the fluid is negative, such that the fluid goes to the left (this is what we observed with one cilium), whereas for $\delta\phi = 0.075\text{s}$ the fluid goes to the right on average. This phenomenon is confirmed in Fig. 14 b), which represents the average fluid velocity with respect to the phase shift. For a phase shift from 0 to approximately 0.05, the fluid moves to the left on average, since its velocity is negative, and from 0.05 to 0.1, the fluid moves to the right. Thus, for the internal activity scenario considered in these simulations, the transport of the fluid by the actives structures strongly depends on the phase shift between their activity.

<figure><img src="image_15.png" /><figcaption>Figure 14: a) Horizontal velocity of the fluid as a function of the height for different values of the phase shift and b) time average of the horizontal fluid velocity in function of the phase shift.</figcaption></figure>

Similarly, we study the influence of the distance between the two cilia on the fluid-structure system. To that aim several simulations are realised by considering different values for the distance between them, from $\delta x_c = 2\mu\text{m}$ to $\delta x_c = 6.5\mu\text{m}$. The other parameters remain constant and, in particular, we consider a phase shift $\delta\phi = 0.075\text{s}$ between the activity scenarios of the two cilia. The velocity profiles presented in Fig. 15 a) are similar, which suggest that the distance between the two cilia is less determinant than the phase shift. Here, the profile with $\delta x_c = 5\mu\text{m}$ correspond to the profile with $\delta\phi = 0.075\text{s}$ in Fig. 14 a). In Fig. 15 b), the average of horizontal fluid velocity is represented with respect to the distance between the two cilia. For all values of the distance, the fluid moves to the right on average. However, the fluid velocity does not seem to be linear with respect to this distance. In particular, a maximum of the transport efficiency is reached at approximately $\delta x_c = 3\mu\text{m}$.

## 4.4 Discussion

The continuum model for active structures that we have presented enables to fully take into account the fluid-structure interactions between the cilia and the surrounding fluid, including the hydrodynamic forces due to the deformations of neighboring cilia. For a given activity scenario of the cilia, we were able to study the influence on the fluid-structure system of the viscosity of the fluid, when one cilium is considered, as well as the influence of the phase shift and the distance between cilia, when two cilia are involved. As a conclusion we showed that, even with an imposed internal activity scenario for the active structures, hydrodynamic forces really are the principal
<figure><img src="image_16.png" /><figcaption>Figure 15: a) Horizontal velocity of the fluid as a function of the height for different values of the distance between the two cilia and b) time average of the horizontal fluid velocity in function of the distance.</figcaption></figure>

ingredient that controls the beating patterns of cilia, since completely different deformations of the structures are observed depending of the different physical and geometrical parameters.

Even though two-dimensional simulations are a good start for the study of such a complex system, it is not quite representative of the real fluid-structure interactions for microorganisms. Actually, in two space dimensions the system is much more constrained than in three space dimensions, since the fluid can not circumvent a cilium. It results that hydrodynamic interactions are less important in three space dimensions and their effects on the emergence of the beating patterns of cilia should be further investigated.

The model as well as the numerical method that we proposed are both compatible with thee-dimensional simulations. As an example, we plot in Fig. 16 the result of the simulation of an active elastic elongated structure beating in a fluid modeled by the Stokes equations in three space dimensions. The activity scenario is the same as the one used for two-dimensional simulations, but we can remark that most of the fluid goes around the cilium instead of above, which leads to slower fluid velocities. In this three space dimensions simulation, the computational cost of the method is much more important and necessitates the development of more sophisticated numerical tools. With Uzawa's algorithm, one should consider the use of preconditioners and parallel finite element solvers to really improve the performance of the method. A different approach would be to consider some methods on non conformal and fixed meshes, which then could be implemented with the use of fast solvers (e.g. solvers using the discrete fast Fourier transform).

## References

[1] G. R. Fulford, J. R. Blake, Muco-ciliary transport in the lung, Journal of Theoretical Biology 121 (4) (1986) 381-402.

[2] C. J. Brokaw, Simulating the effects of fluid viscosity on the behaviour of sperm flagella, Mathematical Methods in the Applied Sciences 24 (17-18) (2001) 1351-1365.

[3] S. M. King, G. J. Pazour, Cilia: motors and regulation, Vol. 92, Academic Press, 2009.

[4] G. Taylor, Analysis of the swimming of microscopic organisms, in: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, Vol. 209, The Royal Society, 1951, pp. 447-461.

[5] G. Taylor, The action of waving cylindrical tails in propelling microscopic organisms, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 211 (1105) (1952) 225-239.

[6] J. R. Blake, Infinite models for ciliary propulsion, Journal of Fluid Mechanics 49 (2) (1971) 209-222.
<figure><img src="image_17.png" /><figcaption>Figure 16: Three-dimensional simulation of a cilium-like structure in a viscous fluid.</figcaption></figure>

[7] M. J. Lighthill, On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers, Communications on Pure and Applied Mathematics 5 (2) (1952) 109–118.

[8] J. R. Blake, A spherical envelope approach to ciliary propulsion, Journal of Fluid Mechanics 46 (1) (1971) 199–208.

[9] J. Gray, G. J. Hancock, The propulsion of sea-urchin spermatozoa, Journal of Experimental Biology 32 (4) (1955) 802–814.

[10] G. J. Hancock, The self-propulsion of microscopic organisms through liquids, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 217 (1128) (1953) 96–121.

[11] M. J. Lighthill, Flagellar hydrodynamics, SIAM Review 18 (2) (1976) 161–230.

[12] R. D. Dresdner, D. F. Katz, S. A. Berger, The propulsion by large amplitude waves of uniflagellar micro-organisms of finite length, Journal of Fluid Mechanics 97 (3) (1980) 591–621.

[13] Y. Ding, J. C. Nawroth, M. J. McFall-Ngai, E. Kanso, Mixing and transport by ciliary carpets: a numerical study, Journal of Fluid Mechanics 743 (2014) 124–140.

[14] L. Lacouture, Modélisation et simulation du mouvement de structures fines dans un fluide visqueux: application au transport mucociliaire, Ph.D. thesis, Université Paris-Saclay, https://tel.archives-ouvertes.fr/tel-01366885/ (2016).

[15] A. Dauptain, J. Favier, A. Bottaro, Hydrodynamics of ciliary propulsion, Journal of Fluids and Structures 24 (8) (2008) 1156 – 1165.

[16] R. Chatelin, Méthodes numériques pour l’écoulement de stokes 3D: fluides à viscosité variable en géométrie complexe mobile; application aux fluides biologiques, Ph.D. thesis, Université de Toulouse, Université Toulouse III-Paul Sabatier, http://thesesups.ups-tlse.fr/2398/ (2013).

[17] S. Chateau, J. Favier, U. D'Ortona, S. Poncet, Transport efficiency of metachronal waves in 3D cilium arrays immersed in a two-phase flow, Journal of Fluid Mechanics 824 (2017) 931–961.
[18] K. E. Machin, Wave propagation along flagella, Journal of Experimental Biology 35 (4) (1958) 796-806.

[19] C. J. Brokaw, Bend propagation along flagella, Nature 209 (5019) (1966) 161-163.

[20] J. Lubliner, J. J. Blum, Model for bend propagation in flagella, Journal of Theoretical Biology 31 (1) (1971) 1-24.

[21] C. J. Brokaw, Computer simulation of flagellar movement: I. demonstration of stable bend propagation and bend initiation by the sliding filament model, Biophysical Journal 12 (5) (1972) 564-586.

[22] M. Hines, J. J. Blum, Bend propagation in flagella. I. Derivation of equations of motion and their simulation, Biophysical Journal 23 (1) (1978) 41-57.

[23] C. J. Brokaw, Computer simulation of flagellar movement. VI. Simple curvature-controlled models are incompletely specified, Biophysical Journal 48 (4) (1985) 633-642.

[24] K. E. Machin, The control and synchronization of flagellar movement, Proceedings of the Royal Society of London. Series B. Biological Sciences 158 (970) (1963) 88-104.

[25] C. J. Brokaw, Computer simulation of flagellar movement IX. Oscillation and symmetry breaking in a model for short flagella and nodal cilia, Cell Motility and the Cytoskeleton 60 (1) (2005) 35-47.

[26] M. Murase, The dynamics of cellular motility, John Wiley and Sons, 1992.

[27] R. H. Dillon, L. J. Fauci, An integrative model of internal axoneme mechanics and external fluid dynamics in ciliary beating, Journal of Theoretical Biology 207 (3) (2000) 415-430.

[28] R. H. Dillon, L. J. Fauci, C. Omoto, Mathematical modeling of axoneme mechanics and fluid dynamics in ciliary and sperm motility, Dynamics of Continuous Discrete and Impulsive Systems Series A 10 (2003) 745-758.

[29] X. Yang, R. H. Dillon, L. J. Fauci, An integrative computational model of multiciliary beating, Bulletin of Mathematical Biology 70 (4) (2008) 1192-1215.

[30] S. Lukens, X. Yang, L. J. Fauci, Using Lagrangian coherent structures to analyze fluid mixing by cilia, Chaos: An Interdisciplinary Journal of Nonlinear Science 20 (1) (2010) 017511.

[31] S. Gueron, K. Levit-Gurevich, A three-dimensional model for ciliary motion based on the internal 9+2 structure, Proceedings of the Royal Society of London. Series B: Biological Sciences 268 (1467) (2001) 599-607.

[32] S. M. Mitran, Metachronal wave formation in a model of pulmonary cilia, Computers and Structures 85 (11) (2007) 763-774.

[33] D. Ambrosi, S. Pezzuto, Active stress vs. active strain in mechanobiology: constitutive issues, Journal of Elasticity 107 (2) (2012) 199-212.

[34] N. P. Smith, A computational study of the interaction between coronary blood flow and myocardial mechanics, Physiological Measurement 25 (4) (2004) 863.

[35] A. V. Panfilov, R. H. Keldermann, M. P. Nash, Self-organized pacemakers in a coupled reaction-diffusion-mechanics system, Physical Review Letters 95 (25) (2005) 258104.

[36] D. K. Bogen, S. A. Rabinowitz, A. Needleman, T. A. Mcmahon, W. H. Abelmann, An analysis of the mechanical disadvantage of myocardial infarction in the canine left ventricle, Circulation Research 47 (5) (1980) 728-741.

[37] S. Pezzuto, D. Ambrosi, A. Quarteroni, An orthotropic active-strain model for the myocardium mechanics and its numerical approximation, European Journal of Mechanics-A/Solids 48 (2014) 83-96.

[38] Y. Payan, J. Ohayon, Biomechanics of living organs: hyperelastic constitutive laws for finite element modeling, Academic Press, 2017.

[39] M. E. J. Holwill, P. Satir, Generation of propulsive forces by cilia and flagella, Cytomechanics (1987) 120-130.

[40] E. M. Purcell, Life at low Reynolds number, American Journal of Physics 45 (1) (1977) 3-11.
[41] F. Vergnet, Structures actives dans un fluide visqueux : modélisation, analyse mathématique et simulations numériques, Ph.D. thesis, Université Paris-Saclay, https://tel.archives-ouvertes.fr/tel-02194265/ (2019).

[42] M. E. Bogovski, Solution of the first boundary value problem for an equation of continuity of an incompressible medium, Doklady Akademii Nauk SSSR 248 (5) (1979) 10371040.

[43] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers, Revue Française d'Automatique, Informatique, Recherche Opérationnelle. Analyse Numérique 8 (R2) (1974) 129-151.

[44] F. Boyer, P. Fabrie, Eléments d'analyse pour l'étude de quelques modèles d'écoulements de fluides visqueux incompressibles, Vol. 52, Springer Science and Business Media, 2005.

[45] A. Logg, K.-A. Mardal, G. Wells, Automated solution of differential equations by the finite element method: The FEniCS book, Vol. 84, Springer Science and Business Media, 2012.

[46] F. Hecht, New development in FreeFem++, Journal of Numerical Mathematics 20 (3-4) (2012) 251-266.

[47] C. Dapogny, C. Dobrzynski, P. Frey, Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems, Journal of Computational Physics 262 (2014) 358-378.

[48] D. J. Smith, E. A. Gaffney, J. R. Blake, Modelling mucociliary clearance, Respiratory Physiology and Neurobiology 163 (1) (2008) 178-188.