Lagrangian mechanics of active systems

Abstract

We present a multi-scale modeling and simulation framework for low-Reynolds number hydrodynamics of shape-changing immersed objects, e.g., biological microswimmers and active surfaces. The key idea is to consider principal shape changes as generalized coordinates and define conjugate generalized hydrodynamic friction forces. Conveniently, the corresponding generalized friction coefficients can be pre-computed and subsequently reused to solve dynamic equations of motion fast. This framework extends Lagrangian mechanics of dissipative systems to active surfaces and active microswimmers, whose shape dynamics is driven by internal forces. As an application case, we predict in-phase and anti-phase synchronization in pairs of cilia for an experimentally measured cilia beat pattern.

Graphic Abstract

Biological hydrodynamics Biology provides ample examples of active shape changes in fluid environments: Bacteria like E. coli rotate helical prokaryotic flagella to swim [1], other bacteria like Spiroplasma propagate twist waves along their flexible body [2], sperm cells and motile algae posses slender cell appendages termed cilia (or eukaryotic flagella), whose regular bending waves propel these cells in a fluid [3, 4]. On epithelial surfaces, collections of beating cilia transport biological fluids such as mucus in airways, cerebrospinal fluid in brain ventricles, and oviduct fluid in the Fallopian tubes [5, 6]. In addition to their important role in self-propulsion and fluid transport, these model systems enable us to learn about internal force generation mechanisms in these cells, such as the collective dynamics of molecular motors inside cilia [7,8,9,10]. On larger scales, the interaction of many shape-changing units leads to the spontaneous formation of spatiotemporal patterns, e.g., in dense suspensions of microswimmers [11], or collections of cilia exhibiting metachronal coordination [12].

These examples represent a class of fluid–structure interaction problems, where shape-changing active structures exert forces on the surrounding fluid, while the surrounding passive fluid exerts hydrodynamic friction forces back on these active structures. These hydrodynamic forces may affect the active shape dynamics; examples include the torque–velocity relationship of rotating prokaryotic flagella [13], the load response of beating cilia and eukaryotic flagella [10, 14], as well as minimal model swimmers [15,16,17]. Closed feedback loops between passive fluids and active structures can lead to emergent dynamics; examples include spontaneous pattern formation in dense microswimmer suspensions [11, 18] or (hydrodynamic) synchronization of beating cilia and flagella [12, 19,20,21,22,23].

Common hydrodynamics methods at low Reynolds numbers At the relevant length and time scales, viscous drag dominates inertia, corresponding to low Reynolds numbers [24,25,26]. In the limit of zero Reynolds numbers, the Navier–Stokes equation of hydrodynamics simplifies to the Stokes equation. Although the Stokes equation is linear, hydrodynamic computations can still be costly, because hydrodynamic interactions are long-ranged [27].

In the past, different computational methods of different degrees of approximation have been used in the community, including resistive force theory for slender filaments, which includes short-range, but not long-range hydrodynamic interactions [28,29,30], the more refined method of slender-body theory, which considers a line distribution of hydrodynamic singularities (point forces) along a filament [31,32,33], or multi-particle collision dynamics, which replaces the continuum description of the Stokes equation by the stochastic dynamics of a large number of “fluid particles” [34,35,36]. Despite its applicability for large-scale problems [37], the stochastic nature of the MPCD algorithm introduces algorithm-specific fluctuations, which can be impractical if one wants to study the role of biological noise. Lattice Boltzmann methods similarly rely on fictitious “fluid particles”, for which in each time step both a streaming and a collision step is performed [38]. Finally, boundary element methods convert the problem of solving the Stokes equation in three-dimensional space to a two-dimensional boundary integral problem of finding a surface distribution of forces on a moving boundary surface. Boundary element methods are similar in spirit to slender-body methods, but less susceptible to issues of regularization, since a two-dimensional distribution of forces is used. Modern algorithms use fast multi-pole methods that solve a tree of hierarchically coarse-grained subproblems instead of solving a single large linear system when computing the force distribution on a surface [39,40,41].

Irrespective of the hydrodynamic computation method used, it can be computationally costly to calculate a solution of the Stokes equation in every time step, while simulating the dynamics of a shape-changing microswimmer or an active surface.

Lagrangian mechanics In this article, we present a multi-scale simulation framework, where the Stokes equation has to be solved only in an initial step for a small set of principal shape modes of a shape-changing surface. The resultant surface distributions of hydrodynamic friction forces define generalized hydrodynamic friction coefficients by a projection method of Lagrangian mechanics [10, 42,43,44,45,46,47,48]. These scalar friction coefficients are independent of the velocity of the moving surface. Once tabulated, these friction coefficients provide a look-up table for subsequent fast simulations of shape dynamics and active motion. Specifically, we view principal shape changes of an active surface as generalized coordinates, for which we compute conjugate generalized friction forces. We obtain effective equations of motion for the generalized coordinates from a force balance between these generalized friction forces and active driving forces. These active driving forces coarse-grain the internal active processes that drive the active shape changes of the surface (such as the collective dynamics of molecular motors). Importantly, these a priori unknown active driving forces can be calibrated for a reference case (e.g., using experimental data) and then used to extrapolate to other application cases of interest. Therefore, our framework extends Lagrangian mechanics of dissipative systems to active surfaces and active microswimmers, whose shape dynamics is driven by active forces.

1 Notation: Stokes equation and hydrodynamic dissipation

Fluid dynamics at the scale of individual biological cells is characterized by low Reynolds numbers, i.e., viscous effects commonly dominate over inertia [24,25,26]. Correspondingly, fluid flow is described by the Stokes equation, which reads for an incompressible Newtonian fluid in the absence of body forces in the bulk [27]

$$\begin{aligned} \varvec{0}= - \varvec{\nabla }p + \mu \,\varvec{\nabla }^2\mathbf {u}, \end{aligned}$$

(1)

with incompressibility condition $\varvec{\nabla }\cdot \mathbf {u}= \varvec{0}$. Here, $\mathbf {u}(\mathbf {x})$ denotes the flow velocity, $p(\mathbf {x})$ the pressure field, and $\mu$ the dynamic viscosity of the fluid.

The total stress tensor $\varvec{\sigma }$ for an incompressible fluid depends on both the hydrostatic pressure p and the symmetrized strain rate tensor $\varvec{\Delta }$ [27]

$$\begin{aligned} \varvec{\sigma }= -p\,\varvec{1}+ 2\mu \,\varvec{\Delta },\quad \varvec{\Delta }= \frac{1}{2} \left[ \varvec{\nabla }\otimes \mathbf {u}+ (\varvec{\nabla }\otimes \mathbf {u})^T \right] . \end{aligned}$$

(2)

Thus, the Stokes equation, Eq. (1) could be equivalently written as $\varvec{0}= \varvec{\nabla }\cdot \varvec{\sigma }$ in the bulk of the fluid. Special conditions apply at boundaries.

No-slip boundary condition for an active surface We consider a surface $\mathcal {S}$ immersed in the fluid that changes its shape as a function of time. For example, $\mathcal {S}$ may represent the outer surface of a shape-changing microswimmer, or even the combined surface for a collection of microswimmers. We introduce the surface velocity $\mathbf {v}(\mathbf {x},t)$ for each point $\mathbf {x}\in \mathcal {S}$ at time t.

We impose a no-slip boundary condition at this surface, i.e., require that the local velocity $\mathbf {u}(\mathbf {x})$ of fluid flow matches the local velocity $\mathbf {v}(\mathbf {x})$ of the surface for each surface point

$$\begin{aligned} \mathbf {u}(\mathbf {x}) = \mathbf {v}(\mathbf {x}) \text { for all } \mathbf {x}\in \mathcal {S}. \end{aligned}$$

(3)

Hydrodynamic friction forces A shape change of the surface $\mathcal {S}$ induces a flow field $\mathbf {u}(\mathbf {x})$ with corresponding stress tensor field $\varvec{\sigma }(\mathbf {x})$. The stress $\sigma (\mathbf {x})$ determines the surface density of forces $\mathbf {f}(\mathbf {x})$ exerted by the surface on the fluid (with units of a stress $\mathrm {N/m^2}$, also called contact force, or traction force density)

$$\begin{aligned} \mathbf {f}(\mathbf {x}) = -\varvec{\sigma }\cdot \mathbf {n}\text { for all } \mathbf {x}\in \mathcal {S}, \end{aligned}$$

(4)

where $\mathbf {n}$ is the surface normal pointing into the fluid. Correspondingly, $-f(\mathbf {x})$ is the surface density of hydrodynamic friction forces exerted by the fluid on the surface. The total force exerted by the surface on the fluid is simply the surface integral of $\mathbf {f}(\mathbf {x})$

