Physical limits of flow sensing in the left-right organizer

Fluid flows generated by motile cilia are guiding the establishment of the left-right asymmetry of the body in the vertebrate left-right organizer. Competing hypotheses have been proposed: the direction of flow is sensed either through mechanosensation, or via the detection of chemical signals transported in the flow. We investigated the physical limits of flow detection to clarify which mechanisms could be reliably used for symmetry breaking. We integrated parameters describing cilia distribution and orientation obtained in vivo in zebrafish into a multiscale physical study of flow generation and detection. Our results show that the number of immotile cilia is too small to ensure robust left and right determination by mechanosensing, given the large spatial variability of the flow. However, motile cilia could sense their own motion by a yet unknown mechanism. Finally, transport of chemical signals by the flow can provide a simple and reliable mechanism of asymmetry establishment. DOI: http://dx.doi.org/10.7554/eLife.25078.001


Introduction
Chirality describes the asymmetry between an object and its mirror image. How chiral asymmetries arise in physics and biology is a fundamental question that has fascinated researchers for many years (Morrow et al., 2017;Wagnière, 2007). A number of independent mechanisms of chirality establishment at multiple scales have been identified in living systems, ranging from subcellular with the establishment of chiral cortical flows (Naganathan et al., 2014) to the cellular scale during the process of cell fate specification (Gó mez-Ló pez et al., 2014) and to the tissue scale during the process of left-right (LR) axis specification (Blum et al., 2014b;Coutelis et al., 2014;Dasgupta and Amack, 2016;Hamada and Tam, 2014;Levin, 2005). Similar concepts also exist in the plant kingdom, suggesting widespread mechanisms of chirality transfer (Wang et al., 2013). Biological symmetry breaking is often associated with the internal handedness or chirality of molecular motors (Ferreira and Vermot, 2017;Inaki et al., 2016;Naganathan et al., 2016). One of the most striking examples of biological symmetry breaking is the mechanism of LR axis determination in the developing embryo, which is crucial for the asymmetric internal organ positioning. In most of the vertebrate species, LR axis determination is set by a symmetry biasing event which is under the control of a directional flow generated by the chiral beating pattern of motile cilia (Hirokawa et al., 2009;Nakamura and Hamada, 2012). To date, even though the molecular mechanisms of LR axis determination are highly conserved from fish, amphibians to mammals (Blum et al., 2014a), many steps of this process are not yet understood (Ferreira and Vermot, 2017;Shinohara and Hamada, 2017;Wolpert, 2014).
Prevailing models of LR specification in most vertebrates Gros et al., 2009) involve groups of cells (the LR organizer, LRO) within the presomitic mesoderm (the segmental plate) coordinating asymmetry establishment through the control of a slow-moving flow (the nodal flow) (Nonaka et al., 1998) and an intercellular amplification of the asymmetric signals within and around the LRO cells (Nakamura et al., 2006). This slow-moving flow is produced by rotation of multiple, motile cilia located at the cell surface of the LRO. This flow leads to a collective cell response that occurs specifically on the left embryonic side of the LRO and is associated with an asymmetric intracellular calcium release (McGrath et al., 2003;Yuan et al., 2015). The LRO is a transient structure and its function is required for a limited period of time. In zebrafish, the LRO is visible between 10 hr after fertilization (3-somite stage) and 16 hr after fertilization (14-somite stage) (Essner et al., 2005;Yuan et al., 2015).
Cilia generate a directed flow by beating in a spatially asymmetric fashion (Satir and Christensen, 2007). Their cycle consists of a working stroke, during which a cilium is stretched away from the surface in order to move the maximum amount of fluid, and a recovery stroke, during which it swipes along the surface, thus reducing the backflow (Marshall and Kintner, 2008). Cilia in the LRO are relatively short and their beating pattern resembles rotation along the mantle of a tilted cone (Hirokawa et al., 2006;Nonaka et al., 2005;Okada et al., 2005). This tilt effectively makes one part of the cycle act as a working stroke and the other part as a recovery stroke. Spatial orientation is, therefore, a key functional feature of motile cilia involved in LR symmetry breaking, as it determines the strength and directionality of the induced flow (Cartwright et al., 2004). The molecular mechanisms that set this orientation involve elements of the Planar Cell Polarity (PCP) pathway (Borovina et al., 2010;Hashimoto et al., 2010;Song et al., 2010) as well as flow itself (Guirao et al., 2010). At the molecular scale, the chirality of biomolecules has been proposed to provide some of the asymmetrical cues for LR symmetry breaking (Brown and Wolpert, 1990;Levin and Mercola, 1998). This paradigm works well with cilia, as the sense of cilia rotation is determined by the chirality of the structure of their protein building blocks (Hilfinger and Jülicher, 2008).
To accurately discriminate left from right, cells need to robustly sense a signal over the noise associated with flow and cilia beat. Two hypotheses have been proposed for asymmetric flow detection. According to the chemosensing hypothesis, the directional flow establishes a LR asymmetric chemical gradient that is detected by signaling systems which leads to LR asymmetric gene expression and cell responses in the LRO (Okada et al., 2005). The mechanosensing hypothesis, on the other hand, proposes that the LRO cells can detect the mechanical effects of flow. It has been suggested that this mechanosensing is mediated by a particular type of sensory cilia that is able to trigger a local, asymmetric response of the so-called crown cells, which are located at the periphery of the node (mouse LRO) (McGrath et al., 2003;Tabin and Vogan, 2003). In the zebrafish LRO (also called Kupffer's vesicle or KV), the model for cilia-mediated mechanotransduction suggests that sensing cilia are immotile and will activate a cellular response due to physical asymmetries generated by the flow (Sampaio et al., 2014). The flow detection apparatus needs to be remarkably efficient, as demonstrated by the observation that, despite wild-type mice having hundreds of motile cilia in their LRO, proper asymmetry establishment occurs even in mutant mice with only two motile cilia in the LRO . Similarly, although wild-type zebrafish usually has around 50-60 motile cilia in the KV, correct asymmetry establishment occurs in zebrafish mutants with only 30 motile cilia (Sampaio et al., 2014). Considering that motile cilia might be chemosensors (Shah et al., 2009), LR symmetry breaking could rely on a combination of both mechanosensory and chemosensory mechanisms that could work in parallel. To date, however, several issues with both hypotheses have emerged: the molecular nature of a possible diffusing compound involved is still not known (Freund et al., 2012;Shinohara and Hamada, 2017) and cilia-mediated mechanotransduction in the mouse node has recently been experimentally challenged (Delling et al., 2016). In many aspects, these debates remain open because the sensitivity of mechanical and chemical detection mechanisms has not yet been assessed quantitatively from a physical standpoint based on properties of cilia and the flows they generate.
To test the sensitivity of the detection mechanisms, we investigated the physical limits of the system to discriminate left and right in the presence of flow irregularities and noise. We analyzed the system from the level of individual cilia to the scale of the entire organ using live imaging to determine the physical features controlling the flow. We used these experimental datasets to calculate the flow in unprecedented detail and assess its robustness. We show that the flow velocities and their local variability impose crucial constraints on potential mechanosensing by cilia and investigate the question whether the small number of non-motile cilia can be sufficient to reliably distinguish between the left and the right side of KV. Furthermore, we use the calculated flow profile to simulate the directed diffusive motion of chemical signals and determine the limit on particle diffusivity for which the mechanism is reliable. Our results show that the physical limits to the reliability of mechanical sensing of flow are stronger than the ones associated with chemical sensing and suggest that chemosensation is the key mechanism for LR axis determination.

Results
Theoretical analysis of cilia-generated flow patterns In the current models, cilia-mediated left-right (LR) symmetry breaking is driven by a chiral flow pattern. In the zebrafish embryo, the directional flow corresponds to an anti-clockwise rotation around the dorsoventral (DV) axis when viewed from the dorsal pole of the Kupffer's vesicle (KV) ( Figure 1A) (Essner et al., 2005;Kramer-Zucker et al., 2005). Motile cilia are the driver of the flow whose directionality mainly depends on two features: cilia density and spatial orientation. Several studies have focused on cilia orientation and density to explain the flow directionality in mouse and fish (Borovina et al., 2010;Cartwright et al., 2004;Montenegro-Johnson et al., 2016;Nonaka et al., 2005;Okabe et al., 2008;Okada et al., 2005;Sampaio et al., 2014;Smith et al., 2008) to understand the principles underlying flow generation. These studies led to contradicting conclusions by proposing different types of cilia orientation, such as posterior tilt (Borovina et al., 2010), dorsal tilt (Supatto et al., 2008), or a mix of the two (Okabe et al., 2008). To solve this issue, we used a theoretical approach and developed a simplified biophysical model that can be tested in vivo.
In order to gain insight into the mechanism of flow generation, we first use a simplified model, which does not consider the small scale inhomogeneities of the flow around individual cilia. This simplification is valid in the limit in which both the cilia length (L) and the characteristic distance between cilia are shorter than the distance at which we observe the flow. The latter distance is characterized by the radius of the KV (R). The cilia layer can then effectively be represented by a net slip velocity at the surface (Vilfan, 2012). With typical parameters L = 6 mm, R = 35 mm and a characteristic distance between cilia of 10 mm, these conditions are roughly satisfied when observing the flows in the center of the KV. Calculations with a detailed hydrodynamic model shown later (section 'Single vesicle analysis reveals a significant variability between embryos') confirm uniform flows in the center, but they also show a significant velocity variability near the KV surface, which is not captured by the simplified model.
The effective surface-slip velocity, which replaces the individual cilia in the simplified model, can be calculated from the cilia parameters as follows. We consider cilia covering a surface with area density r and rotating with the angular velocity w along the mantle of a cone with a semi-cone angle , tilted by the angle in a directionẽ t . The direction of rotation is clockwise as seen from the distal end towards the KV wall. They induce a net flow velocity above the ciliated layer in the direction perpendicular to the tilt direction (Smith et al., 2008;Vilfan, 2012) whereẽ n is a unit vector normal to the vesicle surface, Â is the vector product, C N » 1:2ph denotes the drag coefficient and h the fluid viscosity (see Figure 1C). This contribution results from the fact that a cilium moves more fluid during its working stroke, when it is further away from the surface, than during the recovery stroke, when it is closer. At the same time, even without a tilt (i.e. when the rotation axis is orthogonal to the cell surface), a rotating cilium induces a rotary flow with the ampli-tude~3C N !L 4 = 16ph ð Þ sin 2 ð Þ cos ð Þ. Although a surface (or a cavity) lined with untilted rotating cilia at a uniform density does not produce a long-range flow, non-uniformities in surface density r do lead to an effective slip velocity The flow observed inside the KV, which is, as a first approximation, characterized by uniform rotation (the fluid moves like a rotating rigid sphere) with the angular velocity W and the surface velocitỹ Figure 1B for the definition of the local basisẽ f ;ẽ n ;ẽ m À Á and Figure 1D for the definition of b) can be achieved in two ways (or a combination thereof): 1. Scenario 1: Dorsoventral gradient of cilia density (Figure 2A).
Both scenarios or a combination of both are possible. However, the fluid velocity achieved with scenario 1 is smaller by a factor of 3L=8R ð Þ cos ð Þ= sin 0 ð Þ » 0:1 as compared to scenario 2. This theoretical analysis shows that the directional flow depends on two main topological features of the vesicle: a profile of cilia density and a cilia orientation pattern following a meridional tilt in the KV.

