## Abstract

As a model of ciliary beat, we use two-state oscillators that have a defined direction of oscillation and have strong synchronization properties. By allowing the direction of oscillation to vary according to the interaction with the fluid, with a timescale longer than the timescale of synchronization, we show in simulations that several oscillators can align in a direction set by the geometrical configuration of the system. In this system, the alignment depends on the state of synchronization of the system, and is therefore linked to the beat pattern of the model cilia. By testing various configurations from two to 64 oscillators, we deduce empirically that, when the synchronization state of neighbouring oscillators is in phase, the angles of the oscillators align in a configuration of high hydrodynamic coupling. In arrays of oscillators that break the planar symmetry, a global direction of alignment emerges reflecting this polarity. In symmetric configurations, where several directions are geometrically equivalent, the array still displays strong internal cooperative behaviour. It also appears that the shape of the array is more important than the lattice type and orientation in determining the preferred direction.

## 1. Introduction

Motile cilia are highly conserved organelles, present in a wide range of eukaryotic organisms [1,2]. The periodic beating of these filaments allows very diverse biological processes, from swimming in unicellular organisms such as the algae *Chlamydomonas reinhardtii* to fluid circulation in the human brain ventricles [3], fallopian tubes [4] and airways [5]. Cilia grow outwards from a structure called the basal body, which is anchored to a cell's cytoskeleton [6]. The basal body is itself generated from centrioles, which are subcellular structures responsible for microtubular organization. Once fully grown, a cilium transverses the cell's plasma membrane extending typically several micrometres out of the cell body. The interior structure of motile cilia, the axoneme, is in most cases composed of nine exterior microtubular doublets, and two central microtubules, connected by a well-defined network of motor proteins that slide the microtubules, leading to bending of the entire structure; the mechanisms at work within cilia have been modelled quantitatively [7] and models to reproduce the cycle of oscillation have been proposed [8,9]. In the earlier-mentioned mammalian examples, the cilia belong to multi-ciliated cells, and there are typically of order 200 cilia per cell; each cilium is separated by approximately 200 nm. The same is true and has been explored in the outer surface of various ‘model’ biological organisms, such as *Paramecium* [10] and the algae colony *Volvox* [11] as well as developing embryos of mice [12,13] and *Xenopus* frogs [14]. These ciliated cells, together with mucus-producing cells (goblet cells, in mammals), and ion-regulating cells, form a general tissue type known as mucociliary epithelium.

In fully developed epithelial tissue, the cilia undergo a periodic sequence of power and recovery strokes, with the filament roughly contained in a plane perpendicular to the plane of the cells. The direction of beating is well-defined relative to the organ, for example, it is parallel to the trachea in mammalian airways. This is essential for the generation of fluid flow, mucus clearance away from the lungs in this case, which relies on coordinated beating of cilia to produce transport-efficient metachronal waves [15]. An open question in developmental biology is to find the rules and the cues that enable this fairly complex, and well organized tissue, to be made. The first symmetry to be broken in the development of vertebrates is the anterior/posterior. For example, the planar cell polarity (PCP) pathway sets the initial direction in a developmental stage of the epithelium in *Xenopus* embryos, a tissue that includes multi-ciliate cells [16]. From this point on, there are gradients of a variety of biochemical elements along this axis. Cells can be polarized, both in the intracellular protein localization and in their shape. This process happens before ciliagenesis, and the standard view in biology is that the gradients in biochemical markers control most, if not all, of the subsequent organ development. However, once cilia are generated, they contribute to long-range flows, which can transport chemical factors directionally, or act as a mechanical cue for organization [17–20]. Recent results and reviews have been summarized here, but these questions have been in the minds of researchers for some time [21].