$$\begin{aligned} \mathbf {F}= \int _{\mathcal {S}} \! \mathrm{d}^2\mathbf {x}\, \mathbf {f}(\mathbf {x}). \end{aligned}$$

(5)

Analogously, the total torque (with respect to a reference point $\mathbf {x}_0$) exerted by the surface on the fluid is given by

$$\begin{aligned} \mathbf {M}= \int _{\mathcal {S}} \! \mathrm{d}^2\mathbf {x}\, (\mathbf {x}- \mathbf {x}_0) \times \mathbf {f}(\mathbf {x}) . \end{aligned}$$

(6)

Superposition principle The linearity of the Stokes equation of low-Reynolds number flow, Eq. (1), implies a superposition principle for hydrodynamic friction forces, which will be pivotal for the modeling ansatz presented here. Specifically, we consider a boundary condition with rate of displacement $\mathbf {v}$ that is given as a linear combination of velocity distributions $\mathbf {v}_1$ and $\mathbf {v}_2$ as

$$\begin{aligned} \mathbf {v}=\alpha _1\, \mathbf {v}_1 + \alpha _2\, \mathbf {v}_2 , \end{aligned}$$

(7)

with real coefficients $\alpha _1, \alpha _2\in \mathbb {R}$. Then, the resultant flow field $\mathbf {u}$ is given by $\mathbf {u}=\alpha _1\mathbf {u}_1 + \alpha _2\mathbf {u}_2$, while the surface density of hydrodynamic friction forces $\mathbf {f}$ is $\mathbf {f}=\alpha _1\mathbf {f}_1+\alpha _2\mathbf {f}_2$, where $\mathbf {u}_i$ and $\mathbf {f}_i$ denote the flow field and the surface density of hydrodynamic friction forces corresponding to boundary condition $\mathbf {v}_i$, respectively, for $i=1,2$.

Hydrodynamic dissipation We introduce the rate of work $\mathcal {R}^{(h)}$ exerted by the surface on the fluid

$$\begin{aligned} \mathcal {R}^{(h)} = \int _{\mathcal {S}} \! \mathrm{d}^2\mathbf {x}\, \mathbf {v}(\mathbf {x})\cdot \mathbf {f}(\mathbf {x}) . \end{aligned}$$

(8)

For incompressible Newtonian fluids at zero Reynolds number, $\mathcal {R}^{(h)}$ equals the instantaneous rate of hydrodynamic energy dissipation within the fluid [27]. Indeed, let us consider the local dissipation rate, which is given by $\Phi = 2\mu \, {\varvec{\Delta }}:{\varvec{\Delta }}$, where $\varvec{\Delta }:\varvec{\Delta }=\sum _{i,j} \Delta _{ij}\Delta _{ij}$ denotes tensor contraction. The dissipation rate can be rewritten as $\Phi = \varvec{\nabla }\cdot (\mathbf {u}\cdot \varvec{\sigma })$ using Eqs. (1), (2) and the incompressibility condition $\varvec{\nabla }\cdot \mathbf {u}=\varvec{0}$. Gauss divergence theorem now gives [27] (using $\mathbf {u}(\mathbf {x})=\mathbf {v}(\mathbf {x})$ for $\mathbf {x}\in \mathcal {S}$)

$$\begin{aligned} \underbrace{ \int _{\mathcal {S}} \! \mathrm{d}^2\mathbf {x}\, \mathbf {v}(\mathbf {x})\cdot \mathbf {f}(\mathbf {x}) }_\text {power exerted by surface} = \underbrace{ \int _V \! \mathrm{d}^3\mathbf {x}\, \Phi (\mathbf {x}) }_\text {hydrodynamic dissipation in bulk} . \end{aligned}$$

(9)

Here, V denotes the three-dimensional fluid domain with boundary surface $\mathcal {S}$. At finite Reynolds numbers, $\mathcal {R}^{(h)}$ still equals the rate of work exerted by the surface on the fluid, yet this injected energy would be dissipated as heat with a delay, such that Eq. (9) would only hold for time averages.

2 Lagrangian mechanics: generalized coordinates

We consider a shape-changing surface $\mathcal {S}(t)$. While a description of all possible shape changes of $\mathcal {S}$ would require an infinite number of degrees of freedom, in important application cases, we can restrict ourselves to a constrained set of shape changes characterized by a small number of shape coefficients, or generalized coordinates, $q_1,\ldots ,q_n$.

Examples for minimal model swimmers include undulating sheets with a finite set of admissible wavelengths [49], bead distances as in Najafi’s three-sphere swimmer [50], or lever arm angles in Purcell’s the three-link swimmer [51] and Dreyfus’ rotator [52], see Fig. 1. An example for a biological microswimmer would be the rotation angle $\varphi$ of an idealized rigid helical prokaryotic flagellum. Similarly, the regular traveling bending waves of cilia and eukaryotic flagella can be described by an oscillator phase $\varphi$ that characterizes the current position in a periodic shape cycle [45, 53,54,55]. Elastic degrees of freedom arising from waveform compliance can be incorporated in such a framework as additional amplitude degrees of freedom [10, 48].

We introduce the state vector, $\mathbf {q}=(q_1,\ldots ,q_n)$. The shape dynamics of the active surface $\mathcal {S}(t) = \mathcal {S}[\mathbf {q}(t)]$ is thus entirely described by the dynamics of $\mathbf {q}(t)$. In particular, the local rate of surface displacement depends linearly on the generalized velocities $\dot{q}_i$ as

$$\begin{aligned} \mathbf {v}(\mathbf {x}) = \mathbf {w}_1(\mathbf {x};\mathbf {q})\, \dot{q}_1 + \mathbf {w}_2(\mathbf {x};\mathbf {q})\, \dot{q}_2 + \cdots + \mathbf {w}_n(\mathbf {x};\mathbf {q})\, \dot{q}_n , \end{aligned}$$

(10)

where the normalized velocity fields $\mathbf {w}_i(\mathbf {x}) = \partial \mathbf {x}/ \partial q_i$ depend on $\mathbf {q}(t)$ but not on $\dot{\mathbf {q}}(t)$. In fact, Eq. (10) simply generalizes Eq. (7) to the case of generalized coefficients $\alpha _i = \dot{q}_i$ with units of a generalized velocity. Correspondingly, the surface distribution of hydrodynamic friction forces $\mathbf {f}(\mathbf {x})$ is given as a linear combination

$$\begin{aligned} \mathbf {f}(\mathbf {x}) = \mathbf {g}_1(\mathbf {x};\mathbf {q})\, \dot{q}_1 + \cdots + \mathbf {g}_n(\mathbf {x};\mathbf {q})\,\dot{q}_n , \end{aligned}$$

(11)

where the normalized force densities $\mathbf {g}_i(\mathbf {x};\mathbf {q}) = \mathbf {f}_i(\mathbf {x})/\alpha _i$ correspond to the surface density of hydrodynamic friction forces $\mathbf {f}_i(\mathbf {x})$ induced by the velocity field $\mathbf {v}_i(\mathbf {x}) = \alpha _i\, \mathbf {w}_i(\mathbf {x};\mathbf {q})$. An example of a surface velocity field with corresponding surface density of hydrodynamic friction forces is shown in Fig. 2.

The formalism allows to include also rigid body transformation such as translations and rotations of the surface $\mathcal {S}$ in the set of generalized coordinates. Therefore, the self-propulsion of shape-changing microswimmers can be described using the same formalism, see the section of rigid body transformations below.

Fig. 4

Generalized coordinates: Examples. a Undulating sheet with two wave modes. The amplitudes $q_1$, $q_2$ of the wave modes represent generalized coordinates of the shape-changing surface $\mathcal {S}$. b Rigid body motion of a microswimmer in three-dimensional space is characterized by three translational and three rotational degrees of freedom, corresponding to six generalized coordinates: $q_i$ for translations parallel to the $\mathbf {e}_i$-axis, and $q_{i+3}$ for rotations around the $\mathbf {e}_i$-axis, $i=1,2,3$, respectively. c Najafi’s three-sphere swimmer consists of three collinear spherical beads with time-varying bead distances [50], corresponding to two internal degrees of freedom, $q_{6+1}$ and $q_{6+2}$, in addition to the generalized coordinates of rigid body motion. d Purcell’s three-link swimmer consists of three connected segments [51], whose relative angles $q_{6+1}$ and $q_{6+2}$ can be treated as two generalized coordinates. e Similarly, Dreyfus’ rotator consists of three segments connected at a single joint; the relative angles $q_{6+1}$ and $q_{6+2}$ again define generalized coordinates. This shape-changing microswimmer exhibits pronounced rotation in the plane in addition to translational motion, hence its name. f Simplified geometry of the bacterium E. coli with a single prokaryotic flagellum. A rotary motor inside the cell wall can spin the helical flagellum around its central axis; this internal rotational degree of freedom defines a single generalized coordinate $q_{6+1}$ with periodicity of $2\pi$. g Prototypical flagellar beat pattern of a sperm cell, parameterized by a $2\pi$-periodic phase variable, which defines a generalized coordinate $q_{6+1}$. For the amplitude of regular flagellar bending waves and mean flagellar curvature, we used parameters from [30]