3D-CiliaMap reveals multiple gradients of cilia density
In order to test whether scenario 1 (the cilia density gradient scenario, Figure 2A) or 2 (the meridional tilt scenario, Figure 2B) leads to the directional flow observed in vivo, we developed a live imaging-based method called 3D-CiliaMap ( Figure 1-figure supplement 1 and Video 1). 3D-CiliaMap is designed to quantify experimental features, such as KV size and shape, as well as the spatial distribution, surface density, motility and orientation of cilia ( Figure 1-figure supplement 1). Developing zebrafish embryos can be accurately staged using the number of somites (blocks of presomitic mesoderm tissue that regularly form along both sides of the neural tube during the segmentation period) (Kimmel et al., 1995) ( Figure 1A). We performed our analysis between 8-and 14-somite stages (SS), when the directional flow in the KV is well established (Essner et al., 2005;Kramer-Zucker et al., 2005;Long et al., 2003;Lopes et al., 2010). We analyzed the cilia density in 3D because it has previously been reported to vary along the anteroposterior (AP) axis (Borovina et al., 2010;Kreiling et al., 2007;Okabe et al., 2008;Supatto et al., 2008;Wang et al., 2011Wang et al., , 2012. In order to extract average features, the density maps from 20 vesicles with a total of 1197 cilia were averaged and represented either on the average vesicle spheroid ( Figure 3A) or on a 2D flat map ( Figure 3B). We found that the average density of cilia did not display any significant differences between left and right sides of the KV ( Figure 3B). As expected from previous studies (Wang et al., 2011(Wang et al., , 2012, a steep density gradient was observed along the AP axis, with increasing density towards the anterior pole of the vesicle (about four times denser than posterior) ( Figure 3A,B). Interestingly, a gradient of cilia density was also seen along the DV axis, albeit shallower (dorsal being two times denser than ventral) ( Figure 3B). Together, these results rule out scenario 1, since the dominant AP density gradient would lead to a significant contribution to directional flows around the AP axis, which are not observed (Supatto et al., 2008).

KV motile cilia exhibit a meridional tilt
To test the scenario 2 (the meridional tilt scenario, Figure 2B) in vivo, we mapped cilia orientation in embryos between 8-and 14-SS. We defined the orientation of a cilium with a unit vector along the axis of the conical cilia movement and decomposed it in a local orthogonal basis e f ;ẽ n ;ẽ m À Á on the KV surface as defined in Figure 1B, where: 1. q (tilt) is the angle of the cilium with respect to the KV surface normal (0˚for a cilium orthogonal to the KV surface and 90˚for parallel, Figure 1C). 2. j angle is the orientation of the cilium projected on the KV surface (0˚for a cilium pointing in a meridional direction towards the dorsal pole, Figure 1C).
A meridional tilt would then correspond to q > 0˚and j close to 0˚. Among the 1197 cilia in 20 vesicles, we could determine their orientation and motility status for 86% and 89% of them, respectively. Less than 5% of the cilia were immotile at these stages. Together, we quantified the orientation of 1013 motile cilia (corresponding to 85% of all cilia). We plotted cilia unit vectors in the same local basis (black arrows in Figure 3C) and the average cilium as the 3D mean resultant vector (Berens, 2009) (gray arrow in Figure 3C). The resultant vector length r quantifies the spherical spread (the closer r is to one, the more cilia are concentrated around the mean direction). Despite a relatively broad distribution of cilia orientations (r ¼ 0.76), we found that the q tilt of the average cilium is near 46˚and the average cilium has a 14˚j angle ( Figure 3C). The resultant vector length and angles (r, q, j) from individual vesicles are listed in Table 1. In total, 65% of all motile cilia exhibited a j angle between À45˚and +45˚(meridional quadrant). In addition, no specific spatial distribution of this orientation was detected (red in Figure 3D). These data demonstrate that motile cilia are oriented on average with a meridional tilt and support scenario 2 ( Figure 2B), in which cilia meridional tilt is the dominant mechanism generating the directional flow within the KV ( Figure 2C). Note that the meridional tilt results in dorsal cilia orientation on the equator and in posterior cilia orientation in the anterior part of the dorsal KV. This might explain why several studies reported different orientation of the cilia in the KV (Borovina et al., 2010;Okabe et al., 2008;Supatto et al., 2008). With Table 1. Statistical properties of all KV analyzed.  Variation in cilia surface density over time affects flow amplitude but not its direction Current models for symmetry breaking are not taking into account the dynamics of cilia spatial distribution during the process. Nevertheless, the process of LR patterning occurs in a dynamic organ (Compagnon et al., 2014;Wang et al., 2012;Yuan et al., 2015), which might be associated with changes in the spatial distribution and orientation of cilia during the course of LR patterning. Previous reports have shown that the first signs of asymmetric cell response in the KV are observed between 3-and 8-SS (Francescatto et al., 2010;Sarmah et al., 2005;Yuan et al., 2015). To obtain accurate information about these dynamics, we map cilia positioning within the KV in three pools of embryos: early (3-SS), mid (8-SS) and late (9-14-SS) ( Figure 4A). During this developmental time window, the KV size drastically changes with an average volume increasing from 17 pl at early to 45 pl and 82 pl at mid and late stages. The average number of cilia per vesicle also increases from 30 at early, to 52 and 63 at mid and late stages (Figure 4-figure supplement 1 and Table 1). Not surprisingly, the ratio of motile versus non-motile cilia also changes along with the total number, as previously observed by (Yuan et al., 2015).  Table 1). At mid and late stages, we observed a steep gradient  Figure 5 gðgÞ Cilia distribution, pair correlation + See Figure 5 h Fluid viscosity 0:001 Pa s r Stokes Diffusive particle Stokes radius 0:5 À 10 nm of cilia surface density along the AP axis, which is not yet established at 3-SS ( Figure 4-figure supplement 2B). Instead, a DV gradient of cilia density is present at 3-SS. We next investigated the emergence of cilia motility and orientation by comparing early, mid and late stages. We observed that the spread around the average motile cilium is constant over time (resultant vector length r = 0.74 or 0.75) ( Figure 4A). Additionally, it seems that the average orientation angle q of the motile cilia increases over time (from 32˚to 50˚, Figure 4A). When focusing on 3-SS embryos, we found that motile cilia at this stage already exhibit a clear meridional tilt, as the angle j of the average motile cilium is close to 0˚(57% of all cilia exhibited a j angle in the meridional quadrant [À45˚, +45˚]) and q is high (32˚) ( Figure 4A). Together, these data suggest that the few motile cilia at 3-SS are already well oriented and generate a flow of low amplitude but in the proper direction.
We suspected that the increasing number of motile cilia and the changes in their spatial distribution would significantly alter the flow profile between 3-and 8-SS. To test this, we used our 3D cilia maps to numerically calculate the flow they generate at the different developmental stages ( Figure 4B and Video 2). We first validated our flow simulation by comparing it with experimentally measured flow profiles along the AP axis ( Figure 5). As expected, due to the increase in anterior cilia density, the flow amplitude increases over time ( Figure 4B and Figure 6A) with the most pronounced increase in the anterior region. We quantified the directionality of the flow by calculating the effective angular velocity, defined as the angular velocity of a uniformly rotating sphere with the same angular momentum as the circulating fluid in the vesicle ( Figure 4C and Table 1). Our results consistently show that in 8-SS and 9-14-SS embryos, the rotational flow persistently points in one direction. Most importantly, in 5 out of 8 vesicles, the flow direction is already set at 3-SS ( Figure 4C and Figure 6-figure supplement 1A), suggesting that directional flow can emerge as early as 3-SS, even though the flow is of low amplitude (Video 2). As a consequence, the establishment of an AP gradient of cilia from 3-to 8-SS does not affect the flow direction, which further supports cilia meridional tilt as the dominant mechanism used to generate directional flow within the KV.

Single vesicle analysis reveals a significant variability between embryos
An additional element to consider in the physical mechanisms of symmetry breaking is its robustness. Given that 90-95% of the zebrafish embryos have a properly positioned left axis (Gokey et al., 2016), the mechanism eliciting LR bias has to be highly robust even though the KV size, which affects LR patterning, is variable across embryos (Gokey et al., 2016). We sought to directly probe for the robustness of the biophysical features of the cilia through single vesicle analysis. We first analyzed cilia density in individual vesicles and found that cilia density and orientation are very variable from embryo to embryo ( Figure 6B We conclude that the local flow amplitude itself cannot be a good indicator of the embryonic side and is too variable to serve as a robust predictor of the left and right side of the vesicle. By contrast, the strength of the rotational flow characterized by the effective angular velocity appears robust ( Figure 4C). Thus, while the average vesicle highlights a highly stereotyped organization, single vesicle analysis uncovers a high diversity of densities, orientations and local flow profiles. This identifies the directional flow as the most robust left-right asymmetric feature in the vesicle. Comparing flow-mediated transport of signaling molecules and mechanical flow sensing as mechanisms for symmetry breaking in the KV The hypotheses for the mechanisms of LR symmetry breaking are usually divided into two competing classes, independently of the topological differences amongst species: asymmetric distribution of signaling molecules (Hirokawa et al., 2006;Okada et al., 2005) or asymmetric mechanical influence (Hamada and Tam, 2014;Yoshiba and Hamada, 2014). Making use of the information gained from 3D-CiliaMap and numerically calculated flow profiles, we developed numerical simulations based on specific arrangements of cilia in order to corroborate or refute different hypotheses for symmetry breaking in the LRO. In particular, we took into account the variability between embryos that we observed experimentally.

Mechanosensory mechanism 1: directional flow sensing
We first evaluated the 'Directional flow-sensing mechanism' where the left side is detected by sensing the directionality of the circular flow. Such a mechanism implies that cells need to distinguish between a fluid moving from anterior (A) towards posterior (P) (on the left side) and a fluid moving from P towards A (right side), and the cilium needs to sense the direction of flow. It also has to overcome the following obstacles: (1) The strong temporal variation associated with beating cilia overlays the directional flow to be detected; (2) In addition, the flow fields of adjacent cilia also perturb the static component of the directional flow; (3) The cilia would need a detection threshold sufficiently low to detect the relatively weak flows in the KV. Using our computational modeling, we tested the feasibility of criteria 1-3.
We used the cilia distributions from 20 analyzed vesicles, as well as a larger number of randomly generated vesicles with the same density profiles and interciliary distance distributions (see Materials and methods). In each vesicle, we randomly chose three immotile cilia on each side and calculated the torque acting on them. The torques acting at the base of a cilium are calculated using the numerically determined force distributions along its length. The torque component that could potentially serve as the LR determinant is the meridional one, bending the cilium towards A or P ( Figure 7A). In order to estimate the fraction of cilia that measure the torque in the correct direction above a certain threshold, we plotted cumulative distributions of the meridional component of the torque vector ( Figure 7B,C). The instantaneous torques (dotted lines) show a broad distribution and have a direction opposite to that of the directional flow in about 25% of the cases. These fluctuations are caused by the beating of adjacent cilia, which induces on average an oscillating torque with an r. m.s amplitude of 9 Â 10 À19 Nm. Additionally, the torques are of a similar order of magnitude as the thermal fluctuations acting on the cilium, whose r.m.s. amplitude we estimate as 2 Â 10 À19 Nm (Appendix B). A simple estimate shows that both the thermal and oscillatory noise can be suppressed by temporal averaging (low-pass filtering) the signal with a time constant longer than 2 s. The distribution of time-averaged torques is shown by dashed lines in Figure 7B and C. A sensitivity threshold of 2 Â 10 À19 Nm would be sufficient to achieve a reliability of 95% (less than 5% of the immotile cilia on the left side would not detect the posterior-directed flow) in vesicles with the generated distributions ( Figure 7B, dashed lines), but the greater variability of experimentally characterized vesicles does not allow this level of reliability (about 10% of the cilia are subject to flows of the opposite directionality, Figure 7C, dashed lines). A sufficient reliability can only be achieved by additionally averaging the torques detected on all three cilia on one side (solid lines). In this case the required detection threshold is 10 À19 Nm. For comparison with other flow sensing cilia, we calculated the uniform shear rate required to exert the same torque. We find that for a 6 mm long cilium the threshold torque corresponds to a shear rate of 0.5 s À1 (shear stress 0.5 mPa). In renal cilia, Rydholm and colleagues (Rydholm et al., 2010) observed calcium signals with shear stresses of 20 mPa and higher. These results indicate that the mechanosensory detection of flow would require cilia with an ability of direction-sensitive flow detection with a threshold 1-2 orders of magnitude lower than known comparable mechanosensory cilia. Moreover, many KV have fewer than three immotile cilia on each side, the number that would be needed for ensemble-averaging to overcome spatial inhomogeneities of the flow. The number of immotile cilia is higher at 3-SS, but the weaker and less regular flow excludes reliable side detection at that stage (Figure 7-figure supplement 1). Together, these results suggest that flow sensing in itself is difficult as the flow is weak and masked both by spatial and temporal fluctuations. We thus expect that such a mechanism of sensing would lack the robustness necessary for setting the LR axis accurately.

Mechanosensory mechanism 2: cells sense the motion of their own cilia
An attractive possibility based on mechanosensing is that cells 'detect the motion of their own cilia relative to already established body axes'. In this case, the rotating cilium provides a cell the necessary chirality information, while the directional flow in KV would only appear as an epiphenomenon. The tip velocity of a beating cilium (400 mm/s) is significantly higher than the typical flow velocity (10 mm/s), which implies that in a motile cilium the torques caused by its own motion largely surpass those caused by the directional flow. Using the mobility matrix of a model cilium (see Materials and methods) we calculated the torque components acting on the base of an isolated dorsally tilted motile cilium, positioned either on the left or on the right side ( Figure 7D). The dashed lines show the time averages of the three torque components. The average torque caused by an active cilium's motion is about 20 times higher than the torque caused by the directional flow (Mechanosensory mechanism 1). If a cell could discriminate between a torque towards A or P exerted at the base of its motile cilium, the time-averaged meridional component of the torque vector (blue line), which has a magnitude of about 10 À17 Nm, could serve as a side discriminator. Thus, cells sensing the torque direction generated by their own cilia to is a possible mechanism for an asymmetric response in the KV.

Flow-mediated transport of a signaling molecule
As the third mechanism, we investigated the possibility of 'flow-mediated transport of a secreted signaling molecule' in the KV. It is known that classical motile cilia also contain receptors to detect the external chemical environment (Shah et al., 2009). We propose that cilia on the left and on the  Figure 7 continued on next page right side (45˚around LR axis, representing the areas where the first asymmetric responses have been observed (Francescatto et al., 2010;Sarmah et al., 2005), act as detectors that absorb small particles in contact with their surface. Because the first asymmetric signal was observed on the left (Yuan et al., 2015), we propose that these particles are secreted in the anterior region (30˚around the anterior pole) and that cells in this region do not absorb them ( Figure 7E and Video 3). As a rough estimate, we expect that flow-mediated transport requires a Pé clet number Pe ¼ vR = D > 1, which states that advection dominates over diffusion. With ¼ 0:5 s À1 , we get Pe ¼ R 2 6phr Stokes = k B T ð Þ ¼ r Stokes = 0:4 nm. The condition Pe > 1 is thus fulfilled for particles above nanometer size. To assess the feasibility of flow-mediated transport quantitatively, we simulated the diffusion of small particles in the flow fields calculated before ( Figure 7F). We propose that asymmetry is detected based on the cumulative number of particles detected on each side. Therefore, the time course of particle secretion is not important for our arguments. As a measure of asymmetry, we used the number of particles detected on the left side N left , divided by the total number detected on the left and right combined (N left+right ) ( Figure 7G-L). We estimate that the difference is detectable if the ratio between left and right is at least 2:1, that is, if the fraction is higher than 2/3 ( Figure 7G-L). Computational results on 20 vesicles are shown in Figure 7I. For signaling particles with a Stokes radius of 2 nm, we obtain a robust asymmetric readout in 18 of 20 vesicles. In addition, we performed the computation on a larger number of randomly generated vesicles and see that, for particle sizes of 2-10 nm, more than 95% of the vesicles fulfill the asymmetry requirement. In order to elucidate the requirements on the arrangement of cilia needed for robust asymmetry establishment, we also simulated vesicles with alternative distributions. Figure 7J shows a scenario with cilia distributed homogenously around the whole vesicle, which is equally efficient. It is therefore unclear whether the increased cilia density in the anterior region fulfills a purpose with regard to symmetry breaking. It is possible, however, it becomes beneficial for different secretion scenarios. Since it was previously shown that there is a minimum number of motile cilia required to achieve LR patterning in vivo (Sampaio et al., 2014), we tested the model with a smaller number of motile cilia. With half the number of cilia (N c ¼ 35, Figure 7K) the vesicles still show sufficient average asymmetry, but fail to achieve 95% reliability. Interestingly, the mechanism becomes dysfunctional with 20 motile cilia, which is very close to the minimum observed in vivo by (Sampaio et al., 2014) ( Figure 7L). This indicates that 'flow-mediated transport of a secreted signaling molecule' is plausible in the KV and that the KV contains enough motile cilia to allow robust symmetry establishment through asymmetric transport of signaling molecules, but not much more than needed.

Discussion
Here, we have investigated physical mechanisms that could underlie the left-right (LR) symmetry breaking of the vertebrate body plan. We have combined microscopy and large-scale cilia mapping in living animals with theory to analyze the physical limits of the generation and detection of LR asymmetric flows. The experimental analysis of cilia patterns over time allowed us to generate  Figure 7C, but using cilia maps at 3-somite stage. DOI: 10.7554/eLife.25078.016 Video 3. Simulated transport of signaling molecules in the Kupffer's vesicle (KV). Panoramic view of the KV as seen from the center. The cilia distribution is obtained from vesicle 1 at 9-14-somite stage ( Figure 6-figure  supplement 1C). Cilia shown in blue are motile and those in red immotile or undetermined. Signaling particles (r Stokes ¼ 5 nm, not to scale), secreted from anterior (yellow), are subject to Brownian motion biased by the leftward flow until they are absorbed by a cilium (shown yellow after particle capture). 10 s of video represent 1s real time. Link for 360˚video: https://youtu.be/1caSzBIe5rA DOI: 10.7554/eLife.25078.017 a comprehensive map of complex in vivo cilia behaviors at key stages of LR symmetry breaking. We used this approach to assess the reliability by which asymmetric cues can be detected by cells in the process of the determination of the LR axis of the vertebrate body plan. Recent work in fish and mice has highlighted two possible ways to mechanically sense the directional flow: asymmetric flow velocity sensing (Sampaio et al., 2014) and/or flow direction sensing (McGrath et al., 2003). Our results argue against asymmetric velocity sensing because no robust asymmetries in the flow amplitude were observed in the KV, as previously suggested . The calculated flow pattern obtained when using cilia parameters collected in vivo shows that its directionality is the most robust sign of chirality in the organizer. Another level of complexity, which is shared between mouse and fish, is the presence of beating cilia in the near vicinity of the immotile cilia. It has been suggested that these cilia detect the flow direction . Indeed, the average distance from one cilium to its nearest neighbor is under 8 mm in the KV. At such a distance, the strong local flow perturbations have been shown to be chaotic (Supatto et al., 2008). In addition, the required sensitivity is of the same order of magnitude as thermal fluctuations. Both cilia beating and the thermal fluctuations could theoretically be overcome by time averaging if mechanosensitive cilia acted as low-pass filters. The observed average velocities are still below the detection threshold found in other mechanosensitive cilia (Delling et al., 2016;Goetz et al., 2014;Nauli et al., 2008;Rydholm et al., 2010) and we are not aware of systems that detect static forces with such sensitivity. Hair bundles are known to achieve remarkable sensitivity thresholds close to the level of thermal noise, but the majority of them detect stimuli at higher frequencies and therefore act as high-or band-pass filters (Muller et al., 2016). A notable exception is the utricular otolith, which acts as a low-pass filter, but senses significantly higher forces (Inoue et al., 2013). We conclude that small flow velocities are a major challenge to the reliable sensing of directionality. Moreover, the hypothesis that cilia are capable of directional flow sensing remains untested in vivo (Shinohara and Hamada, 2017).
In addition, we showed that due to spatial inhomogeneities resulting from the cilia distribution, temporal averaging within a single cilium could not be sufficient for robust side determination via detection of flow direction. Signals from several cilia on each side would need to be ensemble-averaged. For robust flow sensing, we found that three cilia on each side are necessary.
We observed in vivo that the number of immotile cilia decreases over time along with an increase in the number of motile cilia, resulting in a stronger flow velocity and directionality. Nevertheless, the stage at which the number of immotile cilia is high corresponds to the stage where flow is the lowest and the least robust. At later stages, when it directionally is established, we observed that only 5% of cilia were immotile, less than three on each side. The number of immotile cilia is therefore too low to sense the flow direction. Combined with the related finding that cilia from isolated mouse cells do not bend significantly in response to flow applied at endogenous flow amplitudes (Delling et al., 2016), we conclude that a mechanism of symmetry breaking based solely on mechanical sensing of directional flow appears incompatible with our in vivo analysis.
These physical limitations led us to reconsider mechanosensing in the framework of LR symmetry breaking. Mechanosensing could, for instance, occur through the sensing of the cilia's own motion. The mean torque exerted by the fluid on a meridionally tilted beating cilium has the opposite direction from the torque exerted by the directional flow on an immotile cilium. A prediction of this autosensing hypothesis is therefore that an artificially induced counter-rotating flow (clockwise in dorsal view) would enhance, rather than reverse, the LR difference. Results in the mouse node (Nonaka et al., 2002) showing that externally imposed flow of opposite directionality can cause asymmetry reversal, contradict this view. Another observation in Xenopus that cannot be explained by the auto-sensing hypothesis is that the laterality mechanism breaks down when the fluid is made viscoelastic, which brings the directional flow to a halt, even though the cilia are still motile (Schweickert et al., 2007).
Conversely, our numerical simulations show that asymmetric transport of signaling molecules in the KV is a much more robust strategy than sensing mechanical cues. We characterize here the physical limit on the size of the signaling molecule for its reliable asymmetric distribution. Our calculations show that such mechanism requires the particle size to be bigger than 2 nm to work. Interestingly, this is in the size range of membrane-bound extra cellular vesicles (ECVs). ECVs play important roles in intercellular communication and may mediate a wide range of physiological and pathological processes (Cocucci et al., 2009;Hogan et al., 2009;Raposo and Stoorvogel, 2013;Wood et al., 2013).
In summary, we analyzed the physical limits of mechanisms that have been proposed for asymmetry establishment in zebrafish. Combining large-scale in vivo imaging with fluid dynamics calculations we were able to map the biophysical features of cilia in the KV and the flows they generate in unprecedented detail. This allowed us to quantitatively test the physical limits of flow detection mechanisms. We show that the small number of immotile cilia found in the KV cannot be sufficient to robustly detect the direction of the flow given its high local variability. Motile cilia could sense the torques exerted by the fluid as a result of their own motion, which largely surpass the influence of the directional flow. However, this mechanism is incompatible with findings in other vertebrate species in which laterality establishment was suppressed in viscoelastic fluids and reversed with an artificial flow (Nonaka et al., 2002;Schweickert et al., 2007). Finally, we show that a chemosensory mechanism in which a LR gradient is established by combining directional flow around the dorsoventral axis with asymmetric particle secretion in the anterior region could explain the observed robust LR asymmetry establishment, provided that the particle size is above the lower limit of about 2 nm. Although the molecular nature of the flow detection mechanism remains obscure, our analysis of physical limitations of two proposed mechanisms rules out directional flow sensing. It also allows us to predict the minimum size of the signaling particle, which will eventually facilitate the search for it.

Zebrafish strains
The zebrafish transgenic line used in the study is actb2:Mmu.Arl13b-GFP (Borovina et al., 2010). Embryos were raised at 32˚C in the dark. For imaging, embryos were soaked in with Bodipy TR (Molecular Probes) for 60 min prior to the desired developmental stage and were subsequently embedded in 0.8% low melting point agarose (Sigma Aldrich) in Danieau solution. Embryos were imaged between 3-and 14-somite stages (SS).

2-photon excitation fluorescence (2PEF) microscopy
To image deep enough into the zebrafish embryo and capture the entire Kupffer's vesicle (KV), each live embryo (n = 28) was imaged using 2PEF microscopy with a TCP SP5 or SP8 direct microscope (Leica Inc.) at 930 nm wavelength (Chameleon Ultra laser, Coherent Inc.) using a low magnification high numerical aperture (NA) water immersion objective (Leica, 25x, 0.95 NA). We imaged the KV of embryos labeled with both Arl13b-GFP and BodipyTR between 3-and 14-SS: 100 Â 100 Â 50 mm 3 3D-stacks with 0.2 Â 0.2 Â 0.8 mm 3 voxel size and 2.4 ms pixel dwell time were typically acquired to maximize the scanning artefact allowing to properly reconstruct cilia orientation in 3D (Figure 1figure supplement 1E) as described in Supatto and Vermot (2011). The fluorescence signal was collected using hybrid internal detectors at 493-575 nm and 594-730 nm in order to discriminate the GFP signal labeling cilia from the signal labeling the KV cell surface. To uncover the orientation of the KV within the body axes, the midline was also imaged. We typically imaged a volume of 600 Â 600 Â150 mm 3 comprising the midline and the KV from top to bottom with a voxel size of 1.15 mm laterally and 5 mm axially.

3D-CiliaMap: quantitative 3D cilia feature mapping
We devised 3D-CiliaMap, a quantitative imaging strategy to visualize and quantify the 3D biophysical features of all endogenous cilia in the 50 to 80 cells constituting the KV in live zebrafish embryos from 3-to 14-SS (Figure 1-figure supplement 1). We used Imaris (Bitplane Inc., RRID:SCR_007370) and custom-made scripts in Matlab (The MathWorks Inc., RRID:SCR_001622) to perform image processing, registration, and analysis, and to extract the following features: KV size, shape and volume, cilia motility, number of cilia per KV, cilia spatial distribution, orientation of rotational axis and surface density (Figure 1-figure supplement 1 and Table 1). Cilia motility, position, and orientation in 3D, as well as the reference frame of the body axes, were obtained from 2PEF images and exported from Imaris to Matlab using ImarisXT (Figure 1-figure supplement 1D,E). Since our analysis relies on the fluorescence signal from the cilia, we discarded embryos with levels of GFP expression too low to analyze them. Similarly, a few cilia per vesicle could be discarded when the signal or the spatial resolution was too low to accurately determine motility (11%) or orientation (14%). Each cilium was defined as a unit vector from its base to its tip (Figure 1-figure supplement 1E). Cilia positions were registered in the body plan reference frame. To estimate the KV surface, we fitted an oblate spheroid to the distribution of cilia bases using the Ellipsoid fit Matlab script by Yury Petrov (Northeastern University, Boston, MA) (see fitting residues in Table 1). We used cilia vector components in the local orthogonal basis (ẽ f ;ẽ n ;ẽ m ) defined at each cilium position on the spheroid surface to quantify cilia orientation angles q (cilium tilt angle respective to the surface normal) and j (cilium orientation on the KV cell surface), as shown in ( Figure 1C). The experimental values were combined from different embryos and displayed in rosette histograms using Matlab. Finally, to estimate the local cilia density we transformed the spheroid into a sphere with surface density conservation and computed a spherical Voronoi diagram of cilia distribution based on the sphere_voronoi Matlab package by John Burkardt (Department of Scientific Computing, Florida State University, https:// people.sc.fsu.edu/~jburkardt/). The scripts are available upon request.

Flow calculation
When calculating the flow in a KV we first approximate it with a sphere of equal volume. We describe each cilium as a chain of 10 spheres (radius a=0.2mm) with a total length of L=6 mm, circling clockwise along a tilted cone with a frequency of 25 Hz. The phases were chosen randomly under the constraint that collisions between cilia were prevented (in rare cases, when inconsistencies in datasets led to unavoidable collisions, a randomly chosen cilium was removed). The mobility matrix M of the system was calculated using the Green's function for point forces inside a spherical cavity as described in Maul and Kim (1994). For the diagonal elements (self-mobility of a particle) we used the expressions where M rr denotes radial and M tt the tangential mobility. x ! i is the position of particle i relative to the center of the KV. In each step the forces on the particles representing points on cilia were calculated by solving the linear equation systemṽ i ¼ P j M ijFj and the fluid velocities subsequently from the Green's function.
As a simple and well-defined measure to characterize the intensity and directionality of the flow, we introduced the effective angular velocitỹ ÂṽðxÞ dV; that is, the angular velocity of a uniformly rotating sphere with the same angular momentum as the fluid in the KV (shown in Figure 4C). Note that the angular momentum is used solely as a velocity measure since the fluid inertia is negligible. can be calculated directly from the force distribution as thus omitting the need for spatial integration. This non-trivial expression can be derived from the following considerations. From symmetry arguments, we know that the effective angular velocity caused by a point forceF acting at pointx inside the cavity can only have the form with an unknown scalar function w r ð Þ. We now consider a distribution of forces on a concentric sphere with radius r i such that the velocity inside is ṽx ð Þ ¼ 0 Âx

1;
jxj r i jxj À3 ÀR À3 r À3 i ÀR À3 ; jxj > r i ( : Using definition (4), the effective angular velocity of this distribution can be obtained by spatial integration with the result At the same time, the force density at the inner sphere that maintains the velocity profile (7), is Inserting this force density into Equation (6) and integrating over the inner sphere gives The expressions (8) and (10) become equivalent when w r i ð Þ ¼ 3 16phR 3 1 À r 2 i R 2 , which leads to the Equation (5) for the effective angular velocity.

Randomly generated vesicles
In addition to the available datasets, we extended our analysis to randomly generated cilia distributions that shared the main features with those observed in real KV. The vesicles were assigned a radius R ¼ 35 m. We randomly distributed N c ¼ 70 cilia with a density function r that had its maximum at a; b ð Þ ¼ 0; 15 ð Þ and a pair correlation function gðgÞ resembling the measured one ( Figure 4C). All cilia were tilted meridionally (towards dorsal) with a tilt angle ¼ 60 Â cos b. The randomly generated distributions allowed us to study the reliability of the proposed flow sensing hypotheses without being limited by the number of vesicles analyzed experimentally. In each simulation with randomly generated vesicles, 200 vesicles were simulated in order to obtain stable results.

Torques
To test the mechanosensing mechanism, we chose three immotile cilia situated in the left region of the KV (up to 45˚away from the left pole), and three cilia situated in the right region of the KV (up to 45˚away from the right pole). In case the number of immotile cilia in one region was insufficient we randomly assigned additional cilia as immotile. To test the case with maximum sensitivity, the passive cilia were set normal to the surface ( ¼ 0). After determining the forces on all cilia, the meridional component (in the direction ofẽ m on the left and Àẽ m on the right) of the torque vector was evaluated around the base of a cilium. Figure 7B,C shows the cumulative distributions of instantaneous values of these torques, their temporal average, as well as the temporal and ensemble average for a group of 3 cilia on each side.

Particle diffusion
We evaluated the model based on diffusion of signaling particles with a Langevin-dynamics simulation in the fluid velocity field evaluated before (Video 3). We assumed that the particles are secreted from random points in a region 30˚around anterior and captured whenever they encounter a cilium elsewhere. The simulation step was 0.001 s and the number of particles traced 1000. The diffusion constant of a particle was determined as D ¼ k B T= 6ph r Stokes ð Þ with the fluid viscosity h ¼ 0:001 Pa s. Particles captured by cilia in the left and right region (up to 45˚away from the left/right direction) were counted and the average ratio N left =N leftþright , as well as its 5th and 95th percentile were plotted.

Statistical analyses
We did not compute or predict the number of samples necessary for statistical differences because the standard deviation of our study's population was not known before starting our analysis. Biological replicate corresponds to the analysis of different embryos of the same stage. Technical replicate corresponds to the analysis of the same embryo imaged the same way. The sample size (replicate and number) to use was as defined by our ability to generate our datasets. We routinely analyze 5 to 10 embryos at each considered stage. Both the mean and the SD (Figure 4-figure supplement 1 and Table 1) were calculated for several of the KV and cilia features measured. For analyses between two groups of embryos, differences were considered statistically significant when the p-value<0.05, as determined using a two-tailed and paired Student's t-test. Circular statistics (resultant vector length r and 95% confidence intervals on the estimation of the mean angle) where computed using the CircStat Matlab toolbox (Berens, 2009). Descriptive statistics (cilia density maps) were displayed using Matlab custom scripts (the scripts are available upon request).