In the specific case of developing orientational cilia order in the airway tissue, there is a hypothesis that flow-induced self-organization might be important. A fraction of the cells in the tissue that will develop into the airway epithelium express a few hundred centrioles, which become basal bodies and grow cilia. At this stage, the cells themselves are already polarized (biochemically and in shape), but the basal bodies when they first appear are not oriented. There have been very recent studies suggesting that the network of microtubules connecting the basal bodies could couple to the cell shape or to the emerging tissue architecture, and possibly orient the cilia [22]. On the other hand, the newly made cilia are exposed to a directional flow from the mucus being produced by other cells. These cilia will also be exposed to the flow that they themselves generate, i.e. they interact with each other through hydrodynamic interaction forces [23]. The question of how are cilia aligned has also been addressed looking at mouse brain ventricles [24], where it was shown that cilia first dock apically with random orientations, and then reorient in a common direction through a coupling between hydrodynamic forces and the PCP protein Vangl2.

This paper is an attempt to explore whether the collective interaction through flow can feasibly explain the very robust alignment of basal bodies with each other inside each cell, and also the alignment to the axis of the developing tissue. It has already been proposed that hydrodynamic interaction might play a central role in the emergence of a collective *synchronized* state of cilia (metachronism) [23,25,26]. Evidence of the role of flows in determining the orientation of cilia is present in the experiments in the *Xenopus* larval skin [17,27] and in mouse brain ventricles [24], but there is a lack of physical models to explain this behaviour. Moreover, these experiments suggest that the positive feedback of the fluid on the orientation of the cilia is due to an average fluid flow. We present here an alternative process to explain the emergence of polar order by hydrodynamic coupling, supported by simulations and a theoretical argument. This model does not require an average net flow. Instead, it involves the time-dependence of the flow generated locally by each cilium, and is therefore related to the relative phase differences between the cilia and hence to the shape of the beating cycle of a single cilium.

In the spirit of a reductionist model, a great number of factors are either neglected or coarse-grained. This gives us a simple (albeit intrinsically nonlinear) rule that describes each cilia element in the model. All the degrees of freedom (the details) that characterize the structure, the semi-flexible nature, and the molecular motor drive of the cilia, are coarse-grained into the rule by which the cilia is moved in the model. The cilium shape is approximated by a sphere for the sake of calculating the induced (and perceived) forces owing to fluid flow; the presence of a solid planar boundary is neglected here, and the liquid is described by a constant viscosity (Newtonian liquid). All these factors can be developed further in future work, and, in fact, to match quantitatively with a given biological scenario, it will be essential at least to adjust the parameter range (e.g. the distance between units) to more specific systems and to add realistic fluid dynamics. The mucus layer in the airway system is spatially stratified [15] and exhibits a viscoelastic response. Despite all these simplifications, the model remains very rich: the basic dynamical rule for each cilium, together with the many-body fluid flow interactions that couple all the active elements with each other, lead to synchronization in a variety of collective dynamical behaviour. In the ‘simple’ system, it is possible to test the consequences on the emergent collective state coming from subtle changes in the forces driving cilia flow, and in the geometrical placement of the active elements investigating the role of local structure and boundary conditions.

Other simple models have been considered in recent years. Synthetic experimental ciliary systems have been assembled based on filaments of magnetic colloidal particles tethered to a surface, and driven by magnetic field gradients [28,29], and these might be a way to drive fluid flow of mixing in microfluidic channels and cavities. These controlled systems are also well suited to exploring the hydrodynamics at the local scale [30]. In other models, the question of synchronization through flow interactions has been explored [31,32]. It is also worth keeping in mind that schools of swimming organisms can certainly align to an external field above a density threshold; self-ordering is also known in swimmers, but its origin is under debate, because there are multiple sources of interaction, including steric and hydrodynamic coupling [33].