Fig. 5

Hydrodynamic friction forces: Example of cilium during power stroke. a Surface velocity $\mathbf {v}_1(\mathbf {x})$ on a shape-changing surface $\mathcal {S}$, here given by a slender cilium (blue) attached to a no-slip boundary surface (gray); the cilium progresses with phase velocity $q_1=\dot{\varphi }_1$ along its periodic beat cycle. For visualization, the three-dimensional shape of the cilium as well as $\mathbf {v}_1(x)$ was projected on the yz-plane (see Fig. 3a for a three-dimensional representation). b Corresponding surface distribution of hydrodynamic friction forces $\mathbf {f}_1(\mathbf {x})$ exerted by the active, shape-changing surface on the surrounding viscous fluid. The force distribution is obtained by solving the Stokes equation, Eq. (1), see also Multi-scale modeling: numerical implementation section for details. From the velocity and force distributions, $\mathbf {v}_1(\mathbf {x})$ and $\mathbf {f}_1(\mathbf {x})$, we can compute a generalized hydrodynamic friction coefficient, $\Gamma _{11}$, which is proportional to the phase-dependent rate of energy dissipation in the surrounding fluid. In the general case of n generalized coordinates $q_1,\ldots ,q_n$, we obtain a $n\times n$ matrix $\Gamma _{ij}$, see also Eq. (15). c Flow field induced by the active shape change of the cilium, here shown as two-dimensional section at $x=0$. The color represents the magnitude $|\mathbf {u}(\mathbf {x})|$ of three-dimensional velocity vectors, whereas white arrows represent the projections of $\mathbf {u}$ on the yz-plane. The flow field was computed as convolution of the fundamental solution of the Stokes equation with the force distribution $\mathbf {f}_1(\mathbf {x})$. Cilium phase corresponding to Fig. 3a: $\varphi _1 = 1.4\,\pi$, cilium beat frequency: $\omega _0 /(2\pi ) = 32 \,\mathrm {Hz}$ [12], dynamic viscosity of fluid: $\mu = 10^{-3}\, \text {Pa}\,\text {s}$ (corresponding to viscosity of water at $20\,^\circ \mathrm {C}$)

3 Generalized hydrodynamic friction forces

We introduce generalized hydrodynamic friction forces $P_i$ conjugate to the generalized coordinates $q_i$, following the convention of Lagrangian dynamics of dissipative systems [42], see also [43, 46]

$$\begin{aligned} P_i = \int _{\mathcal {S}} \!\mathrm{d}^2\mathbf {x}\, \mathbf {w}_i(\mathbf {x}) \cdot \mathbf {f}(\mathbf {x}) ,\quad i=1,\ldots ,n . \end{aligned}$$

(12)

The superposition principle for the shape changes $\mathbf {w}_i(\mathbf {x})$, Eq. (10), allows us to rewrite the total hydrodynamic dissipation rate $\mathcal {R}^{(h)}$ as a sum of products of generalized velocities times their conjugate generalized friction force

$$\begin{aligned} \mathcal {R}^{(h)} = \sum _i P_i\, \dot{q}_i. \end{aligned}$$

(13)

Note that the different generalized coordinates $q_i$ may have different physical units, in which case also all derived quantities will have different units; nonetheless, all vector and matrix operations of the formalism are consistent unit-wise.

In the special case, where some of the $q_i$ denote a rigid body transformation of an immersed microswimmer, i.e., a rigid body translation or rotation, the conjugate generalized force simply corresponds to the respective components of the total force $\mathbf {F}$ or total torque $\mathbf {M}$ exerted by the swimmer on the fluid, respectively, see the section on rigid body motion below.

Generalized hydrodynamic friction coefficients Using the superposition principle of Stokes flow, we can conveniently express the generalized hydrodynamic friction forces as linear function of the generalized velocities $\dot{q}_i$ by introducing generalized hydrodynamic friction coefficients

$$\begin{aligned} P_i = \sum _{j=1}^n \Gamma _{ij}\, \dot{q}_j ,\quad i=1,\ldots ,n . \end{aligned}$$

(14)

The generalized friction coefficients can be computed as scalar products between the (normalized) velocity profiles $\mathbf {w}_i(\mathbf {x})$, and the (normalized) force profiles $\mathbf {g}_j(\mathbf {x})$, see also Fig. 2

$$\begin{aligned} \Gamma _{ij} = \int _{\mathcal {S}} \! \mathrm{d}^2\mathbf {x}\, \mathbf {w}_i(\mathbf {x}) \cdot \mathbf {g}_j(\mathbf {x}) ,\quad i,j=1,\ldots ,n . \end{aligned}$$

(15)

Alternatively, we could express $\Gamma _{ij}$ in terms of partial derivatives with respect to the generalized velocities $\dot{q}_i$ as $\Gamma _{ij} = \int _{\mathcal {S}} \! \mathrm{d}^2\mathbf {x}\, (\partial \,\mathbf {v}(\mathbf {x})/\partial \,\dot{q}_i) \cdot (\partial \,\mathbf {f}(\mathbf {x})/\partial \,\dot{q}_j)$. We refer to diagonal elements $\Gamma _{ii}$ of the generalized hydrodynamic friction matrix $\varvec{\Gamma }$ as self-friction coefficients. Off-diagonal elements $\Gamma _{ij}$, $i\ne j$, or cross-friction coefficients, characterize a coupling between different degrees of freedom (e.g., a coupling between translational and rotational degrees of freedom for chiral objects or direct hydrodynamic interactions between different sub-objects that can, in principle, move independently).

The rate of hydrodynamic dissipation can thus be expressed as a quadratic form in the generalized velocity $\dot{\mathbf {q}}$

$$\begin{aligned} \mathcal {R}^{(h)} = \dot{\mathbf {q}}\cdot \varvec{\Gamma }\cdot \dot{\mathbf {q}} = \sum _{i,j} \Gamma _{ij} \, \dot{q}_i\dot{q}_j . \end{aligned}$$

(16)

The hydrodynamic dissipation rate $\mathcal {R}^{(h)}$ plays the role of a Rayleigh dissipation function for Lagrangian mechanics of dissipative systems [42]. Specifically, we could have equivalently defined the generalized forces as $2P_i = \partial \mathcal {R}^{(h)} / \partial \dot{q}_i$. (Following standard notation, the Rayleigh dissipation function is actually $\mathcal {R}^{(h)}/2$ [42].)

The matrix $\varvec{\Gamma }$ is symmetric, which represents a special case of Onsager reciprocity [56]. The proof follows directly from the Lorentz reciprocal theorem [27]: Let $\mathbf {v}_i$, $\varvec{\sigma }_i$ and $\mathbf {v}_j$, $\varvec{\sigma }_j$ denote the flow field and stress tensor associated with a change of only $q_i$ with rate $\dot{q}_i$, or a change of only $q_j$ with rate $\dot{q}_j$, respectively, while all other generalized coordinates are kept constant; then,

$$\begin{aligned} \dot{q}_i\, \Gamma _{ij}(\mathbf {q})\, \dot{q}_j&{=}&-\int _\mathcal {S}\! \mathrm{d}^2\mathbf {x}\, \mathbf {v}_i\cdot \varvec{\sigma }_j\cdot \mathbf {n}\nonumber \\&{\mathop {=}\limits ^{(*)}} -\int _\mathcal {S}\! \mathrm{d}^2\mathbf {x}\, \mathbf {v}_j\cdot \varvec{\sigma }_i\cdot \mathbf {n}= \dot{q}_j\, \Gamma _{ji}(\mathbf {q})\, \dot{q}_i, \end{aligned}$$

(17)

where we used the Lorentz reciprocal theorem at $(*)$.

The matrix $\varvec{\Gamma }$ is also positive semi-definite, consistent with the fact that the rate of energy dissipation should be nonnegative. (In fact, $\varvec{\Gamma }$ should be positive definite, except maybe at singular points $\mathbf {q}$ in configuration space, where $\mathbf {w}_i=\partial \mathbf {x}/\partial q_i$, $i=1,\ldots ,n$ are linearly dependent.)

In addition to hydrodynamic dissipation as characterized by $\mathcal {R}^{(h)}$, internal dissipative processes can be included in our framework, provided the corresponding dissipation function is likewise a quadratic form of the generalized velocity [10, 48].

4 Equation of motion

