Lagrangian Mechanics of Active Systems

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 re-used 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.

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 propagates 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 spatio-temporal 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 steps 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 sub-problems 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 methods manuscript, 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. Thereby, our framework extends Lagrangian mechanics of dissipative systems to active surfaces and active microswimmers, whose shape dynamics is driven by active forces.

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] with incompressibility condition ∇ · u = 0. Here, u(x) denotes the flow velocity, p(x) the pressure field, and µ the dynamic viscosity of the fluid. The total stress tensor σ for an incompressible fluid depends on both the hydrostatic pressure p and the symmetrized strain rate tensor ∆ [27] Thus, the Stokes equation, Eq. (1) could be equivalently written as 0 = ∇ · σ in the bulk of the fluid. Special conditions apply at boundaries. No-slip boundary condition for an active surface. We consider a surface S immersed in the fluid that changes its shape as a function of time. For example, 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 v(x, t) for each point x ∈ S at time t.
We impose a no-slip boundary condition at this surface, i.e., require that the local velocity u(x) of fluid flow matches the local velocity v(x) of the surface for each surface point Hydrodynamic friction forces. A shape change of the surface S induces a flow field u(x) with corresponding stress tensor field σ(x). The stress σ(x) determines the surface density of forces f (x) exerted by the surface on the fluid (with units of a stress N/m 2 , also called contact force, or traction force density) where n is the surface normal pointing into the fluid. Correspondingly, −f (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 f (x) Analogously, the total torque (with respect to a reference point x 0 ) exerted by the surface on the fluid is given by 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 v that is given as a linear combination of velocity distributions v 1 and v 2 as with real coefficients α 1 , α 2 ∈ R. Then, the resultant flow field u is given by u = α 1 u 1 + α 2 u 2 , while the surface density of hydrodynamic friction forces f is f = α 1 f 1 + α 2 f 2 , where u i and f i denote the flow field and the surface density of hydrodynamic friction forces corresponding to boundary condition v i , respectively, for i = 1, 2.
Hydrodynamic dissipation. We introduce the rate of work R (h) exerted by the surface on the fluid For incompressible Newtonian fluids at zero Reynolds number, 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 Φ = 2µ ∆ : ∆, where ∆ : ∆ = i,j ∆ ij ∆ ij denotes tensor contraction. The dissipation rate can be rewritten as Φ = ∇ · (u · σ) using Eqs. (1), (2) and the incompressibility condition ∇ · u = 0. Gauss divergence theorem now gives [27] (using Here, V denotes the three-dimensional fluid domain with boundary surface S. At finite Reynolds numbers, 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.

LAGRANGIAN MECHANICS: GENERALIZED COORDINATES
We consider a shape-changing surface S(t). While a description of all possible shape changes of S would require an infinite number of degrees of freedom, in important applications, we can restrict ourselves to a constrained set of shape changes characterized by a small number of shape coefficients, or generalized coordinates, q 1 , . . . , 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 ϕ 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 ϕ 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, q = (q 1 , . . . , q n ). The shape dynamics of the active surface S(t) = S[q(t)] is thus entirely described by the dynamics of q(t). In particular, the local rate of surface displacement depends linearly on the generalized velocitiesq i as where the normalized velocity fields w i (x) = ∂x/∂q i depend on q(t) but not onq(t). In fact, Eq. (10) simply generalizes Eq. (7) to the case of generalized coefficients α i =q i with units of a generalized velocity. Correspondingly, the surface distribution of hydrodynamic friction forces f (x) is given as a linear combination where the normalized force densities g i (x; q) = f i (x)/α i correspond to the surface density of hydrodynamic friction forces f i (x) induced by the velocity field v i (x) = α i w i (x; 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 S in the set of generalized coordinates. Thereby, the self-propulsion of shape-changing microswimmers can be described using the same formalism, see the section of rigid body transformations below. Undulating sheet with two wave modes. The amplitudes q1, q2 of the wave modes represent generalized coordinates of the shape-changing surface 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: qi for translations parallel to the ei-axis, and qi+3 for rotations around the ei-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, q6+1 and q6+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 q6+1 and q6+2 can be treated as two generalized coordinates. (E) Similarly, Dreyfus' rotator consists of three segments connected at a single joint; the relative angles q6+1 and q6+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 q6+1 with periodicity of 2π. (G) Prototypical flagellar beat pattern of a sperm cell, parametrized by a 2π-periodic phase variable, which defines a generalized coordinate q6+1. For the amplitude of regular flagellar bending waves and mean flagellar curvature, we used parameters from [30].

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] The superposition principle for the shape changes w i (x), Eq. (10), allows us to rewrite the total hydrodynamic dissipation rate R (h) as a sum of products of generalized velocities times their conjugate generalized friction force 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 F or total torque 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 velocitiesq i by introducing generalized hydrodynamic friction coefficients The generalized friction coefficients can be computed as scalar products between the (normalized) velocity profiles w i (x), and the (normalized) force profiles Alternatively, we could express Γ ij in terms of partial derivatives with respect to the generalized velocitiesq i as . We refer to diagonal elements Γ ii of the generalized hydrodynamic friction matrix Γ as self-friction coefficients. Off-diagonal elements Γ ij , i = 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 velocityq The hydrodynamic dissipation rate 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 = ∂R (h) /∂q i . (Following standard notation, the Rayleigh dissipation function is actually R (h) /2 [42]). The matrix Γ is symmetric, which represents a special case of Onsager reciprocity [27,56] Γ ji (q) = Γ ij (q) The matrix is also positive semi-definite, consistent with the fact that the rate of energy dissipation should be nonnegative. (In fact, Γ should be positive definite, except maybe at singular points q in configuration space, where w i = ∂x/∂q i , i = 1, . . . , n are linearly dependent.) In addition to hydrodynamic dissipation as characterized by 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].

EQUATION OF MOTION
Balance of generalized forces. We introduce active driving forces Q i , i = 1, . . . , 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 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 , . . . , q n , and in the choice of the active driving forces Q 1 , . . . , 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 velocitiesqq where Q = (Q 1 , . . . , 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 = i, i.e., we may assume Q i = Q i (ϕ 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 the application section below, 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.

RIGID BODY MOTION OF A SELF-PROPELLED MICROSWIMMER
The above formalism includes the important application case of shape-changing microswimmers and their selfpropulsion 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 S and introduce a material frame of this microswimmer consisting of a reference point x 0 and a set of orthonormal vectors e 1 , e 2 , e 3 .
A rigid body motion of the swimmer is characterized by a translation of its reference point,ẋ 0 = v 0 = v 1 e 1 + v 2 e 2 + v 3 e 3 , and a rotation of its material frame withė i = ε ijk Ω j e k , where ε ijk denotes the Levi-Cevita symbol and we use Einstein summation convention. The components v 1 , v 2 , v 3 and Ω 1 , Ω 2 and Ω 3 of the translational and the rotational velocity vector with respect to the basis e 1 , e 2 , e 3 , respectively, represent the six degrees of freedom of rigid body motion and satisfy We choose these velocity components as the six generalized velocitieṡ The coordinates q 1 , . . . , q 6 are elements of the Lie group se(3) = R 3 × 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 F and the total hydrodynamic friction torque M, respectively P 1 = F · e 1 , P 2 = F · e 2 , P 3 = F · e 3 , P 4 = M · e 1 , P 5 = M · e 2 , P 6 = M · e 3 .
In this case, the 6 × 6-matrix of generalized hydrodynamic friction coefficients Γ 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 x ∈ S on the surface 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 r 0 of the swimmer is mapped to the origin 0 ∈ 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 x of a point x ∈ S under this transformation, i.e., the coordinates of the surface after it has been brought into a reference condition [62]. Eq.
The superposition principle of low-Reynolds number flow, Eq. (11), implies that the surface density f (x) of hydrodynamic friction forces can be written as a superposition of contributions due to rigid body motion and a contribution f act (x) due the active shape change where f act (x) depends only on the shape changeẋ , but not on the translational velocity v 0 nor the rotational velocity Ω.
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, F = F ext , M = M 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 Eq. (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, f act (x), gives rise to a contribution F act = S d 2 x f act (x) to the total force, as well as an analogous contribution M 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, v 0 and Ω.
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 F ext , or external torques M ext are captured by corresponding external forces Q ext i . Examples include gravitational force for a non-buoyant swimmer or torques exerted by an external rotating magnetic fields on an artificial microswimmer with non-zero 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 R (h) , where P iqi equals the rate of work required for the shape change with rateq 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 , . . . , q 3 representing rigid body translation, enforcingq 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 M = P 4 e 1 + P 5 e 2 + P 6 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 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,q i = −u 0 · e i , i = 1, 2, 3. Correspondingly, the total hydrodynamic friction force F = P 1 e 1 + P 2 e 2 + P 3 e 3 represents the constraining force required to clamp the microswimmer in such an external flow.

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 Γ for a set of reference configurations q and save this as a look-up table; the friction matrix Γ(q) for arbitrary 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 S. Assume the surface S is represented by a triangulated mesh. The triangular faces (or 'elements') shall be enumerated by k ∈ I with midpoints x k and respective areas A k .
In a first step, we compute a (normalized) surface distribution of velocities w i (x k ), k ∈ I for each generalized coordinate i = 1, . . . , n, either by computing the derivative w i (x k ) = ∂x k (q)/∂q i analytically, or by evaluating the finite difference quotient for each midpoint x k , k ∈ I, where ∆ i is the unit vector whose components are all zero, except the i th component, and ε is a small number. We can use boundary element methods to numerically compute a surface density of hydrodynamic friction forces f (x k ) with physical units of a stress, given an arbitrary surface distribution of velocities v(x k ) specified at each midpoint x k , k ∈ 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 f j (x k ) = α j g j (x k ) of hydrodynamic friction forces, corresponding to the velocity distribution v j (x k ) = α j w j (x k ). Here, α j is an arbitrary constant to ensure proper physical units of a velocity for v j . We thus obtain n surface distributions of (normalized) hydrodynamic friction forces g j (x k ), j = 1, . . . , n, one for each generalized coordinate q j . These force distributions g j (x k ) depend on q, but not onq. Finally, we compute the components Γ ij of the generalized hydrodynamic friction matrix Γ by taking the scalar product of the i th (normalized) velocity distribution w i (x k ), and the j th (normalized) force distribution g j (x k ) where A k was the area of the k th triangle. We can interpret α j g j (x k )A k at the total force acting on the k th element (with proper physical units of a force) if the generalized coordinate q i would change at a rate α i . Importantly, it suffices to compute the generalized hydrodynamic friction matrix Γ 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, . . . , n for m n different choices of q. By exploiting symmetries, as well as translational and rotational invariance for individual microswimmers, this number can be reduced further. The friction matrix Γ(q) for arbitrary q can then be found by interpolation. For example, 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 f (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 relieve 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].

APPLICATION: PAIR OF INTERACTING CILIA
We demonstrate our LAMAS modeling framework using the example of hydrodynamic synchronization in pairs of cilia. We thereby 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 parametrized by a 2π-periodic phase variable ϕ, see Fig. 3A. For unperturbed beating, the phase speed equals the angular frequency of the cilia beat,φ(t) = ω 0 .  [12] as reported in [68], shown as sequence of three-dimensional shapes parameterized by a 2π-periodic phase variable ϕ (color code). Spacing of square grid: 2 µm. (B) We consider a pair of cilia with respective phases ϕ1 and ϕ2, whose base points are separated by a distance d along a direction that encloses an angle ψ with the x axis (where the y axis is set by the direction of the effective stroke of both cilia). (C) Self-friction coefficient Γ11(ϕ1) of a single cilium as function of its phase variable ϕ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, Γ11 is proportional to the phase-dependent active cilia driving force Q1(ϕ1). (D) Generalized hydrodynamic friction coefficient Γ12(ϕ1, ϕ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 µm, orientation angle ψ = π/3. (E) The magnitude of hydrodynamic interactions, here quantified by the L2-norm of Γ12, decay as ∼ 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 (ψ = 0: dark-blue, ψ = π/3: teal, ψ = 2π/3: light-green; also indicated by the direction arrows.) (F) We characterize the stability of the in-phase synchronized state, defined by ϕ1(t) = ϕ2(t), for different relative orientations of the two cilia by a Lyapunov exponent λ, see Eq. (30). Colored dots at respective positions in the xy-plane represent the value of λ 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, λ < 0), while positive values imply that in-phase synchronization is linearly unstable (red colors, λ > 0). (G) We determined the steady-state phase difference δ * between the two cilia for different relative cilia positions, analogous to panel F. While δ * = 0 for cilia orientations with stable in-phase synchronization (cyan), we observe anti-phase synchronization with δ * ≈ π for cilia orientations with λ > 0 (red colors). For relative cilia orientation aligned with the direction of the effective stroke of the cilia beat (ψ = π/2), we observed cases of multi-stability (bi-colored dots). (H) Consistent with the far-field scaling of hydrodynamic interactions as shown in panel E, we find that also the Lyapunov exponent λ, 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: ω0/(2π) = 32 Hz [12], fluid viscosity: µ = 10 −3 Pa s.

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 ϕ 1 and ϕ 2 fully characterize the dynamics of the two beating cilia, and represent a set of generalized coordinates with state vector q = (ϕ 1 , ϕ 2 ).
For our example, the force balance equation, Eq. (18) takes the form This equation can be further simplified. The symmetry relation Eq. (17) implies Γ 12 (ϕ 1 , ϕ 2 ) = Γ 21 (ϕ 1 , ϕ 2 ). Numerical computation of Γ 11 (ϕ 1 , ϕ 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), Γ 11 (ϕ 1 , ϕ 2 ) and Γ 22 (ϕ 1 , ϕ 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φ 1 = ω 0 , where ω 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 (ϕ 1 ). Specifically, for a single cilium, the force balance equation reads, Q 1 (ϕ 1 ) = Γ 11 (ϕ 1 )φ 1 . We conclude Q 1 (ϕ 1 ) = ω 0 Γ 11 (ϕ 1 ); Fig. 3C displays the phase-dependence of Γ 11 (ϕ 1 ). Since both cilia are assumed identical with same intrinsic beat frequency ω 0 , this also specifies the active driving force Q 2 (ϕ 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φ Eq. (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 Γ ij (ϕ 1 , ϕ 2 ) can be found in the appendix. An example of the generalized friction coefficient Γ 12 (ϕ 1 , ϕ 2 ), which characterizes hydrodynamic interactions between the two cilia, is shown in Fig. 3D; Fig. S2 shows Γ 12 (ϕ 1 , ϕ 2 ) for additional cilia orientations.
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 Γ 12 (ϕ 1 , ϕ 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].
We define a dimensionless Lyapunov exponent as which characterizes whether the initial perturbation decays or grows. The in-phase synchronized state is linearly stable if |δ 1 | < |δ 0 | (hence λ < 0), and linearly unstable if |δ 1 | > |δ 0 | (hence λ > 0). Fig. 3F shows λ 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 λ > 0 for direction angles ψ = π/6 and π/3, 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 δ * that correspond to stable periodic solutions, see Fig. 3G. As a technical point, δ(t) may weakly oscillate during each cycle; we therefore define δ * as the phase difference at times for which ϕ is an integer multiple of 2π.
When the in-phase synchronized state is linearly stable for a given cilia configuration, we obviously have δ * = 0. If, however, the in-phase synchronized state is linearly unstable, we approximately find δ * ≈ π, corresponding to anti-phase synchronization. For a few cilia configurations, we observe multi-stability, characterized by two different values of the phase differences δ * that correspond to stable periodic solutions; these configurations are indicated as bi-colored half circles in Fig. 3G.
The magnitude |λ| 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.

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 R (h) , which defines generalized friction forces P i conjugate to each generalized coordinate q i as a partial derivative 2P i = ∂R (h) /∂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 to efficiently explore biological fluid dynamics problems at low Reynolds numbers 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, which previously had been either studied using detailed simulations with many degrees of freedom [70,71], or using minimal models [69,[72][73][74]. A key simplifying assumption will be to approximate manybody 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,75]. This active noise jeopardizes synchronization of cilia by hydrodynamic interactions. Additionally, noise randomizes the motion of biological microswimmers. While thermal noise causes noticeable rotational diffusional of micrometer-sized bacteria such as E. coli [76], amplitude fluctuations of cilia beating affect the swimming of ten-fold 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-ranged 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 visco-elastic 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. FIG. S2: Hydrodynamic interaction as function of cilia phases for different cilia orientations. Generalized hydrodynamic friction coefficient Γ12(ϕ1, ϕ2) as in Fig. 3D for different cilia orientation angles: left: ψ = π/2 (direction of effective stroke), middle: ψ = π/3, right: ψ = 0 (perpendicular to direction of effective stroke). Cilia distance: d = 18 µm. Note different color scale compared to Fig. 3D.