Fluid flow at low Reynolds' number (*Re*) is reviewed in [34]. In this regime of lengthscales and velocities, there is no inertia, and an object's velocity is directly proportional to the instantaneous net force acting on it. The flow field around a spherical object decays as the inverse distance from the object, at far field and in the absence of boundaries. As a model for motile cilia, we develop a simple phase oscillator that we studied in previous work to determine the threshold for hydrodynamic synchronization [35,36] and the properties of the dynamical collective state [37–40]. In this model, a configuration-dependent force (‘geometrical switch’) is applied to a colloidal sphere, generating an oscillation of fixed amplitude but free phase. In this work, we develop this model, by adding freedom to orient the direction of the oscillations, in order to look at possible alignment of several oscillators. This work differs from the theoretical model studied in [41], where it was shown that alignment can occur in the direction of an average fluid flow that can be generated by the cilia themselves. Here, we show that an external flow is not a necessary condition for orientational alignment. Alignment can emerge spontaneously for certain spatial arrangements of oscillators. We note here that cells in some systems as, for example, *Paramecium* show very regularly spaced cilia, possibly defining axes of alignment, whereas in the airways, the cilia distribution is more random within cells, and the distribution of the multi-ciliated cells themselves is also not completely regular.

## 2. Model of active oscillators with orientational freedom

The motion of a single cilium is modelled by a spherical colloidal particle driven by an external potential *U* = *k*(*x* − *x*_{0})* ^{α}*. In order to actively drive the particle, the centre of the trap

*x*

_{0}is switched alternately between two positions, represented by the grey lines in figure 1

*a*. When the trap

*B*is turned on, the bead is driven by the right (blue) potential. When it reaches a distance

*ξ*from the trap's centre, trap

*B*is switched off, and trap

*A*is turned on, driving the bead in the opposite direction. The same switching rule applies to the new trap. This is a realization of a one-dimensional free-phase oscillator with a constant amplitude 2

*R*

_{g}. We have implemented this oscillator experimentally with optical tweezers in previous work by using harmonic driving potentials [35,38,42] and also with variable potential shapes, described by a variable power exponent

*α*[40]. It can be assumed here that two such oscillators display strong synchronization—through hydrodynamic coupling—provided that

*α*≠ 1; the emerging dynamical steady state of a pair is in phase if

*α*< 1 and in antiphase if α > 1 [40,43].

In order to add freedom in the direction of oscillation of the above-described oscillator, we allow the particle and the traps to move in the direction orthogonal to the direction defined by the centre of trap *A* and trap *B* (figure 1*b*): in the two-dimensional plane, the bead can move in a disc of radius *R _{g}* and the traps can move on a circle of radius

*R*+

_{g}*ξ*. The angular freedom is implemented such that the trap moves by following the angular position of the bead only if (i) the bead is outside the disc of radius

*R*and (ii) it is on the side of the active trap. When condition (i) is satisfied, the orthogonal trapping force is set to 0; otherwise, it is set to a high value such that the bead is strongly confined along the direction defined by the traps. The radius

_{f}*R*defining the boundary between angular freedom and confinement is a parameter that controls how much angular freedom is allowed. Condition (ii) ensures that an orthogonal displacement of the particle when it is approaching the active trap will not be cancelled by a possible orthogonal displacement in the opposite direction just after the trap switched. Figure 1

_{f}*b*shows the track of the particle position of a single oscillator over a few cycles of oscillation. The direction of oscillation is slightly changing because of the freedom of rotation when conditions (i) and (ii) are satisfied. Because the particle is not coupled, in this figure, to any other oscillator, the direction is only changing because of free diffusion of the particle. Figure 1

*c*represents the direction of the active trap. The square-like shape is due to the switching between traps

*A*and

*B*, and each jump corresponds to an angular change of

*π*. The angular freedom acts before a trap switch, as a slight change of the trap angle. Because this work focuses on the alignment of oscillators, all angles in later figures will be plotted modulo

*π*in order to remove the square-like shape in the graphs.

## 3. Alignment of two oscillators

In this section, we study the behaviour of two oscillators, set apart by a distance *d* = 10 µm. The simulations shown here use the same methods as in [40]. In particular, the hydrodynamic interaction between the spheres is described by the Oseen tensor, hence assuming a bulk fluid, and the simulations include thermal noise. Figure 2 shows a simulation for *α* = 0.7 and *R _{f}*/