Balance of generalized forces We introduce active driving forces $Q_i$, $i=1,\ldots ,n$ that coarse-grain internal processes that drive the active shape changes of the active surface. Previous minimal models of flagella synchronization considered spheres moving along circular orbits driven by a tangential force [44, 57,58,59]. Our active driving forces $Q_i$ generalize the active driving forces considered in these models.

We postulate a balance of generalized forces between driving forces and hydrodynamic friction forces

$$\begin{aligned} Q_i = P_i,\ i=1,\ldots ,n . \end{aligned}$$

(18)

We emphasize that Eq. (18) is simply an instance of Newton’s second law and thus does not involve any new assumptions. Simplifying modeling assumptions have only been made in constraining the shape dynamics to a finite number of degrees of freedom $q_1,\ldots ,q_n$ and in the choice of the active driving forces $Q_1,\ldots ,Q_n$.

From the force balance equation, Eq. (18), and Eq. (12) expressing the generalized forces $P_i$, we obtain equations of motion for the generalized velocities $\dot{\mathbf {q}}$

$$\begin{aligned} \dot{\mathbf {q}} = \varvec{\Gamma }^{-1}\cdot \mathbf {Q}, \end{aligned}$$

(19)

where $\mathbf {Q}= (Q_1,\ldots ,Q_n)^T$ is the vector of active forces.

Calibration of active driving forces Because each driving force $Q_i$ characterizes internal processes, it is plausible to assume that $Q_i$ only depends on the corresponding degree of freedom $q_i$, but not on the other $q_j$, $j\ne i$, i.e., we may assume $Q_i=Q_i(\varphi _i)$. This assumption will hold in particular in applications, where the index i enumerates different microswimmers or different cilia. In principle, $Q_i$ may additionally depend on the friction force $P_i$ itself, i.e., if the internal active processes may change under load [17]. In this case, Eq. (18) becomes a self-consistency equation that has to be solved using methods for implicit equations. For a number of biological application cases, it was sufficient to assume that $Q_i$ is independent of load [10, 45, 48]. In this case, the active driving forces can be uniquely calibrated from a reference dynamics, ideally known from experiments. Once this is done, one can extrapolate to alternative dynamic scenarios.

As an example for this calibration procedure, previous work used experimental data of in-phase synchronized beating in the biflagellate green alga Chlamydomonas, which allowed to predict the response to perturbations of this synchronized state [45]. Similarly, measured cilia beat patterns in the absence of external flow have been used to calibrate active driving forces and predict the response to external flow [10]. In section “Application: pair of interacting cilia”, we consider the dynamics of an isolated cilium with constant phase speed to calibrate its active driving force. We then use this model to predict synchronization dynamics for a pair of cilia. In all these cases, the driving forces $Q_i$ coarse-grain internal active processes.

Additionally, the formalism allows to incorporate internal elastic degrees of freedom $q_i$ and the corresponding elastic restoring forces $Q_i$ in a formally equivalent manner. An example includes the waveform compliance of flagellar bending waves [10, 48]. Similarly, one can include external forces acting on self-propelled shape-changing microswimmers, as discussed in the next section.

5 Rigid body motion of a self-propelled microswimmer

The above formalism includes the important application case of shape-changing microswimmers and their self-propulsion in a viscous fluid. For that aim, we introduce rigid body transformation and include these in the set of generalized coordinates.

Specifically, we consider a microswimmer with outer surface $\mathcal {S}$ and introduce a material frame of this microswimmer consisting of a reference point $\mathbf {x}_0$ and a set of orthonormal vectors $\mathbf {e}_1$, $\mathbf {e}_2$, $\mathbf {e}_3$.

A rigid body motion of the swimmer is characterized by a translation of its reference point, $\dot{\mathbf {x}}_0 = \mathbf {v}_0 = v_1\,\mathbf {e}_1 + v_2\,\mathbf {e}_2 + v_3\,\mathbf {e}_3$, and a rotation of its material frame with $\dot{\mathbf {e}}_k = \varepsilon _{ijk} \Omega _j \mathbf {e}_i$, where $\varepsilon _{ijk}$ denotes the Levi-Cevita symbol and we use Einstein summation convention. The components $v_1$, $v_2$, $v_3$ and $\Omega _1$, $\Omega _2$ and $\Omega _3$ of the translational and the rotational velocity vector with respect to the basis $\mathbf {e}_1$, $\mathbf {e}_2$, $\mathbf {e}_3$, respectively, represent the six degrees of freedom of rigid body motion and satisfy $v_1=\mathbf {v}_0\cdot \mathbf {e}_1$, $v_2=\mathbf {v}_0\cdot \mathbf {e}_2$, $v_3=\mathbf {v}_0\cdot \mathbf {e}_3$, and $\Omega _1=\dot{\mathbf {e}}_2\cdot {\mathbf {e}}_3=-\dot{\mathbf {e}}_3\cdot {\mathbf {e}}_2$, $\Omega _2=\dot{\mathbf {e}}_3\cdot {\mathbf {e}}_1=-\dot{\mathbf {e}}_1\cdot {\mathbf {e}}_3$, $\Omega _3=\dot{\mathbf {e}}_1\cdot {\mathbf {e}}_2=-\dot{\mathbf {e}}_2\cdot {\mathbf {e}}_1$.

We choose these velocity components as the six generalized velocities

$$\begin{aligned} \dot{q}_1 = v_1,\ \dot{q}_2 = v_2,\ \dot{q}_3 = v_3,\ \dot{q}_4 = \Omega _1,\ \dot{q}_5 = \Omega _2,\ \dot{q}_6 = \Omega _3 . \end{aligned}$$

(20)

Formally, the coordinates $q_1, \ldots , q_6$ are elements of the Lie group $\mathfrak {se}(3)=\mathbb {R}^3\times \mathfrak {so}(3)$ of rigid body transformation [60].

For the special case, where the generalized velocities represent rigid body motion as in Eq. (20), the conjugate generalized hydrodynamic friction forces defined in Eq. (12) are simply given by the components of the total hydrodynamic friction force $\mathbf {F}$ and the total hydrodynamic friction torque $\mathbf {M}$, respectively

$$\begin{aligned}&P_1 = \mathbf {F}\cdot \mathbf {e}_1,\ P_2 = \mathbf {F}\cdot \mathbf {e}_2,\ P_3 = \mathbf {F}\cdot \mathbf {e}_3,\nonumber \\&P_4 = \mathbf {M}\cdot \mathbf {e}_1,\ P_5 = \mathbf {M}\cdot \mathbf {e}_2,\ P_6 = \mathbf {M}\cdot \mathbf {e}_3 . \end{aligned}$$

(21)

In this case, the $6\times 6$ matrix of generalized hydrodynamic friction coefficients $\varvec{\Gamma }$ reduces to the well-known hydrodynamic friction matrix (inverse mobility matrix) of an arbitrary-shaped rigid object. For a collection of rigid objects (e.g., a collection of rigid spheres as considered in [61]), we recover the inverse of the grand mobility matrix.

We can describe active shape changes of the microswimmer using coordinates $x'_1$, $x'_2$, $x'_3$ relative to the swimmer’s material frame for each point $\mathbf {x}\in \mathcal {S}$ on the surface

$$\begin{aligned} \mathbf {x}= \mathbf {x}_0 + x'_1\, \mathbf {e}_1 + x'_2\, \mathbf {e}_2 + x'_3\, \mathbf {e}_3 . \end{aligned}$$

(22)

We introduce the time-dependent rigid body transformation that maps the material frame of the swimmer to the laboratory frame, such that the reference point $\mathbf {r}_0$ of the swimmer is mapped to the origin $\varvec{0}\in \mathbb {R}^3$, and the material frame vectors are mapped to the standard unit vectors, respectively. The coordinates $x'_1$, $x'_2$, $x'_3$ are then just the coordinates of the image $\mathbf {x}'$ of a point $\mathbf {x}\in \mathcal {S}$ under this transformation, i.e., the coordinates of the surface after it has been brought into a reference condition [62]. Equation (22) allows us to decompose the displacement velocity $\mathbf {v}(x)$ of the surface into a contribution stemming from the rigid body motion and a contribution stemming only from any shape change

$$\begin{aligned} \mathbf {v}(\mathbf {x}) = \dot{\mathbf {x}} = \underbrace{ \dot{\mathbf {x}}_0 + \sum _{j=1}^3 x'_j\,\dot{\mathbf {e}}_j }_\text {rigid body motion} + \underbrace{ \sum _{j=1}^3 \dot{x}'_j\,\mathbf {e}_j }_\text {shape change} . \end{aligned}$$

(23)

The superposition principle of low-Reynolds number flow, Eq. (11), implies that the surface density $f(\mathbf {x})$ of hydrodynamic friction forces can be written as a superposition of contributions due to rigid body motion and a contribution $\mathbf {f}_\mathrm {act}(\mathbf {x})$ due the active shape change

