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

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-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-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–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-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-30], the more refined method of slender-body theory, which considers a line distribution of hydrodynamic singularities (point forces) along a filament [31-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-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-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–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.

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-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 . Here, denotes the flow velocity, 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 in the bulk of the fluid. Special conditions apply at boundaries.

No-slip boundary condition for an active surface We consider a surface immersed in the fluid that changes its shape as a function of time. For example, 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 for each point at time t.

We impose a no-slip boundary condition at this surface, i.e., require that the local velocity of fluid flow matches the local velocity of the surface for each surface pointHydrodynamic friction forces A shape change of the surface induces a flow field with corresponding stress tensor field . The stress determines the surface density of forces exerted by the surface on the fluid (with units of a stress , also called contact force, or traction force density)where is the surface normal pointing into the fluid. Correspondingly, 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 Analogously, the total torque (with respect to a reference point ) exerted by the surface on the fluid is given bySuperposition 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 that is given as a linear combination of velocity distributions and aswith real coefficients . Then, the resultant flow field is given by , while the surface density of hydrodynamic friction forces is , where and denote the flow field and the surface density of hydrodynamic friction forces corresponding to boundary condition , respectively, for .

Hydrodynamic dissipation We introduce the rate of work exerted by the surface on the fluidFor incompressible Newtonian fluids at zero Reynolds number, equals the instantaneous rate of hydrodynamic energy dissipation within the fluid [27]. Indeed, let us consider the local dissipation rate, which is given by , where denotes tensor contraction. The dissipation rate can be rewritten as using Eqs. (1), (2) and the incompressibility condition . Gauss divergence theorem now gives [27] (using for )Here, V denotes the three-dimensional fluid domain with boundary surface . At finite Reynolds numbers, 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 . While a description of all possible shape changes of 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, .

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–55]. Elastic degrees of freedom arising from waveform compliance can be incorporated in such a framework as additional amplitude degrees of freedom [10, 48].

Fig. 4

Generalized coordinates: Examples. a Undulating sheet with two wave modes. The amplitudes , of the wave modes represent generalized coordinates of the shape-changing surface . 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: for translations parallel to the -axis, and for rotations around the -axis, , 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, and , 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 and can be treated as two generalized coordinates. e Similarly, Dreyfus’ rotator consists of three segments connected at a single joint; the relative angles and 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 with periodicity of . g Prototypical flagellar beat pattern of a sperm cell, parameterized by a -periodic phase variable, which defines a generalized coordinate . For the amplitude of regular flagellar bending waves and mean flagellar curvature, we used parameters from [30]

We introduce the state vector, . The shape dynamics of the active surface is thus entirely described by the dynamics of . In particular, the local rate of surface displacement depends linearly on the generalized velocities aswhere the normalized velocity fields depend on but not on . In fact, Eq. (10) simply generalizes Eq. (7) to the case of generalized coefficients with units of a generalized velocity. Correspondingly, the surface distribution of hydrodynamic friction forces is given as a linear combinationwhere the normalized force densities correspond to the surface density of hydrodynamic friction forces induced by the velocity field . An example of a surface velocity field with corresponding surface density of hydrodynamic friction forces is shown in Fig. 2.

Fig. 5

Hydrodynamic friction forces: Example of cilium during power stroke. a Surface velocity on a shape-changing surface , here given by a slender cilium (blue) attached to a no-slip boundary surface (gray); the cilium progresses with phase velocity along its periodic beat cycle. For visualization, the three-dimensional shape of the cilium as well as was projected on the yz-plane (see Fig. 3a for a three-dimensional representation). b Corresponding surface distribution of hydrodynamic friction forces 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, and , we can compute a generalized hydrodynamic friction coefficient, , which is proportional to the phase-dependent rate of energy dissipation in the surrounding fluid. In the general case of n generalized coordinates , we obtain a matrix , see also Eq. (15). c Flow field induced by the active shape change of the cilium, here shown as two-dimensional section at . The color represents the magnitude of three-dimensional velocity vectors, whereas white arrows represent the projections of on the yz-plane. The flow field was computed as convolution of the fundamental solution of the Stokes equation with the force distribution . Cilium phase corresponding to Fig. 3a: , cilium beat frequency: [12], dynamic viscosity of fluid: (corresponding to viscosity of water at )