*R*= 0.5. The other parameters are set to values characteristic of cilia and are given in figure 2. The relative state of the oscillators has two features. First, the oscillations synchronize in phase within a few cycles (figure 2

_{g}*b*; typically four cycles, with the parameters used in this work). The synchronized state agrees with [43] and [40]. Furthermore, the directions of the two oscillators converge to 0 rad, as shown in figure 2

*c*. The characteristic time to converge is, however, much higher than the time to synchronize: about 12 cycles for

*R*/

_{f}*R*= 0.5.

_{g}When varying *α* and *R _{f}*, different behaviours can be obtained. This is summarized in figure 3 showing the distribution of the angles (

*θ*

_{1},

*θ*

_{2}) of the active traps for each oscillator, when running 3000 s long simulations. For

*α*= 0.7, leading to in-phase oscillations, the angles converge to (0, 0) for all values of

*R*/

_{f}*R*except 0. When

_{g}*R*/

_{f}*R*= 0, thermal fluctuations tend to randomize the angle of the oscillators each time they cross their centre

_{g}*O*

_{1}and

*O*

_{2}, resulting in a uniform (

*θ*

_{1},

*θ*

_{2}) distribution. This also broadens the peak at (0, 0) when

*R*/

_{f}*R*≈ 0. The case of

_{g}*α*= 2, which leads to oscillations in antiphase, is more complex. For the same reason as before, the angle's distribution is uniform for

*R*/

_{f}*R*= 0. When this ratio is increased, the distribution shows peaks at (0, 0) (

_{g}*R*/

_{f}*R*= 0.05), a locus of angles avoiding the (0, 0) alignment (0.2 ≤

_{g}*R*/

_{f}*R*≤ 0.8) and again convergence to (0, 0) (

_{g}*R*/

_{f}*R*= 0.95). In the middle-range of angular freedoms 0.2 ≤

_{g}*R*/

_{f}*R*≤ 0.8, the two oscillators do not align. Instead, the system tends to stay in configurations that minimize the hydrodynamic coupling between the particles (e.g. the (0,

_{g}*π*/2) configuration).

## 4. An explanation for the synchronization of two oscillators

We derive here a model to explain the alignment properties observed in the simulations above. For the sake of simplicity, we neglect the Brownian fluctuations to reduce the question to a problem of convergence of the angles of the two oscillators. The distributions in figure 3 in the middle-range of *R _{f}*/

*R*can be explained for both

_{g}*α*> 1 and

*α*< 1 by a simple treatment in which the synchronization of the oscillators and their alignment are considered separately.