$$\begin{aligned} \mathbf {f}(\mathbf {x})= & {} v_1 \,\mathbf {g}_1(\mathbf {x}) + v_2\, \mathbf {g}_2(\mathbf {x}) + v_3\, \mathbf {g}_3(\mathbf {x}) + \Omega _1\,\mathbf {g}_4(\mathbf {x}) \nonumber \\&+ \Omega _2\,\mathbf {g}_5(\mathbf {x}) + \Omega _3\,\mathbf {g}_6(\mathbf {x}) + \mathbf {f}_\mathrm {act}(\mathbf {x}) , \end{aligned}$$

(24)

where $f_\mathrm {act}(\mathbf {x})$ depends only on the shape change $\dot{\mathbf {x}}'$, but not on the translational velocity $\mathbf {v}_0$ nor the rotational velocity $\varvec{\Omega }$.

Since inertia is assumed negligible, the total force and total torque acting on a microswimmer must equal any external force or torque acting on the swimmer, $\mathbf {F}= \mathbf {F}^\mathrm {ext}$, $\mathbf {M}=\mathbf {M}^\mathrm {ext}$ [27]. It follows that a microswimmer that is free from external forces or torques does not exert any net force or torque on the surrounding fluid itself

$$\begin{aligned} \mathbf {F}= \mathbf {0} ,\quad \mathbf {M}= \mathbf {0} . \end{aligned}$$

(25)

Equation (25) holds in particular for a neutrally buoyant biological microswimmer (a good approximation for many biological microswimmers).

The surface density of hydrodynamic friction forces due to active shape changes, $\mathbf {f}_\mathrm {act}(\mathbf {x})$, gives rise to a contribution $\mathbf {F}_\mathrm {act}= \int _\mathcal {S}\! \mathrm{d}^2\mathbf {x}\, \mathbf {f}_\mathrm {act}(\mathbf {x})$ to the total force, as well as an analogous contribution $\mathbf {M}_\mathrm {act}$ to the total torque. The force and torque balance equations, Eq. (25), thus provide an inhomogeneous system of six linear equations for the six components of the translational and rotational velocity, $\mathbf {v}_0$ and $\varvec{\Omega }$.

We emphasize that Eq. (18) is very general, and includes the following application cases of microswimmer motion:

• External forces or torques: For example, external forces $\mathbf {F}^\mathrm {ext}$ or external torques $\mathbf {M}^\mathrm {ext}$ are captured by corresponding external forces $Q_i^\mathrm {ext}$. Examples include gravitational force for a non-buoyant swimmer or torques exerted by an external rotating magnetic field on an artificial microswimmer with nonzero magnetic dipole moment.

• Prescribed shape dynamics: For a prescribed shape dynamics, say of shape coordinate $q_i$ with prescribed driving protocol $q_i(t)$, one would omit the corresponding force balance equation $Q_i=P_i$ from the set of equations Eq. (18), and solve for the equation of motion of the other coordinates with prescribed $q_i(t)$. The conjugate hydrodynamic friction force $P_i$ nonetheless appears in the formula for the total hydrodynamic dissipation rate $\mathcal {R}^{(h)}$, where $P_i\dot{q}_i$ equals the rate of work required for the shape change with rate $\dot{q}_i$. A number of classical theory publications on self-propelled biological microswimmers considered prescribed shape dynamics [28, 49,50,51,52, 62].

• Constrained motion: Several applications considered constrained swimmers, for example, biological microswimmers clamped in micropipettes constrained from translational motion [19, 20, 22]. Formally, this is a special case of a coordinate $q_i$ with prescribed dynamics for the coordinates $q_1,\ldots ,q_3$ representing rigid body translation, enforcing $\dot{q}_i=0$. The conjugate hydrodynamic friction force $P_i$ equals the external constraining force required to impose the constraint. Similarly, to constrain a microswimmer from rotational motion requires a constraining torque $\mathbf {M}= P_4\,\mathbf {e}_1 + P_5\,\mathbf {e}_2 + P_6\,\mathbf {e}_3$. As a historical note, in their classical 1955 paper, Gray & Hancock considered self-propulsion of sperm cells with constrained rotational motion to simplify the calculation [28].

  Finally, clamped microswimmers exposed to uniform external flow with flow velocity $\mathbf {u}_0$ far from the swimmer as considered in [10] can be incorporated into our formalism by switching to a co-moving reference frame in which the fluid is at rest. In the co-moving frame, the clamped swimmer is “dragged” through the fluid, corresponding to a constraint for rigid body translation, $\dot{q}_i = -\mathbf {u}_0\cdot \mathbf {e}_i$, $i=1,2,3$. Correspondingly, the total hydrodynamic friction force $\mathbf {F}= P_1\,\mathbf {e}_1 + P_2\,\mathbf {e}_2 + P_3\,\mathbf {e}_3$ represents the constraining force required to clamp the microswimmer in such an external flow.

6 Multi-scale modeling: numerical implementation

To solve for the dynamics of an active surface according to Eq. (19), it suffices to compute the generalized hydrodynamic friction matrix $\varvec{\Gamma }$ for a set of reference configurations $\mathbf {q}$ and save this as a look-up table; the friction matrix $\varvec{\Gamma }(\mathbf {q})$ for arbitrary $\mathbf {q}$ can then be found by interpolation. This allows to solve the equation of motion Eq. (19) fast, using pre-computed hydrodynamic friction coefficients. We outline the numerical implementation of this general program.

While Eq. (15) may look abstract, all quantities can be directly obtained from numerical computations for arbitrary surfaces $\mathcal {S}$. Assume the surface $\mathcal {S}$ is represented by a triangulated mesh. The triangular faces (or “elements”) shall be enumerated by $k\in \mathcal {I}$ with midpoints $\mathbf {x}_k$ and respective areas $A_k$.

In a first step, we compute a (normalized) surface distribution of velocities $\mathbf {w}_i(\mathbf {x}_k)$, $k\in \mathcal {I}$ for each generalized coordinate $i=1,\ldots ,n$, either by computing the derivative $\mathbf {w}_i(\mathbf {x}_k) = \partial \mathbf {x}_k(\mathbf {q}) / \partial q_i$ analytically or by evaluating the finite difference quotient

$$\begin{aligned} \mathbf {w}_i(\mathbf {x}_k) = \frac{ \mathbf {x}_k (\mathbf {q}+ \varepsilon \,\Delta _i) - \mathbf {x}_k(\mathbf {q}) }{\varepsilon } , \end{aligned}$$

(26)

for each midpoint $\mathbf {x}_k$, $k\in \mathcal {I}$, where $\Delta _i$ is the unit vector whose components are all zero, except the ith component, and $\varepsilon$ is a small number.

We can use boundary element methods to numerically compute a surface density of hydrodynamic friction forces $\mathbf {f}(\mathbf {x}_k)$ with physical units of a stress, given an arbitrary surface distribution of velocities $\mathbf {v}(\mathbf {x}_k)$ specified at each midpoint $\mathbf {x}_k$, $k\in \mathcal {I}$. Specifically, in the application example below, we use the fast multi-pole boundary element method fastBEM [39, 40].

In the next step, we compute the surface density $\mathbf {f}_j(\mathbf {x}_k) = \alpha _j\,\mathbf {g}_j(\mathbf {x}_k)$ of hydrodynamic friction forces, corresponding to the velocity distribution $\mathbf {v}_j(\mathbf {x}_k)=\alpha _j\,\mathbf {w}_j(\mathbf {x}_k)$. Here, $\alpha _j$ is an arbitrary constant to ensure proper physical units of a velocity for $\mathbf {v}_j$. We thus obtain n surface distributions of (normalized) hydrodynamic friction forces $\mathbf {g}_j(\mathbf {x}_k)$, $j=1,\ldots ,n$, one for each generalized coordinate $q_j$. These force distributions $\mathbf {g}_j(\mathbf {x}_k)$ depend on $\mathbf {q}$, but not on $\dot{\mathbf {q}}$. Finally, we compute the components $\Gamma _{ij}$ of the generalized hydrodynamic friction matrix $\varvec{\Gamma }$ by taking the scalar product of the ith (normalized) velocity distribution $\mathbf {w}_i(\mathbf {x}_k)$, and the jth (normalized) force distribution $\mathbf {g}_j(\mathbf {x}_k)$

$$\begin{aligned} \Gamma _{ij} = \sum _{k\in \mathcal {I}} \mathbf {w}_i(\mathbf {x}_k)\cdot \mathbf {g}_j(\mathbf {x}_k)\, A_k ,\quad i,j=1,\ldots ,n , \end{aligned}$$

(27)

where $A_k$ was the area of the kth triangle. We can interpret $\alpha _j\mathbf {g}_j(\mathbf {x}_k)A_k$ at the total force acting on the kth element (with proper physical units of a force) if the generalized coordinate $q_i$ would change at a rate $\alpha _i$.