The formalism allows to include also rigid body transformation such as translations and rotations of the surface 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. 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 -periodic phase variable (color code). Spacing of square grid: . b We consider a pair of cilia with respective phases and , 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 of a single cilium as function of its phase variable , 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, is proportional to the phase-dependent active cilia driving force . d Generalized hydrodynamic friction coefficient 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 , orientation angle . e The magnitude of hydrodynamic interactions, here quantified by the -norm of , decay as , consistent with the theoretical scaling expected from the Blake tensor [65]. Different curves correspond to different separation directions between the two cilia (: dark blue, : light green, : teal; also indicated by the direction arrows.) f We characterize the stability of the in-phase synchronized state, defined by , 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, ), while positive values imply that in-phase synchronization is linearly unstable (red colors, ). g We determined the steady-state phase difference between the two cilia for different relative cilia positions, analogous to panel F. While for cilia orientations with stable in-phase synchronization (cyan), we observe anti-phase synchronization with for cilia orientations with (red colors). For relative cilia orientation aligned with the direction of the effective stroke of the cilia beat (), 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 , which represents an effective synchronization strength, likewise decays as with distance d between the two cilia. Different curves represent different separation directions, analogous to panel E. Frequency of cilia beat: [12], fluid viscosity:

Generalized hydrodynamic friction forces

We introduce generalized hydrodynamic friction forces conjugate to the generalized coordinates , following the convention of Lagrangian dynamics of dissipative systems [42], see also [43, 46]The superposition principle for the shape changes , Eq. (10), allows us to rewrite the total hydrodynamic dissipation rate as a sum of products of generalized velocities times their conjugate generalized friction forceNote that the different generalized coordinates 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 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 or total torque 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 by introducing generalized hydrodynamic friction coefficientsThe generalized friction coefficients can be computed as scalar products between the (normalized) velocity profiles , and the (normalized) force profiles , see also Fig. 2Alternatively, we could express in terms of partial derivatives with respect to the generalized velocities as . We refer to diagonal elements of the generalized hydrodynamic friction matrix as self-friction coefficients. Off-diagonal elements , , 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 The hydrodynamic dissipation rate 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 . (Following standard notation, the Rayleigh dissipation function is actually [42].)

The matrix is symmetric, which represents a special case of Onsager reciprocity [56]. The proof follows directly from the Lorentz reciprocal theorem [27]: Let , and , denote the flow field and stress tensor associated with a change of only with rate , or a change of only with rate , respectively, while all other generalized coordinates are kept constant; then,where we used the Lorentz reciprocal theorem at .

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 in configuration space, where , are linearly dependent.)

In addition to hydrodynamic dissipation as characterized by , 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 , 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–59]. Our active driving forces generalize the active driving forces considered in these models.

We postulate a balance of generalized forces between driving forces and hydrodynamic friction forcesWe 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 and in the choice of the active driving forces .

From the force balance equation, Eq. (18), and Eq. (12) expressing the generalized forces , we obtain equations of motion for the generalized velocities where is the vector of active forces.

Calibration of active driving forces Because each driving force characterizes internal processes, it is plausible to assume that only depends on the corresponding degree of freedom , but not on the other , , i.e., we may assume . This assumption will hold in particular in applications, where the index i enumerates different microswimmers or different cilia. In principle, may additionally depend on the friction force 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 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 coarse-grain internal active processes.

Additionally, the formalism allows to incorporate internal elastic degrees of freedom and the corresponding elastic restoring forces 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 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 and introduce a material frame of this microswimmer consisting of a reference point and a set of orthonormal vectors , , .

A rigid body motion of the swimmer is characterized by a translation of its reference point, , and a rotation of its material frame with , where denotes the Levi-Cevita symbol and we use Einstein summation convention. The components , , and , and of the translational and the rotational velocity vector with respect to the basis , , , respectively, represent the six degrees of freedom of rigid body motion and satisfy , , , and , , .