The state of an oscillator is described by the position *r* of the particle, the angle of the oscillator (in [0, *π*[) and a variable *σ* = ±1 indicating the direction along *θ* where the trap is. The first and third parameters can be merged into a single parameter, describing the ‘geometrical phase’, *ϕ* (in [0, 2*π*[). Therefore, the state of the oscillator is fully described by two parameters: *ϕ* and *θ*.

The evolution in time of *ϕ*_{1} and *θ*_{1} (left oscillator), *ϕ*_{2} and *θ*_{2} (right oscillator) is a complex problem in which the four variables are coupled. However, in the middle-range of *R _{f}*/

*R*, the phase difference

_{g}*ϕ*

_{2}−

*ϕ*

_{1}converges with a much lower relaxation time compared with the characteristic relaxation time of the angles

*θ*

_{1}and

*θ*

_{2}. When looking at the alignment properties of the oscillators, at timescales at which the directions move, we can assume that the synchronization of

*ϕ*

_{1}and

*ϕ*

_{2}occurs instantaneously: in-phase, if

*α*< 1 or in antiphase, if

*α*> 1, as described by the theory for the one-dimensional version of the oscillator [43]. The solving of the evolution of the system is reduced to the following question: how do (

*θ*

_{1},

*θ*

_{2}) evolve in time when the oscillations are assumed in phase (if

*α*< 1) or in antiphase (if

*α*> 1)?

Starting from an initial condition (*θ*_{1}, *θ*_{2}) of the traps' positions when the beads pass *O*_{1} and *O*_{2}, the angles after the traps switched for both beads become (*θ*_{1} + *Δ**θ*_{1}, *θ*_{2} + *Δ**θ*_{2}). It is convenient to introduce here three frames in figure 2*a*: , and . When the particle *i* is approaching its active trap, it will undergo the change of angle *Δ**θ*_{i} when the particle moves from the position *R _{f}* to

*R*.

_{g}*Δ*

*θ*

_{i}is related to the velocity of the particle along the orthogonal direction and to the time

*τ*it takes to move from

*R*to

_{f}*R*in the radial direction: 4.1where

_{g}*t*

_{f}_{,i}and

*t*

_{g}_{,i}are the instants at which the particle reaches a radial position

*R*and

_{f}*R*, respectively.

_{g}*σ*

_{i}is the variable describing which of the two traps that is driving the oscillator

*i*is active:

*σ*

_{i}= 1 if the active trap is in the direction of , or

*σ*

_{i}= −1 if it is in the opposite direction. Between these two positions, the orthogonal trapping force is zero. Neglecting thermal fluctuations, the only contribution to is from the hydrodynamic coupling. We describe the coupling with the Oseen tensor that relates the velocities

**v**

*of the particles to the driving forces*

_{i}**F**

*acting on them: 4.2where 4.3Here,*

_{j}*η*is the viscosity,

*a*the beads' radii and

**I**is the unit tensor. Equation (4.3) assumes that the two particles are far away compared with the amplitude of their oscillations 2

*R*, so that the separation between the particles can be approximated to the constant .

_{g}Introducing radial driving forces in equation (4.2) and using equation (4.1) leads to:
4.4with
4.5and
4.6Here, *ε* = 3*a*/(4*d*), *γ* = 6*π*η*a* and *F*(*r*) is the force from the driving potential, which has the same shape from *U*(*r*) ∼ *r*^{α} for the two oscillators, so that
4.7

To write equations (4.4)–(4.6), we assumed that the variations in angle in a half-cycle are small, so that the function *h _{rθ}* can be put out of the integral in equation (4.1). Writing the equations to the highest order, the

*t*

_{f}_{,i}and

*t*

_{g}_{,i}integration boundaries can also be replaced by generic variables

*t*and

_{f}*t*corresponding to the times at which an uncoupled oscillator would be respectively at positions

_{g}*R*and

_{f}*R*, and

_{g}*r*=

*r*(

*t*) ≈

*r*(

_{i}*t*) corresponds to the radial position of that uncoupled oscillator. With these approximations,

*Δ*

*θ*

_{0}does not depend on the angles

*θ*

_{1}and

*θ*

_{2}, nor on the oscillator. This is a constant that depends on how the angular freedom is implemented in the model of oscillator. The

*h*(

_{rθ}*θ*

_{1},

*θ*

_{2}) term, however, is related to the variation in the hydrodynamic coupling of a bead moving along the radial direction with the velocity of the other bead along the orthogonal direction, depending on the angles of the oscillators.

To go further in the calculation, we need to include that oscillations are either in phase, if *α* < 1, or in antiphase if *α* > 1. We introduce a variable *δ*, equal to 1 for oscillations in phase and −1 for oscillations in antiphase. The state of synchronization *δ* determines the sign of the product *σ*_{1}*σ*_{2} [38]. More precisely,
4.8with
4.9The quantity *γεh*_{rr} represents the coupling force on bead 1 and projected along , coming from a radial force acting on bead 2 along . Equations (4.4) and (4.8) lead to the iterative map:
4.10

This system of equations can be studied by linear stability analysis for *δ* = −1 and *δ* = 1 (*Δ**θ*_{0} being positive), by separating the regions of different signs. Instead, in figure 4, we simply plot equation (4.10) in the (*θ*_{1}, *θ*_{2}) plane by representing the evolution of an initial condition (*θ*_{1}, *θ*_{2}) as an arrow centred on (*θ*_{1}*θ*_{2}) and of direction (*Δ**θ*_{1}, *Δ**θ*_{2}). For *α* < 1 (*δ* = 1, figure 4*a*), the system converges from any initial condition to (*θ*_{1}, *θ*_{2}) = (0, 0), in agreement with the simulations in figure 3. For *α* > 1, in figure 4*b*, the system moves towards a position on the locus of points defined by *h _{rr}*(

*θ*

_{1},

*θ*

_{2}) = 0 (solid lines). However, because the angular speed does not converge to zero around this line, the system will jump from one side to the other side of the line. This results in a change in sign of

*h*. Therefore, the oscillations that were in the synchronized state in antiphase before the jump, become in phase. The system will tend to return to the stable state of oscillations in phase, but, because, the

_{rr}*rr*-coupling is close to 0 nearby the line, this takes several cycles of oscillations. Therefore, when crossing the zero-

*rr*-coupling line, the angles will move away from the line for several cycles, following, for a while, arrows in the opposite direction as the ones indicated in figure 4

*b*. Once the system has returned in the synchronized state in antiphase, the angles will follow the convergence map for

*α*> 1 again. Therefore, the system constantly oscillates between the two sides of the zero-

*rr*-coupling line, with a large amplitude, related to how fast the oscillations converge to the synchronized state. This agrees again with the simulations in figure 3 in the middle-range for

*R*/

_{f}*R*. When

_{g}*R*∼

_{f}*R*, the assumption that oscillations are exactly in phase or in antiphase becomes wrong, as thermal fluctuations and coupling introduce little delays between the switches of the traps of the two oscillators. When the delays become of the order of the time an oscillator spends in the 0.2 ≤

_{g}*R*/

_{f}*R*≤ 0.8 region, equation (4.8) does not apply.

_{g}To summarize, when *α* < 1, the oscillators align to the direction of highest synchronization strength (or highest *rr*-coupling), whereas for *α* > 1, they take orientations that minimize the synchronization and the system is barely synchronized.

## 5. Linear array of oscillators

When increasing the number of oscillators, various geometrical configurations can be studied numerically; the resultant dynamics of the system can be represented by plotting the angle of each oscillator as a function of time.

Figure 5 shows simulations of 60 oscillators equally spaced by a distance *d* = 10 µm along a line (‘chain’ configuration, (*a*–*d*)) and a circle (‘ring’ configuration, (*e*–*h*)). In the chain configuration, the angles are measured from the direction of the line of oscillators, whereas in the ring configurations, the angle of an oscillator is measured from the tangent to the circle at the position of the oscillator. As for two oscillators, the cases *α* = 0.7 (typical for *α* < 1) and *α* = 2 (typical for *α* > 1) lead to different behaviours. For *α* < 1, two neighbouring oscillators tend to be parallel. In the chain configuration, the system shows a strong alignment of the oscillators along the direction of the line (*θ*_{i} = 0). In the more symmetric ring configuration, the oscillators tend to align tangentially to the circle (*θ*_{i} = 0 again). The alignment is weaker in the ring configuration, as shown by the width of the peak in the distribution of the angle displayed as insets in figure 5 (distributions over time and oscillator index). In both chains and rings, the oscillators align in the configuration that maximizes the hydrodynamic *rr*-coupling. For *α* > 1, the graphs are more ‘granulated’: neighbouring oscillators tend to minimize their coupling, leading to different angles between consecutive oscillators. However, the average distribution shows a single, very wide peak at *π*/2, orthogonal to the direction of the line in the chain and orthogonal to the tangent to the circle in the ring. In all cases, the thermal noise has little effect on the width of the peaks in the angle distribution for the parameters used here.

## 6. Two-dimensional arrays of oscillators

Cilia in biological systems are usually arranged on two-dimensional carpets rather than chains or rings. To capture this, rectangular *N _{x}* ×

*N*arrays of about 64 oscillators with square and hexagonal lattices are simulated. Figure 6 shows phase angles for two lattices, with

_{y}*T*= 296 K. The distance between neighbouring oscillators is kept the same between the two lattices, and only simulations for

*α*< 1 are shown, because

*α*> 1 does not display strong alignment properties. The angles are measured from the

*x*-axis, and the oscillators are indexed as indicated in figure 6

*a,d*. The first result is that the oscillators show a cooperative behaviour, and tend to align with the same angle at a given time in the four configurations studied. This is even visible in the highly symmetric 8 × 8 square lattice simulation that shows no preferred average direction of oscillation. The second result, shown by the angle distributions in the insets, is that except for the 8 × 8 square lattice simulation (figure 6

*b*), the system has a preferred direction of either 0 or

*π*/2 rad. The preference for one direction could have two origins: the boundary conditions (shape of the surface covered by the array), or the type (square or hexagonal) and orientation of the lattice. The surface of the array is an

*N*×

_{x}d*N*rectangle for the square lattice and an rectangle for the hexagonal lattice. Therefore, the rectangular surface is stretched along

_{y}d*y*in (

*c*) and (

*f*), stretched along

*x*in (

*e*) and a square in (

*b*). It follows from the angles distributions that the oscillators align along the direction of stretching of the array, which is a configuration of higher

*rr*-coupling than the other axis of the rectangle. This behaviour is inherent to the higher hydrodynamic coupling in the Oseen tensor in the

*x*-direction (3

*a*/(2

*d*)) than

*y*direction (3

*a*/(4

*d*)), for two particles at positions on the axis. The type of the lattice has a negligible effect on the alignment compared with the aspect ratio of the array.

## 7. Discussion and conclusions

We implemented a model oscillator that has the ability to change its direction of oscillation (the beating plane), and aims to describe the collective behaviour of multiple cilia subject to hydrodynamic interaction. The oscillation direction is free: on each cycle, it adjusts by an angle proportional to the velocity of the fluid flow projected orthogonally to the oscillator. While we do not discuss how this freedom is allowed in real cilia, many models can lead to a such response of the angle to an external flow. We chose a specific model in order to calculate its emergent properties, but the details of how angular freedom is implemented are not crucial to the results: the quantity *Δ**θ*_{0} in equation (4.6) depends on the details, but this enters as a constant in equation (4.4), setting the amplitude of the response to the flow. It is the structure of this equation that underlies our results, so other ways to implement flexibility would lead to similar conclusions as the ones we draw in this paper.

The alignment properties of a system of oscillators depend on the details of the driving force, which are matched very simply to the synchronized state of a system of two oscillators. In this paper, the synchronized state is tuned by the parameter *α*, which characterizes the driving force, and is either in phase or in antiphase. An empirical rule, that emerges from all the simulations in this paper, is that a given configuration of oscillators tends to put itself in angular configurations of highest (for in-phase synchronization) or lowest (in antiphase synchronization) coupling (in absolute value) between the oscillators, projected along their directions of oscillation (*rr*-coupling). This rule is confirmed by a theoretical model in the case of two oscillators, for which, the locus of maxima of *|h _{rr}|* is reduced to a single point (

*θ*

_{1},

*θ*

_{2}) when

*α*< 1, and the locus of minima is a curve when

*α*> 1. For that reason, and because synchronization tends to be lost near the zero-

*rr*-coupling line, only oscillators with an in-phase stable state display strong alignment. This study provides a link between the beating pattern of the cilia (widely believed to determine the synchronized state) and the orientation in arrays of cilia. It is interesting to note here that ‘real’ cilia actually beat in phase with their neighbours. This makes the case

*α*< 1 more relevant for cilia, and the study shows that it is the case that leads to alignment. We also suggest that the loss of polar order in motile mutant cilia could be due to a modified driving potential acting on them.

While we do not define precisely what ‘highest’ and ‘lowest’ coupling means in a system of more than two oscillators, simulations with large number of oscillators seem to confirm a similar rule that could be used to predict the state of alignment of a system. We speculate that the rule is related to how the relaxation times of the normal modes of oscillation depend on the angles, as the normal modes already play a key role in determining the state of synchronization, when no angular freedom is added [38]. In addition, in an array of cilia, the synchronized state as a metachronal wave can be seen as a state minimizing the energy required for the beating [25]. Similarly, the rule aligning the oscillators could therefore be related to a problem of energy minimization or maximization.

Alignment properties in large arrays is described both in terms of cooperative behaviour between oscillators and in estimations of the preferred direction of oscillation (if any). The rule sets conditions on both properties. States with all oscillators aligned in the same direction at a given time have high *rr*-coupling and are therefore seen in all the simulations for *α* < 1, leading to high cooperativity between the oscillators. Non-symmetric configurations such as chains of oscillators or rectangular arrays also tend to confine the angles in a particular direction, which is the direction in which the system is elongated. In two-dimensional arrays, the type of lattice could also affect the alignment of oscillators. However, it appears that in arrays of below approximately 100 oscillators, the effect of the shape of the array is more important than the lattice to determine the preferred direction of oscillation.

In this paper, we have used parameters close to motile cilia conditions. The size of the array varies a lot with the biological system: thousands of cilia can be arranged in a dense array such as in the alga *Volvox carteri* [11], or in *Paramecium* [44]; or a few hundred cilia can be packed on the surface of a multi-ciliated airway cell, this cluster interacting with the clusters on other cells in the tissue [18]. In all these cases, the boundary conditions could have a determinant role in the choice of the direction of alignment.

A *quantitative* connection to the biological question of how airway epithelium gains its full organization in development (once cells are elongated) will clearly need further development of this or related models. In this biological system, the cells are elongated along the proximal/distal axis, conferring polarity to the airway tissue as a whole, and the cilia within each cell become aligned to each other. Furthermore, they become aligned with the tissue axis. The main aspects that have been neglected in our model are the likely viscoelastic properties of the fluids, and the effects on the fluid flow due to the planar surface on which the cilia are anchored. Adding a surface changes the form of the velocity decay away from a cilium (and hence the hydrodynamic coupling between pairs) from 1/*r* in bulk to 1/*r*^{3} with *r* the distance between the two particles. While this is not expected to change much our results on the alignment of two oscillators, because the distance between the beads stays almost constant, it does change the hydrodynamic coupling from a long-range interaction in bulk to a shorter range interaction in the presence of a surface. In addition, as a consequence of the small size of cilia, the Brownian motion leads to non-negligible fluctuations of the phase of the oscillators. However, the simulations with a large number of oscillators suggest that the thermal noise is irrelevant when looking at the alignment properties of the oscillators.

The current model and the behaviour highlighted in this paper show that hydrodynamic coupling is able, at least in principle and qualitatively, to lead to orientation of active oscillators, and that the resulting collective dynamical state responds to the symmetry of the shape of the array (symmetry of the boundary condition). The case of two-dimensional arrays with cilia described by a driving force with *α* < 1 could explain how cilia distributed on the surface of an elongated cell align with each other, and pick the axis of stroke from the cell's elongation. In this speculative scenario, the microtubule networks which link the basal bodies within the cells, and are seen to correlate with polarity [16], might act to freeze in place the orientation of the basal bodies.

The study here proves the possibility of spontaneous alignment, in the absence of an external flow. Various ciliated tissues (described in §1) initially develop in the absence of flow. The presence of a directed external flow (which is possible in some systems, such as the developing airways) would be an even stronger aligning factor. In general, the importance at various stages of development of flow and mechanical force transduction, which are undeniably present between cilia, remain an open question clearly requiring multidisciplinary approaches.

## Funding statement

This work is supported by the Marie Curie grant ITN-COMPLOIDS (FP7-PEOPLE-ITN-2008 no. 234810).

## Acknowledgements

We thank M. Cosentino Lagomarsino for stimulating this work and providing a critical review.

- Received June 26, 2013.
- Accepted July 2, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.