Importantly, it suffices to compute the generalized hydrodynamic friction matrix $\varvec{\Gamma }$ only for a set of reference configurations and save this as a look-up table. If m discrete values are used for each of the n generalized coordinates, the Stokes equation needs to be solved a total of $n\,m^n$ times, as we need to change each of the $q_j$, $j=1,\ldots ,n$ for $m^n$ different choices of $\mathbf {q}$. By exploiting symmetries, as well as translational and rotational invariance for individual microswimmers, this number can be reduced further. The friction matrix $\varvec{\Gamma }(\mathbf {q})$ for arbitrary $\mathbf {q}$ can then be found by interpolation. For example, cubic spline interpolation, low-order polynomials, and (double) Fourier series were used in previous applications [10, 45, 48].

In principle, different hydrodynamic simulation methods could be used to solve the Stokes equation and compute the force distribution $\mathbf {f}(\mathbf {x}_k)$. Deterministic lattice Boltzmann solvers may be suitable, provided the effective Reynolds numbers are sufficiently small. An early application represented the surface of a microswimmer not by a triangulated mesh, but as a collection of equally sized spheres, and computed the grand mobility matrix for these spheres using the hydrolib package [63]. In the application example below, we employ the fast multi-pole boundary element method fastBEM [39, 40], available for download at [64]. The open-source implementation of the fast boundary element method STKFMM directly incorporates the fundamental solution of the Stokes equation close to a no-slip boundary wall [65] and thus relieves the need for an explicit representation of the boundary as a triangulated mesh, yet currently only supports the computation of velocity fields from force distributions [66, 67].

7 Application: pair of interacting cilia

We demonstrate our LAMAS modeling framework using the example of hydrodynamic synchronization in pairs of cilia. We therefore reconsider the question of in-phase and anti-phase synchronization previously addressed by Vilfan et al. [57], yet, instead of a minimal model of spheres orbiting on circular trajectories, we employ in our simulations a realistic cilia beat pattern obtained from previous experiments.

We digitalized and reconstructed three-dimensional shapes of a beating cilium on the surface of the unicellular ciliated protist Paramecium [12] as presented in [68]. The cilia beat is periodic, and we can thus describe the shape of the cilia centerline as a periodic shape sequence parameterized by a $2\pi$-periodic phase variable $\varphi$, see Fig. 3a. For unperturbed beating, the phase speed equals the angular frequency of the cilia beat, $\dot{\varphi }(t) = \omega _0$.

Fig. 6

In-phase and anti-phase synchronization in a pair of interacting cilia. a Cilia beat pattern from unicellular Paramecium [12] as reported in [68], shown as sequence of three-dimensional shapes parameterized by a $2\pi$-periodic phase variable $\varphi$ (color code). Spacing of square grid: $2\,\upmu \mathrm {m}$. b We consider a pair of cilia with respective phases $\varphi _1$ and $\varphi _2$, whose base points are separated by a distance d along a direction that encloses an angle $\psi$ with the x axis (where the y axis is set by the direction of the effective stroke of both cilia). c Self-friction coefficient $\Gamma _{11}(\varphi _1)$ of a single cilium as function of its phase variable $\varphi _1$, obtained by solving the Stokes equation of three-dimensional flow (blue dots), as well as continuous representation as Fourier series (orange line). In the case of a single cilium, $\Gamma _{11}$ is proportional to the phase-dependent active cilia driving force $Q_1(\varphi _1)$. d Generalized hydrodynamic friction coefficient $\Gamma _{12}(\varphi _1,\varphi _2)$ characterizing hydrodynamic interactions from the second cilium to the first cilium, see also Eq. (28). Positive values (red colors) imply that the motion of the second cilium causes the first cilium to beat slower, while negative values (blue colors) imply that the first cilium beats faster. Cilia distance $d = 18\,\upmu \mathrm {m}$, orientation angle $\psi = 5 \pi / 6$. e The magnitude of hydrodynamic interactions, here quantified by the $L_2$-norm of $\Gamma _{12}$, decay as $\sim 1/d^3$, consistent with the theoretical scaling expected from the Blake tensor [65]. Different curves correspond to different separation directions between the two cilia ($\psi =0$: dark blue, $\psi =\pi /3$: light green, $\psi =2\pi /3$: teal; also indicated by the direction arrows.) f We characterize the stability of the in-phase synchronized state, defined by $\varphi _1(t)=\varphi _2(t)$, for different relative orientations of the two cilia by a Lyapunov exponent $\lambda$, see Eq. (30). Colored dots at respective positions in the xy-plane represent the value of $\lambda$ if the second cilium is positioned at the position of the dot and the first cilium is located at the origin. Negative values imply that in-phase synchronization is linearly stable (green colors, $\lambda <0$), while positive values imply that in-phase synchronization is linearly unstable (red colors, $\lambda >0$). g We determined the steady-state phase difference $\delta ^*$ between the two cilia for different relative cilia positions, analogous to panel F. While $\delta ^*=0$ for cilia orientations with stable in-phase synchronization (cyan), we observe anti-phase synchronization with $\delta ^*\approx \pi$ for cilia orientations with $\lambda >0$ (red colors). For relative cilia orientation aligned with the direction of the effective stroke of the cilia beat ($\psi =\pi /2$), we observed cases of multi-stability (bicolored dots). h Consistent with the far-field scaling of hydrodynamic interactions as shown in panel E, we find that also the Lyapunov exponent $\lambda$, which represents an effective synchronization strength, likewise decays as $1/d^3$ with distance d between the two cilia. Different curves represent different separation directions, analogous to panel E. Frequency of cilia beat: $\omega _0 /(2\pi ) = 32 \,\mathrm {Hz}$ [12], fluid viscosity: $\mu = 10^{-3}\, \text {Pa}\,\text {s}$

7.1 Equation of motion for a pair of cilia

We consider two identical cilia beating in the same direction attached to a no-slip boundary wall, see Fig. 3a and b. We describe each cilium by a single phase variable that parameterizes its periodic sequence of centerline shapes. The two phase variables $\varphi _1$ and $\varphi _2$ fully characterize the dynamics of the two beating cilia, and represent a set of generalized coordinates with state vector $\mathbf {q}=(\varphi _1,\varphi _2)$.

For our example, the force balance equation, Eq. (18), takes the form

$$\begin{aligned} Q_1(\varphi _1)&= \Gamma _{11}(\varphi _1, \varphi _2) \dot{\varphi }_1 + \Gamma _{12}(\varphi _1, \varphi _2) \dot{\varphi }_2 \nonumber \\ Q_2(\varphi _2)&= \Gamma _{21}(\varphi _1, \varphi _2) \dot{\varphi }_1 + \Gamma _{22}(\varphi _2, \varphi _2) \dot{\varphi }_2 . \end{aligned}$$

(28)

This equation can be further simplified. The symmetry relation Eq. (17) implies $\Gamma _{12}(\varphi _1,\varphi _2) = \Gamma _{21}(\varphi _1,\varphi _2)$ for the hydrodynamic interactions. Numerical computation of $\Gamma _{11}(\varphi _1, \varphi _2)$ shows that this self-friction coefficient of the first cilium is virtually independent of the phase of the second cilium, and almost does not change when the other cilium is not present at all. An analogous statement holds for the second cilium. Therefore, we can replace the two self-friction coefficients in Eq. (28), $\Gamma _{11}(\varphi _1, \varphi _2)$ and $\Gamma _{22}(\varphi _1, \varphi _2)$, by the self-friction coefficient for a single cilium to very good approximation. This approximation allows us to define the active driving forces using the case of a single cilium.

Calibration of active driving force We require that a single cilium should beat at a constant phase speed $\dot{\varphi }_1 = \omega _0$, where $\omega _0$ denotes the intrinsic beat frequency of the cilium if there are no interactions with other cilia. This requirement uniquely determines the active driving force $Q_1(\varphi _1)$. Specifically, for a single cilium, the force balance equation reads, $Q_1(\varphi _1) = \Gamma _{11}(\varphi _1)\,\dot{\varphi }_1$. We conclude $Q_1(\varphi _1) = \omega _0\, \Gamma _{11}(\varphi _1)$; Fig. 3c displays the phase dependence of $\Gamma _{11}(\varphi _1)$. Since both cilia are assumed identical with same intrinsic beat frequency $\omega _0$, this also specifies the active driving force $Q_2(\varphi _2)$ of the second cilium.

Equation of motion Using the force balance equation, Eq. 28, and the calibrated driving force, we obtain the equation of motion