We choose these velocity components as the six generalized velocitiesFormally, the coordinates are elements of the Lie group 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 and the total hydrodynamic friction torque , respectivelyIn this case, the 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 , , relative to the swimmer’s material frame for each point on the surfaceWe introduce the time-dependent rigid body transformation that maps the material frame of the swimmer to the laboratory frame, such that the reference point of the swimmer is mapped to the origin , and the material frame vectors are mapped to the standard unit vectors, respectively. The coordinates , , are then just the coordinates of the image of a point 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 of the surface into a contribution stemming from the rigid body motion and a contribution stemming only from any shape changeThe superposition principle of low-Reynolds number flow, Eq. (11), implies that the surface density of hydrodynamic friction forces can be written as a superposition of contributions due to rigid body motion and a contribution due the active shape changewhere depends only on the shape change , but not on the translational velocity 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, , [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 itselfEquation (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, , gives rise to a contribution to the total force, as well as an analogous contribution 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, 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 or external torques are captured by corresponding external forces . 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 with prescribed driving protocol , one would omit the corresponding force balance equation from the set of equations Eq. (18), and solve for the equation of motion of the other coordinates with prescribed . The conjugate hydrodynamic friction force nonetheless appears in the formula for the total hydrodynamic dissipation rate , where equals the rate of work required for the shape change with rate . A number of classical theory publications on self-propelled biological microswimmers considered prescribed shape dynamics [28, 49–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 with prescribed dynamics for the coordinates representing rigid body translation, enforcing . The conjugate hydrodynamic friction force equals the external constraining force required to impose the constraint. Similarly, to constrain a microswimmer from rotational motion requires a constraining torque . 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 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, , . Correspondingly, the total hydrodynamic friction force 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 and save this as a look-up table; the friction matrix for arbitrary 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 . Assume the surface is represented by a triangulated mesh. The triangular faces (or “elements”) shall be enumerated by with midpoints and respective areas .

In a first step, we compute a (normalized) surface distribution of velocities , for each generalized coordinate , either by computing the derivative analytically or by evaluating the finite difference quotientfor each midpoint , , where is the unit vector whose components are all zero, except the ith component, and is a small number.

We can use boundary element methods to numerically compute a surface density of hydrodynamic friction forces with physical units of a stress, given an arbitrary surface distribution of velocities specified at each midpoint , . 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 of hydrodynamic friction forces, corresponding to the velocity distribution . Here, is an arbitrary constant to ensure proper physical units of a velocity for . We thus obtain n surface distributions of (normalized) hydrodynamic friction forces , , one for each generalized coordinate . These force distributions depend on , but not on . Finally, we compute the components of the generalized hydrodynamic friction matrix by taking the scalar product of the ith (normalized) velocity distribution , and the jth (normalized) force distribution where was the area of the kth triangle. We can interpret at the total force acting on the kth element (with proper physical units of a force) if the generalized coordinate would change at a rate .

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 times, as we need to change each of the , for different choices of . By exploiting symmetries, as well as translational and rotational invariance for individual microswimmers, this number can be reduced further. The friction matrix for arbitrary 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 . 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].

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 -periodic phase variable , see Fig. 3a. For unperturbed beating, the phase speed equals the angular frequency of the cilia beat, .

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 and fully characterize the dynamics of the two beating cilia, and represent a set of generalized coordinates with state vector .

For our example, the force balance equation, Eq. (18), takes the formThis equation can be further simplified. The symmetry relation Eq. (17) implies for the hydrodynamic interactions. Numerical computation of 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), and , 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 , where denotes the intrinsic beat frequency of the cilium if there are no interactions with other cilia. This requirement uniquely determines the active driving force . Specifically, for a single cilium, the force balance equation reads, . We conclude ; Fig. 3c displays the phase dependence of . Since both cilia are assumed identical with same intrinsic beat frequency , this also specifies the active driving force of the second cilium.

Equation of motion Using the force balance equation, Eq. 28, and the calibrated driving force, we obtain the equation of motionEquation (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 can be found in the appendix. An example of the generalized friction coefficient , which characterizes hydrodynamic interactions between the two cilia, is shown in Fig. 3d; Fig. 5 shows for additional cilia orientations.

Fig. 8

Hydrodynamic interaction as function of cilia phases for different cilia orientations. Generalized hydrodynamic friction coefficient as in Fig. 3d for different cilia orientation angles: left: (direction of effective stroke), middle: (oblique to direction of effective stroke), right: (perpendicular to direction of effective stroke). Cilia distance: . Note the different color scale compared to Fig. 3d

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

Hydrodynamic interactions decay as . For large separation distances d between the two cilia, hydrodynamic interactions between the two cilia as characterized by decay as , 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 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 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 during one beat cycle. Specifically, we integrated Eq. (29) with the initial condition and for a small perturbation up to time T defined by (corresponding to the completion of a full beat cycle), and recorded the new phase difference .

We define a dimensionless Lyapunov exponent aswhich characterizes whether the initial perturbation decays or grows. The in-phase synchronized state is linearly stable if (hence ) and linearly unstable if (hence ).

Figure 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 for direction angles and , 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, may weakly oscillate during each cycle; we therefore define as the phase difference at times for which is an integer multiple of .

When the in-phase synchronized state is linearly stable for a given cilia configuration, we obviously have . 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 bicolored half circles in Fig. 3g.

The magnitude of the Lyapunov exponents decreases as 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.

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 , which defines generalized friction forces conjugate to each generalized coordinate as a partial derivative . 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–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.