$$\begin{aligned} \dot{\varphi _1}&= \omega _0 - C_1(\varphi _1,\varphi _2) \,\dot{\varphi }_2 ,\quad C_1(\varphi _1,\varphi _2) = \frac{\Gamma _{12}(\varphi _1, \varphi _2)}{\Gamma _{11}(\varphi _1)} \nonumber \\ \dot{\varphi _2}&= \omega _0 - C_2(\varphi _1,\varphi _2) \,\dot{\varphi }_1 ,\quad C_2(\varphi _1,\varphi _2) = \frac{\Gamma _{12}(\varphi _1, \varphi _2)}{\Gamma _{11}(\varphi _2)} . \end{aligned}$$

(29)

Equation (29) describes a pair of coupled phase oscillators.

In the following, we use Eq. (29) and pre-computed friction coefficients to analyze in-phase and anti-phase synchronization of the two cilia depending on their relative position. Details on the numerical computation of $\Gamma _{ij}(\varphi _1,\varphi _2)$ can be found in the appendix. An example of the generalized friction coefficient $\Gamma _{12}(\varphi _1, \varphi _2)$, which characterizes hydrodynamic interactions between the two cilia, is shown in Fig. 3d; Fig. 5 shows $\Gamma _{12}(\varphi _1,\varphi _2)$ for additional cilia orientations.

7.2 Results: in-phase and anti-phase synchronization as function of direction

Hydrodynamic interactions decay as $1/d^3$. For large separation distances d between the two cilia, hydrodynamic interactions between the two cilia as characterized by $\Gamma _{12}(\varphi _1,\varphi _2)$ decay as $1/d^3$, see Fig. 3e. This asymptotic scaling is consistent with the expected leading-order singularity of the flow field induced by a single cilium. Specifically, the flow field induced by a point force parallel to a no-slip boundary wall is given by the Blake tensor [65] and decays as $1/d^3$ for points at a constant height from the boundary, which is the relevant case for the interaction between cilia [69].

Linear stability analysis Since both cilia were assumed identical, the in-phase synchronized state with $\varphi _1(t)=\varphi _2(t)$ is always a periodic solution of Eq. (29). To assess the linear stability of this in-phase synchronized state, we monitored the evolution of a small perturbation of the phase difference $\delta (t) = \varphi _2(t)-\varphi _1(t)$ during one beat cycle. Specifically, we integrated Eq. (29) with the initial condition $\varphi _1(t=0)=-\delta _0/2$ and $\varphi _2(t=0)=\delta _0/2$ for a small perturbation $|\delta _0|\ll 1$ up to time T defined by $\overline{\varphi }(T) = [\varphi _1(T)+\varphi _2(T)]/2=2\pi$ (corresponding to the completion of a full beat cycle), and recorded the new phase difference $\delta _1=\delta (T)$.

We define a dimensionless Lyapunov exponent as

$$\begin{aligned} \lambda = \log |\delta _1 / \delta _0|, \end{aligned}$$

(30)

which characterizes whether the initial perturbation decays or grows. The in-phase synchronized state is linearly stable if $|\delta _1| < |\delta _0|$ (hence $\lambda < 0$) and linearly unstable if $|\delta _1| > |\delta _0|$ (hence $\lambda > 0$).

Figure 3f shows $\lambda$ as function of relative cilia position. Here, the first cilium is located at the origin, while the second cilium is located at the position of the respective colored dots.

The symmetry of Eq. (29) implies that the synchronization behavior is invariant under a point reflection, which swaps cilia 1 and 2. Whether in-phase synchronization is stable or not only depends on the direction of the separation vector between the two cilia (where $\lambda >0$ for direction angles $\psi =2\pi /3$ and $5\pi /6$, in which case the cilia synchronize anti-phase, as discussed next).

Additionally, we analyzed the steady-state dynamics of Eq. (29) and identified phase differences $\delta ^*$ that correspond to stable periodic solutions, see Fig. 3g. As a technical point, $\delta (t)$ may weakly oscillate during each cycle; we therefore define $\delta ^*$ as the phase difference at times for which $\overline{\varphi }$ is an integer multiple of $2\pi$.

When the in-phase synchronized state is linearly stable for a given cilia configuration, we obviously have $\delta ^*=0$. If, however, the in-phase synchronized state is linearly unstable, we approximately find $\delta ^*\approx \pi$, corresponding to anti-phase synchronization. For a few cilia configurations, we observe multi-stability, characterized by two different values of the phase differences $\delta ^*$ that correspond to stable periodic solutions; these configurations are indicated as bicolored half circles in Fig. 3g.

The magnitude $|\lambda |$ of the Lyapunov exponents decreases as $1/d^3$ with distance d between the two cilia, see Fig. 3h, consistent with the far-field scaling of hydrodynamic interactions shown in panel E. This suggests that short-range interactions between close-by cilia may dominate emergent behavior in carpets of many cilia.

8 Discussion

Summary We presented a multi-scale modeling and simulation framework for active surfaces immersed in viscous fluids. This includes self-propulsion of shape-changing microswimmers as a special case. The key idea is to constrain the shape dynamics to a small number of principal deformation modes. These modes represent generalized coordinates, for which generalized hydrodynamic friction coefficients are defined according to the formalism of Lagrangian mechanics. To actually compute these friction coefficients, the Stokes equation is solved for an infinitesimal change of each generalized coordinate in an initial step. For subsequent dynamic simulations, a generalized force balance between hydrodynamic friction forces and active driving forces is solved in each time step. This is sufficiently fast since this second step does not involve any hydrodynamic computations, but uses the pre-computed hydrodynamic friction coefficients.

Our formalism generalizes classical Lagrangian dynamics of dissipative systems [42] to active systems. The rate of work exerted by the active surface on the surrounding fluid provides a Rayleigh dissipation function $\mathcal {R}^{(h)}$, which defines generalized friction forces $P_i$ conjugate to each generalized coordinate $q_i$ as a partial derivative $2P_i = \partial \mathcal {R}^{(h)}/\partial q_i$. Numerically, the generalized friction forces are computed from a surface density of hydrodynamic friction forces using a Lagrangian projection method. Active driving forces coarse-grain internal active processes, such as the dynamics of molecular motors inside cilia and flagella. These active driving forces can be calibrated from a reference data set, for which the dynamics is already known or prescribed.

Our approach shares the idea of multi-scale modeling with recent developments of reduced-order models, which likewise propose a decomposition of biological hydrodynamics problems with multiple queries into an initial setup phase during which the Stokes equation needs to be solved for example configurations (“offline phase”), and an subsequent phase of parameter space exploration (“online phase”) [41]. However, our approach does not require an affine dependence of hydrodynamic quantities on model parameters.

Potential applications We applied our general framework to a model example of mutual synchronization between two cilia, using an experimentally measured cilia beat pattern. Future work will generalize this approach to cilia carpets with many cilia [70], which previously had been either studied using detailed simulations with many degrees of freedom [71, 72], or using minimal models [69, 73,74,75]. A key simplifying assumption will be to approximate many-body hydrodynamic interactions between many cilia as a superposition of pairwise interactions. A similar approach can be applied to study self-organized pattern formation in suspension of shape-changing microswimmers, using the approximation of pairwise interactions between microswimmers, which is valid for dilute suspensions.

An important feature of our modeling framework is that biological noise can be incorporated in a natural way. Beating cilia are known to exhibit both phase fluctuations (frequency jitter), as well as amplitude fluctuations [20, 53, 76]. This active noise jeopardizes synchronization of cilia by hydrodynamic interactions. Additionally, noise randomizes the motion of biological microswimmers. While thermal noise causes noticeable rotational diffusion of micrometer-sized bacteria such as E. coli [77], amplitude fluctuations of cilia beating affect the swimming of tenfold larger eukaryotic swimmers [47]. In our framework, active noise is incorporated by using stochastic active driving forces. In previous work, adding additive Gaussian white noise with noise strengths calibrated from experiments was sufficient to account for effective diffusion of swimming sperm cells, or noisy synchronization of coupled cilia [53]. For simulations accounting for biological noise, it is beneficial to use a deterministic solver for the Stokes equation as done here, in order to not confound physically relevant noise and fluctuations from a stochastic hydrodynamic simulation method.

Next, our modeling framework helps to conceptualize the load response of cilia and flagella, which beat slower if the hydrodynamic load opposing their beat increases. Classical work showed how cilia and flagella beat slower at increased viscosity of the surrounding fluid [12, 14]; likewise, external flows change the speed of cilia beating [10, 17, 23]. The load response of cilia and flagella is a prerequisite for putative mechanisms of synchronization by hydrodynamic or mechanical interactions. We propose that the generalized hydrodynamic friction force defined here can serve as a proxy for the effective hydrodynamic load acting on an actively shape-changing structure such as a beating cilium.

Limitations Our approach is restricted to the limit of zero Reynolds numbers, because it crucially relies on the superposition principle for Stokes flow. In a laminar flow regime at finite Reynolds numbers, we expect that computations of self-friction are still accurate, but long-range hydrodynamic interactions become increasingly less accurate with increasing distance if the Stokes equation is used. Nonetheless, our approach should still serve as a reasonable approximation, since any long-range hydrodynamic interactions that are incorrectly predicted by the Stokes equation will be very weak already.

In principle, a similar framework could be developed using the linearized Navier–Stokes equations instead of the Stokes equation used here, but only in Fourier space. The linearized Navier–Stokes equation provides a refined approximation for long-range hydrodynamic interactions if the Reynolds number for oscillatory motion becomes appreciable. In this case, a superposition principle applies for time-periodic flows. However, working in frequency space instead of the time domain will make the practical solution of dynamic problems more difficult.

Another limitation of our approach is that it is inherently restricted to Newtonian fluids. While certain important biological fluid dynamics problems involve viscoelastic fluids, the lack of a superposition principle in this case implies that other methods need to be used.

Conclusion Our modeling and simulation framework LAMAS can be complimentary to existing methods. Our approach is particularly suited to screen extensive parameter ranges, provided the modified parameters concern the dynamical model (such as active driving forces or effective elastic constants [48]), and do not require re-computation of the generalized hydrodynamic friction coefficients. Likewise, our approach allows to compute multiple stochastic realizations of the same problem fast.

Appendix: numerical methods

Fig. 7

Triangulated mesh for pair of cilia attached to boundary surface. a Entire mesh consisting, of two cilia represented as bent cylinders (red, cyan), as well as a thin disk of radius $60\,\upmu \mathrm {m}$ representing the boundary surface (gray), corresponding to approximately 1000 triangular elements per cilium and 7500 elements for the boundary surface. b Close-up view on a single cilium. Nodes on the bottom of the surface are hidden from view. Cilia distance $d = 18\,\upmu \mathrm {m}$, orientation angle $\psi = 2\pi /3$

Fig. 8

Hydrodynamic interaction as function of cilia phases for different cilia orientations. Generalized hydrodynamic friction coefficient $\Gamma _{12}(\varphi _1, \varphi _2)$ as in Fig. 3d for different cilia orientation angles: left: $\psi = \pi /2$ (direction of effective stroke), middle: $\psi = 2\pi /3$ (oblique to direction of effective stroke), right: $\psi = \pi$ (perpendicular to direction of effective stroke). Cilia distance: $d = 18\,\upmu \mathrm {m}$. Note the different color scale compared to Fig. 3d

We present additional details on the numerical computations for the application case of a pair of interacting cilia.

Mesh generation We generated a triangulated mesh for the combined surface $\mathcal {S}$ of cilia and boundary surface using a custom-built Python package (available upon request), see also Fig. 4. We represent the digitalized shapes of the cilia centerline as a family of space curves $\mathbf {r}(s,\varphi )$ parameterized by arclength s with $0\le s\le L$, where $L=10\,\upmu \mathrm {m}$ is the length of the cilium, and a $2\pi$-periodic phase variable $\varphi$. The centerline shapes of the two cilia are thus given by $\mathbf {r}_{0,1} + \mathbf {r}(s,\varphi _1)$ and $\mathbf {r}_{0,2} + \mathbf {r}(s,\varphi _2)$, where the base points $\mathbf {r}_{0,1}$ and $\mathbf {r}_{0,2}$ have a distance d. The separation vector $\mathbf {r}_{0,2}-\mathbf {r}_{0,1}$ encloses an angle $\psi$ with x-axis (where y axis is set by the effective stroke of both cilia), see Fig. 3b.

We generate a triangulated mesh for the each cilium by treating the cilium as a bent cylinder of radius $0.125\,\upmu \mathrm {m}$, using 8 node points in azimuthal direction, and 61 nodes in longitudinal direction, as well as one apex node at the proximal and distal ends, respectively. For numerical stability, the proximal apices of each cilium mesh have a distance of $0.25\,\upmu \mathrm {m}$ from the boundary surface. A smaller distance virtually does not change the computed friction coefficients, but can cause convergence issues.

The hydrodynamic solver fastBEM requires closed surfaces, which prompted us to use a circular disk of finite thickness (radius $60\,\upmu \mathrm {m}$, thickness $1.5\,\upmu \mathrm {m}$) instead of a plane surface. Initial simulations confirmed that using a larger disk radius virtually did not change results. Disk faces were meshed using the Python triangle package (minimum triangle angle $20^\circ$, maximum triangle area $2\,\upmu \mathrm {m}^2$ on the upper surface of the disk up to a distance of $50\,\upmu \mathrm {m}$ from the disk center, $20\,\upmu \mathrm {m}^2$ otherwise). Additionally, to improve the convergence of the solver, we refined the mesh in a small area below the proximal apices of the cilia (maximum triangle area $0.04\,\upmu \mathrm {m}^2$ up to a distance of $0.625\,\upmu \mathrm {m}$ from cilia base points). In total, each meshed cilium consists of 975 triangles, while the meshed disk consists of approximately 7600 triangles.

Hydrodynamic computations We employed a fast multipole boundary element method termed fastBEM [39, 40, 64] to solve the inverse problem of finding the surface distributions of hydrodynamic friction forces $\mathbf {f}(\mathbf {x})$ for given velocity fields $\mathbf {v}(\mathbf {x})$ on the combined surface $\mathcal {S}$ of both cilia and the boundary surface. To solve this inverse problem, the algorithm employs an iterative linear GMRES solver (tolerance parameter used here, $\mathrm {tol}=5 \cdot 10^{-4}$). In principle, the solver would allow also for mixed boundary conditions that specify a combination of forces and velocities on different parts of the surface, which is, however, not needed here.

Initial tests showed that the self-friction is virtually independent of the phase the other cilium, allowing us to approximate $\Gamma _{11}(\varphi _1,\varphi _2) \approx \Gamma _{11}(\varphi _1)$ and $\Gamma _{22}(\varphi _1,\varphi _2) \approx \Gamma _{11}(\varphi _2)$, where $\Gamma _{11}(\varphi _1)$ corresponds to the simulation result for a single cilium. For the smallest distance tested here, $14\,\upmu \mathrm {m}$, the difference was at most $2\%$. Thus, computation of $\Gamma _{11}$ required m hydrodynamic computations for $m=20$ equidistant phase values. The symmetry relation Eq. (17) gives $\Gamma _{21}(\varphi _1,\varphi _2) = \Gamma _{12}(\varphi _1,\varphi _2)$; thus, it is sufficient to compute only $\Gamma _{12}$ (i.e., perform only computations where cilium number 2 moves, while cilium number 1 is static). To compute $\Gamma _{12}(\varphi _1, \varphi _2)$, we performed $m^2 = 400$ hydrodynamic computations for $m^2$ pairs of phase values on a equidistant $(\varphi _1,\varphi _2)$-grid, with mean CPU time of about $10^2$ seconds per computation. We repeated these computations for 42 different relative positions of cilia as shown in Fig. 3f.

Interpolation From the generalized friction coefficients computed for a discrete set of $(\varphi _1,\varphi _2)$-values, we obtained in a final step a continuous representation in the form of a (double) Fourier series truncated after order 4 (corresponding, e.g., to $(2 \cdot 4+1)^2=81$ Fourier terms for $\Gamma _{12}(\varphi _1,\varphi _2)$).

Dynamical equation The system of coupled ordinary differential equations, Eq. (28), was solved with Python (method scipy.integrate.solve_ivp, tolerance $10^{-8}$, scipy version 1.5.0). In each time step, we compute the inverse of matrix $\mathbf {\Gamma }$.

Lyapunov exponents For the computation of Lyapunov exponents $\lambda$ shown in Fig. 3f, we used a small perturbation $\delta _0 = 10^{-2}$ of the in-phase synchronized state. Preliminary simulations using a smaller perturbation $\delta _0 = 10^{-3}$ gave virtually identical results.

Steady-state phase difference We computed the steady-state phase difference $\delta ^*$ between the two cilia as a fixed point of the Poincaré map $\mathcal {L}:\delta _0\rightarrow \delta _1$. Specifically, we computed $\mathcal {L}(\delta _0)$ for 30 equidistant values of $\delta _0$ in the interval $[0,2\pi )$ by integrating Eq. (28) using initial conditions $\varphi _1(t=0)=-\delta _0/2$ and $\varphi _2(t=0)=+\delta _0/2$. We then numerically solved for fixed points $\mathcal {L}(\delta ^*) = \delta ^*$, using monotonic cubic spline interpolation of $\mathcal {L}$. The periodic solution corresponding to a steady-state phase difference $\delta ^*$ is stable if the numerical derivative $d \mathcal {L}/ d\delta _0|_{\delta _0=\delta ^*}$ is smaller than 1.