In order to understand the mechanisms driving CSF flow in the central canal (CC), we first developed a thorough in vivo quantification of CSF flow in the CC. The obtained flow profiles will later serve as a basis to compare to our theoretical model. To establish quantitative and reproducible flow profiles along the dorsoventral (D-V) axis of the CC, we carried out an improved automated analysis inspired by our previous manual kymograph measurements (Cantaut-Belarif et al., 2018; Sternberg et al., 2018). By co-injecting TexasRed-Dextran (3,000 MW) with 20 nm-diameter green fluorescent beads in the diencephalic/mesencephalic ventricle (DV/MV), we measured both the volume of cavities and particle trajectories. The dextran rapidly diffuses throughout all fluid-filled cavities, including the CC (Figure 1A), enabling its visualization and assessment of the injection quality. Fluorescent nanobeads follow the bidirectional CSF flow (Figure 1B) but also undergo random Brownian motion due to their small size (Video 1). Bead trajectories also show the presence of a few regions of recirculation, commonly referred to as vortices (Figure 1B, Video 1). From the videos of fluorescent beads trajectories (Video 1, Figure 1 C1), we generated kymographs at various D-V positions in the CC (Figure 1 C2, C3, C4, and Video 2). The typical velocity profile is bimodal with positive flow (rostral to caudal) in the ventral side and negative flow (caudal to rostral) in the dorsal side, and reverses close to the center of the CC (relative D-V position = 0.53, Figure 1D from 110 wild type (WT) embryos each imaged in two rostro-caudal positions).

Video 1

Download asset

Video 2

Download asset

Figure 1

Download asset Open asset

Automated kymographs show better reproducibility compared to particle tracking velocimetry (PTV) or particle imaging velocimetry (PIV), as discussed in Materials and methods section. Velocity profiles exhibit a maximal speed of 4.78 ± 0.79 μm.s⁻¹ (relative D-V position = 0.29 on average) ventrally and a minimal speed of - 4.80 ± 0.82 μm.s⁻¹ (relative D-V position = 0.82 on average) dorsally (Figure 1B1 and B2). The profiles are nearly symmetrical (relative D-V position at the zero-crossing = 0.53 on average). Similarly, the average speed throughout the CC is null (- 0.02 ± 0.15 μm.s⁻¹). This implies that the average CSF flow rate through a section of CC is null as well. According to Poiseuille's relationship, this refutes the hypothesis of a global pressure gradient as the unique driver of the flow in the CC.

To build a theoretical model of CSF flow, essential parameters such as the canal geometry and the local cilia dynamics have to be quantified. The CC geometry, as estimated by TexasRed-Dextran (3,000 MW) injections, mostly corresponds to a cylinder (Figure 2A1 and A2) of diameter 8.9 ± 0.9 μm (129 WT embryos, Figure 2A3). The dynamics of motile cilia in the CC are investigated in 30 hpf Tg(β-actin:Arl13b-GFP) embryos in which GFP labels all cilia (Borovina et al., 2010; Cantaut-Belarif et al., 2018; Sternberg et al., 2018). Density and motility of cilia are higher in the ventral CC as previously reported (Borovina et al., 2010; Cantaut-Belarif et al., 2018; Sternberg et al., 2018). Movies recording the cilia beating illustrate the high density of cilia and diversity in ciliary length and beating frequency (Figure 2B, Video 3 and Figure 2—figure supplement 1). While some cilia beat at high frequency (up to45 Hz, see Materials and methods) with a single frequency pattern, others beat at lower frequency (~13 Hz) with more complex and irregular patterns. An automated quantification of cilia beating frequency using temporal Fourier transform (Figure 2B1) revealed distinct spots of constant frequencies (Figure 2B2). These regions likely correspond to the area traveled by a single cilium during one beating cycle. For each spot, we extracted the main frequency, orientation, length and height (Figures 2C1, 3, 4). We separated dorsal and ventral motile cilia based on their relative position to the center of CC (see Materials and methods). Motile ventral cilia were four times more numerous than the dorsal ones. Ventral motile cilia were typically 5.8 μm long and beat at 38.1 Hz with a caudal tilt of 62.6° towards the tail (median of the absolute value, see Materials and methods). This cilia length corresponds to a height of 2.7 μm in the CC, representing a bit less than half of the CC on the ventral side. We observed a caudal tilt, consistently with previous reports (Borovina et al., 2010; Cantaut-Belarif et al., 2018). The beating frequencies we measured were higher than the 12–15 Hz previously reported in zebrafish CC (Kramer-Zucker et al., 2005) (see Materials and methods for a possible explanation), but are consistent with beating frequencies of ependymal cilia located in the brain ventricles of rats (Smith et al., 2012) and zebrafish larvae (Olstad et al., 2018; Smith et al., 2012). Dorsal motile cilia were sparse and evenly distributed along the rostro-caudal axis: they beat with a slightly smaller central frequency (~32 Hz, Figure 2C2), but exhibited similar properties to ventral cilia otherwise (see Figure 2C3, 4, caudal polarization of 59.7 °, length of 5.5 μm).

Video 3

Download asset

Figure 2 with 1 supplement see all

Download asset Open asset

Figure 3 with 1 supplement see all

Download asset Open asset

Figure 4 with 1 supplement see all

Download asset Open asset

The CC is roughly shaped like a cylinder (Figure 2A) and contains motile cilia that line primarily the ventral wall (Figure 2B). The broad distribution of cilia beating frequency and cilia length, as well as the high density of cilia in the CC that could interact with each other and thereby introduce non linearities, would make it extremely challenging to model the effect of each individual cilium beating. As a first approximation, since the CSF bidirectional flow does not vary over time and over rostro-caudal axis (Figure 1 C4), we built a two-dimensional model that homogenizes cilia contribution as a constant force. The model divides the CC into two regions of equal thickness $h=d/2$ (Figure 3 A1, A2). It is assumed that the influence of dorsal cilia can be neglected, while active ventral cilia generate a constant volume force $f_v$ within the ventral region. Dimensional analysis gives ${\displaystyle f_v=\alpha \mu f/h}$, where $f$ is the cilia mean beating frequency (40 Hz), μ is the CSF dynamic viscosity and α is a numerical factor of order unity. The Stokes equation controls flow at small scales as inertial terms are negligible in comparison to viscous terms, and simplifies in dorsal and ventral regions as:

where $x$ and y represent the coordinates respectively in the rostro-caudal and dorso-ventral directions. $dP/dx$ is the common pressure gradient applying in both regions, and ${\partial }_y^2{v}_{dorsal}$ is the second y-derivative of the rostro-caudal velocity $v_{dorsal}\left(y\right)$. The derived flow obeys two parabolic profiles:

where $y_A$ and $y_B$ are two integration constants to be determined, and $dP/dx$ is the pressure gradient. Matching velocity and shear rate at the boundary between dorsal and ventral regions ($y=h$), two equations can link these three unknowns. To solve this system of equations, a third relation is provided by the null-flux condition across the diameter of the central canal as observed in vivo:

In the particular case of central canal geometry, the ventral ciliary region occupies about half of the diameter that is $h=d/2$. This enables to find: $y_A=d/4$, $y_B=3d/4$ and $dP/dx=f_v/2$, and the velocity profile adopts this symmetrical biparabolic form:

More details concerning the derivation of the velocity profile can be found in the Appendix 1. Note that this symmetry in the velocity profile, represented in Figure 3A3, is only obtained for the particular case $h=d/2$.

With a unique fitting parameter α that controls the amplitude of the constant volumetric force chosen as 0.5, (leading to $f_v=4000$ N/m³), we obtain a symmetric bimodal flow with extreme velocities equal to ±5 μm/s (Figure 3B), and an axial pressure gradient $dP/dx=2000$ N/m³. We additionally solved the 2D Stokes equation by Finite Element Method numerical simulations (Comsol Multiphysics software). Theoretical and numerical predictions are in good agreement with the 57 velocity profiles measured in vivo and tested here (Figure 3B).

The model is further refined to account for the spatial heterogeneity in the cilia distribution in the ventral region. We designed spots of constant volume force $f_v$ regularly distributed along the ventral wall of CC with a period $w$ (Figure 3C1), and added some regions without cilia in the ventral side. With this 'sparse cilia' model, numerical simulations show that recirculation regions appear between these spots (Figure 3C2). The amplitude of the recirculation regions is more pronounced as $w$ increases. For large passive regions of rostro-caudal extension $w$ larger than the diameter of the CC (${\displaystyle w>d}$), all the streamlines are located within the vortices, and no directed flow is observed between the active regions. This model thus predicts that vortices appear due to the alternation of active and passive regions. Note that the presence of passive regions also decreases the pressure in the CC proportionally to the fraction of passive regions (Figure 3—figure supplement 1A), a direct consequence of the linearity of the Stokes equation. We can further show that differences in beating frequency between active neighboring cilia lead to the emergence of the same recirculation regions (Figure 3—figure supplement 1D).

To verify that asymmetric cilia beating is solely responsible for the bidirectional flow in vivo as proposed by our model, we aimed at measuring the CSF flow in several mutant embryos for which either ciliogenesis or ciliary motility is affected; The traf3ip^tp499 mutation, referred to as elipsa (Omori et al., 2008), affects a protein involved in the intraflagellar transport (IFT), which leads to early cilia degenerescence at ~30 hpf. The dnaaf1^tm317b mutation, corresponding to the gene formally called lrrc50 (van Rooijen et al., 2008), leads to a complete loss of cilia motility by destabilizing the ciliary axoneme (van Rooijen et al., 2008). The cfap298^tm304 mutation, referred to as kurly (Brand et al., 1996), affects in a thermosensitive manner a cytoplasmic protein necessary for cilia motility and polarity (Jaffe et al., 2016). This protein is partly inactivated at the temperature at which our experiments were conducted (28°C, Jaffe et al., 2016). The foxj1a^nw2 newly-generated mutation (see Materials and methods and characterization in Figure 4 – figure supplement 1) reported here for the first time affects a transcription factor, which does not affect primary cilia but cilia whose tubulin is glutamylated (Figure 4—figure supplement 1E1, E2). Cilia with glutamylated-tubulin in the CC of WT embryos are caudally-polarized and ventrally-located and most likely correspond to motile cilia (Figure 4—figure supplement 1E1, E2).

In all these mutants with defective cilia, we injected beads in the DV/MV, but we could not observe the beads transported inside the CC. In anaesthetized mutant embryos, no particles entered the CC even several hours post injection, as illustrated here in elipsa embryos (Figure 4A1 and A2). This observation is consistent with previous reports (Cantaut-Belarif et al., 2018; Zhang et al., 2018) showing that long-range transport of beads down the central canal is altered in mutants with defective cilia. To account for this reduced transport, we verified whether the CC geometry was affected in embryos with defective cilia. Surprisingly, we observed that the CC diameter was dramatically reduced in all mutants tested: elipsa, lrrc50, foxj1a^nw2, and kurly (Figure 4B,C). Although the CC forms in these mutants, it is flattened along the D-V axis and not as straight as control siblings (Figure 4B1, B2 and C). In order to evaluate whether these features could be due to the convex body curvature of the embryos, we analyzed the CC structure of mutant embryos scospondin^icm15/icm15, which are also curled down but exhibit normal cilia activity (Cantaut-Belarif et al., 2018). For these mutants, we did not measure any difference in the CC diameter. This observation suggests that cilia integrity and motility are critical to properly form the CC and instruct its diameter and straightness.

As exogenous beads are not transported down the CC in anaesthetized mutants with defective cilia, we could not measure the local CSF flow, and hence test whether the cilia beating is directly responsible for CSF bidirectional flow. However, we noticed that allowing mutant embryos to spontaneously twitch (by not anesthetizing them before imaging) facilitated the transport of beads down the CC (see next section for the characterization of the effect of muscle contraction on CSF flow). By using the effect of muscle contractions, we were able to get beads down the CC in mutants with a moderate reduction of the CC diameter, such as elipsa and kurly. Fluorescent beads in the CC of elipsa and kurly mostly undergo random Brownian motion, but some allow to measure a residual flow with many recirculation zones (Video 4), in good agreement with the predictions of the sparse cilia model.

Video 4

Download asset

It is also predicted that the bidirectional CSF flow in the CC is generated locally and does not rely on a global pressure gradient generated in the ventricles. This is confirmed by multiple observations: i) even after large regions without cilia and local flow, the CSF flow in more caudal regions can show a bidirectional flow, as long as motile cilia beat locally in the ventral region (as predicted by Figure 3C). ii) When a portion of the CC of around 200 μm in length is isolated by using photoablation to block flow on both extremities (see Materials and methods and Figure 4D1), the canal kept its diameter at the center and the CSF still flowed bidirectionally (Video 5) with a similar velocity profile (Figure 4D2).

Video 5

Download asset

So far, we mainly restricted our study to anaesthetized zebrafish embryos. However, muscle contractions can also generate several biological flows (Blatter et al., 2016). Recent results demonstrated that muscle contraction can create a massive transient flow in the brain ventricles of zebrafish embryos (Olstad et al., 2018). Consequently, we investigated CSF flow and transport on embryos without anesthesia, which have spontaneously contractions of the body wall muscles (spontaneous twitches). Animals were restrained in agar and contractions were identified from the notochord displacement (white arrow) on the transmitted image (Figure 5A). Following each contraction, rapid displacements of CSF (too fast to be properly characterized) were observed in either the caudal or rostral direction, followed by a slower counter-flow in the opposite direction, whose velocity decayed until the original bidirectional flow is recovered (Figure 5B and C, Video 6, Video 7). The average CSF flow throughout the entire CC section is correlated with the instantaneous muscle contractions (Figure 5D) and showed a sudden increase after contractions. By comparing the increase in CSF flow over 5 s after each contraction with the contraction strength, we observed a robust flow rate increase after each contraction (Figure 5E; 204 contractions from 22 WT embryos). However, we did not find any correlation between the flow direction and estimates of the contraction direction or strength. Body contractions also enhance transport of beads in mutants with defective cilia, as mentioned in the previous section (Video 4). As no significant difference in the CC diameter was detected during the muscle contractions, we propose that the reported enhancement of CSF flow is not due to a local variation of the CC diameter, but instead due to a change in the volume of the brain ventricles that occurs during muscle contraction at larval stage (Olstad et al., 2018).

Video 6

Download asset

Video 7

Download asset

Figure 5

Download asset Open asset

We characterized above the local CSF flow in the CC, but did not investigate its ability to transport particles over long distances, despite its importance for embryonic development (Chang et al., 2016). The rest of the study will focus on long-range circulation of CSF, first in the CC, and then in the entire CSF circulatory system. At the micrometric scale of the central canal, Brownian diffusion can be sufficient to spread molecules without requiring a flow advection. However, larger particles (of diameter larger than about 10 nm) are less sensitive to Brownian motion and the presence of a flow can dramatically increase their transport efficiency throughout the embryo body. It is at the crossover between these two qualitatively different transport phenomena that the influence of the bidirectional flow in the CC will be characterized.

Here, we first demonstrate experimentally that CSF long-range transport of nanoparticles has characteristics of both passive diffusion and active transport. The particle transport follows an apparent diffusive law, as the propagation front of nanoparticles progresses with the square root of time (Zhang et al., 2018) and this propagation is found to be slower for larger particle size (Figure 6A1 and A2). On the other hand, the propagation front of 20 nm-diameter fluorescent beads is accelerated by muscle contractions, when the local CSF flow velocity is increased (Figure 6A2 and A3). An estimation of orders of magnitude also suggests that CSF long-range transport originates from a combination of pure diffusion and flow advection. Indeed, 20 nm-diameter fluorescent beads injected in the DV/MV are found to travel through the entire length of CC (~2 mm) in about two hours (Figure 6A2), far from the timescales expected for passive diffusion (over one day) or pure ballistic transport at the speed of the flow (with V = 5 μm/s, this would lead to less than 10 min).

Figure 6

Download asset Open asset

To account for these observations, we developed theoretical and numerical calculations according to CC geometry. We modeled the transport of particles by including an advection term in the diffusion equation. In a bi-dimensional canal (Figure 3A2), the transport equation can be written as:

where $c$ is the local concentration of particles, and ${\partial }_t$,${\partial }_x$ and ${\partial }_y$ are respectively the time and spatial derivatives in the rostro-caudal ($x$) and dorso-ventral ($y$) directions. $v\left(y\right)$ is the velocity profile (Equations 6, 7), and $D$ is the Brownian diffusion coefficient of particles in CSF. For a spherical particle of radius $r$, the Stokes-Einstein relation gives:

 where $T$ is the temperature and $k_B$ is Boltzmann’s constant. As described by Taylor (1953) and Aris (1956), the presence of velocity gradients increases the effective diffusivity of solutes due to a coupling between classical Brownian diffusion and flow advection. Over long distances, transport still follows a diffusive process, but with an effective diffusivity enhanced by the shear flow (Aris, 1956; Taylor, 1953). From Equation 8, we performed the theoretical derivation of this effective diffusivity $D_{eff}$ for particles of arbitrary size in the presence of a bidirectional flow (in Appendix 2, and Appendix 2—figure 1), based on the derivation done by Bruus for a unidirectional Poiseuille flow (Bruus, 2008). This leads to:

where $\mathrm{P}\mathrm{é}=Vd/D$ is the dimensionless Péclet number, comparing the characteristic times for spanwise diffusion $d^2/D$ and convection $\mathrm{d}/V$. The particle transport process in the CC thus follows a diffusion law whose diffusion coefficient is increased with the flow maximal velocity $V$. Note that the prefactor is larger here than for a unidirectional flow of similar velocity (Appendix 2).

This theoretical prediction is confirmed by numerical simulations, comparing the long-range transport of particles of diameter 40 nm (in the range of commonly observed embryonic extracellular secreted vesicles in the CSF [Bachy et al., 2008], without flow, and in the presence of a bidirectional flow (Figure 6B) of variable velocity. The simulations show that at times longer than the characteristic diffusion time $d^2/D$ = 10 s, a diffusive-like gaussian concentration profile is recovered with a wider extent in the presence of a flow (Figure 6B). This diffusive transport is numerically obtained for three different maximal speeds of the flow: 0, 6 μm/s, and 20 μm/s (Figure 6C1), and shows that this effective diffusion is enhanced by the flow speed, as predicted by our theoretical derivation. Another interesting feature of this diffusivity enhancement is its weaker dependence on the size of the transported particles, as opposed to classical diffusion for which the diffusivity is inversely proportional to the particle size (Figure 6C2, Equation 7).

In vivo front tracking experiments illustrate this flow-induced accelerated transport of 20 nm nanoparticles in paralyzed and contracting zebrafish embryos (Figure 6D). While both contracting and paralyzed embryos exhibit a diffusive-like transport, the effective diffusivity is higher in the case of contracting embryos. It confirms that muscle contractions can increase the long-range transport efficiency (an increase of the steady bidirectional velocity would have similar effects, Figure 6C1). Interestingly, the transport is predicted to be enhanced in both rostro-caudal and caudo-rostral directions, which is confirmed by in vivo experiments. When injected locally in the caudal part of the CC, fluorescent nanobeads are indeed transported rostrally towards the brain ventricles (Figure 6E1). Similarly to rostral-to-caudal transport, the reverse transport is sped up incontracting embryos (Figure 6E2).

To understand how the CSF carries particles between CC and the brain ventricles at early stages, we imaged all the CSF-filled cavities by injecting the small dye Texas Red-Dextran (3,000 MW) in live 30 hpf embryos (Figure 7). At this stage, cilia in the brain ventricles are mostly non-motile (Olstad et al., 2018). At this stage, CSF flow in the brain ventricles is mostly Brownian, although the heart beat is responsible for a pulsatile component to the flow at the heart beat frequency (Video 8), as previously observed (Fame et al., 2016; Olstad et al., 2018). Our observations indicate that before 30 hpf CSF flow in the brain ventricles does not contribute to the overall CSF circulation. After the brain ventricles, we investigated how CSF flows in between the brain ventricles and the CC. During embryogenesis in rats (Sevc et al., 2009) and zebrafish (Kondrychyn et al., 2013, Ribeiro et al., 2017; Guo et al., 2018), the progressive closing of the neural plate gives rise to a ‘primitive lumen’(pL) that changes shape and volume over time. At 30 hpf in zebrafish embryos, we investigated here three transient paths where CSF flows in between the brain ventricles and the CC (Figure 7). First, the CSF circulates in a small ventral opening of the RV pL, in continuity with the ventral larger opening of the the CC pL. (as shown in histology in Fame et al., 2016). We measured that the diameter of this ventral opening is smaller than the CC diameter (7.1 ± 0.7 μm) and is 233 ± 24 μm long (Figure 7A and B – green arrow, 7C,11 embryos, Video 9) and that few motile cilia enable there a unidirectional flow from rostral to caudal, that is towards the CC (Video 10). This fast and unidirectional flow follows a parabolic Poiseuille-like velocity profile with a maximum speed of 15 μm / s (Figure 7D1 and D2). We also observed that CSF circulates in a ventrally-located small opening of the DV/MV pL connected with the TV (Figure 7B- red arrow, Video 9). Finally, we observed a small dorsal opening of the CC pL with a very narrow diameter (~2.5 μm, Figure 7A and E –blue arrow, Video 11). The CC pL is indeed composed of two main openings, with a very narrow but continuous connection between the two (Kondrychyn et al., 2013, Ribeiro et al., 2017). We observed that the apical domains of the two openings of the CC pL, marked with the tight junctions protein ZO-1, are closely-spaced, so that CSF flow is independent in each opening (Figure 7E). This is the reason why we neglected the influence of this dorsal opening of the CC pL to model the CSF flow in the larger opening, and we called the CC only its larger ventral opening. Interestingly, we found that the CC ventral and dorsal openings converge into a single opening, dense in motile cilia, in the caudal part of the spinal cord (Figure 7 E2, Video 12). During development, as it has already been partly described in rodents and zebrafish (Sevc et al., 2009; Kondrychyn et al., 2013, Ribeiro et al., 2017), dorsal and ventral openings of the CC pL converge ventrally at three dpf (Figure 7F) to finally merge into a single central canal lumen at five dpf (Figure 7G). Further histological investigations will be needed to clarify the geometry of all cavities where CSF flows in the embryo.

Video 8

Download asset

Video 9

Download asset

Video 10

Download asset

Video 11

Download asset

Video 12

Download asset

Figure 7

Download asset Open asset

In previous sections, we showed that long-range transport of particles in the CSF is accelerated by a fast and unidirectional flow between the ventricles and the CC, and by a slower but efficient bidirectional flow in the CC. We now investigate whether particle transport between the brain and the spinal cord is important for the development, or whether the ventricular and spinal CSF systems are independent at this stage. We performed 2-photon mediated photo-ablations at the beginning of the CC to disconnect the CC from the brain ventricles (Figure 8A1 and A2). We found that 30 hr after ablation, half of the ablated larvae were shorter in length than siblings that underwent a control ablation of the yolk extension or had no ablation (Figure 8B1 and B2). Larvae with the opening ablated had an increase in body-axis curvature (Figure 8C1 and C2). Similar effects on development were observed when CSF from the brain ventricles was drained (Chang et al., 2016), which affected the transport of CSF, including RBP4, a transport protein involved in the synthesis of Vitamin A. These experiments indicate that global CSF circulation contributes to embryonic growth, and suggest that the transport of particles between brain and spinal cord is a critical process for development.

Figure 8

Download asset Open asset

In order to find particles that could be efficiently transported by the bidirectional flow, we investigated the content of the CSF and identified abundant extracellular vesicles using the label-free imaging technique, full field optical coherence tomography (Figure 8D and Video 13). Numerous vesicles with high lipids or proteins content filled the CC, in agreement with previous observations (Sternberg et al., 2018). Fluorescent labeling of CD63 (Video 14, Materials and methods) and electron microscopy with immunogold labeling of exosomes suggested that some of these particles may correspond to exosomes in the CSF (Figure 8E, Materials and methods) (Verweij et al., 2019). Electron microscopy also reveals a variety of other extracellular vesicles, such as lipoparticles, microvesicles. The diameter of these vesicles ranges between 50 and 250 nm (Figure 8E), typically in a regime where Taylor-Aris dispersion is important. The bidirectional flow of endogenous extracellular vesicles can be put in evidence in the CC using full field optical coherence tomography as well as fluorescence live imaging (Video 13, Video 14). Altogether, our works suggests that the effective CSF circulation contributes to embryonic growth via the transport of extracellular vesicles between brain and spinal cord.

Video 13

Download asset

Video 14

Download asset

All procedures were performed on zebrafish embryos and larvae before 5 days post fertilization in accordance with the European Communities Council Directive (2010/63/EU) and French law (87/848) and approved by the Brain and Spine Institute (Institut du Cerveau et de la Moelle épinière). All experiments were performed on Danio rerio embryos of the AB strain, and were mitfa^-/- (nacre) background. Animals were raised at 28.5°C under a 14/10 light/dark cycle until the start of the experiment.

All values are represented as median ± 0.5 standard deviation (s.d.) to evaluate the spread of each distribution, except for the error bars in the flow profiles showing the standard error of the mean (SEM), in which case we considered the average flow at one position. The speed distribution of particles depends on experimental errors and Brownian motion, and is therefore not considered. All statistics were performed using MATLAB. Statistical details related to sample size and p values, are reported in the figure legends. In figure panels, asterisks denote the statistical significance calculated by a two-tailed t test for samples with unequal variance (Welch t test): ^∗, p<0.05; ^∗∗, p<0.01; ^∗∗∗, p<0.001; ns, p>0.05.

30 hpf embryos were manually dechorionated, mounted in 1.5% low-melting point agarose in the lateral position. A few nanoliters of a mix containing 20 nm beads, 3,000 MW TexasRed dextran, and alpha-bungarotoxin (TOCRIS) were injected in the center of the diencephalic ventricle in order to label CSF and anaesthetize the fish in a single injection. The injections were performed using a Picospitzer device (World Precision Instruments). 20 nm carboxylate FluoSpheres of center wavelength 505/515 nm (yellow/green, F8888, Molecular Probes) were diluted to a 2% concentration in artificial CSF (containing in mM: 134 NaCl, 2.9 KCl, 1.2 MgCl2, 10 HEPES, 10 glucose, 2 CaCl2; 290 mOsM + /- 3 mOsm, pH adjusted to 7.8 with NaOH), and then sonicated for 10 s. 3,000 MW Texas Red dextran (Excitation/Emission wavelength: 595/615/ThermoFisher Scientific D3329) was diluted at 0.2% mass/volume concentration and was used both to serve as a marker of injection, and to assess the boundaries of the ventricles and central canal. A final concentration of 100 μM of alpha-bungarotoxin was used. We could visually assess that such concentration, as well as the ventricle injection of alpha-bungarotoxin, although not conventional, were sufficient to paralyze the embryos. For experiments in contracting fish, the alpha-bungarotoxin was replaced by artificial CSF. Following the injection, we waited one hour before performing experiments to allow the particles to be transported down the central canal. Beads larger than 100 nm-diameter were restrained to brain ventricles, while beads of 20 and 45 nm diameter could enter the central canal.

For CSF flow experiments, time-lapse images were acquired at room temperature on an inverted Leica DMI8 spinning disk confocal microscope equipped with a Hamamatsu Orca Flash 4.0 camera, using a 40X water immersion objective (N.A. = 0.8, pixel size 189 nm). 300 images of the beads mostly in the rostral part of central canal were acquired at a frame rate of 10 Hz for 30 s using Metamorph software (www.moleculardevices.com). Some experiments in contracting fish were performed at a frame rate of 40 Hz for 20 s to more accurately follow the muscle contraction.

A MATLAB Particle Tracking Velocimetry (PTV) algorithm was adapted from a classical PTV algorithm (Crocker and Grier, 1996) and used to obtain instantaneous particle velocities. This algorithm was only used to produce Figure 1B as it only worked consistently in the best SNR conditions. The first crucial step was to segment the particles inside the region of interest corresponding to the central canal. This region was automatically inferred based on an intensity threshold on the low pass filtered mean image. Particle segmentation was more efficiently performed using a machine learning-based segmentation software (iLastik Sommer et al., 2011). Particle fine positioning was performed by fitting a 2D gaussian to each segmented region (particle) and by extracting the center position. A large matrix of particle position (X, Y, time) was obtained at this point. Then, the PTV algorithm itself was applied. The algorithm connects the particles found at successive times in order to minimize the total squared distance traveled by all beads. Note that this algorithm is optimized for purely Brownian displacements and limits the number of particles that one can track before saturating its computer memory, as the number of calculation scales as N²_particles. Once the particles were tracked, each trace was obtained and kept only if it contained more than five time points. The instantaneous speed was calculated and plotted as a vector on the mean intensity image. An additional color code ranged from blue to red, corresponding to the minimal and maximal speed.

First, time lapses were automatically rotated and cropped with an extra manual step for fine adjustments to isolate a portion of the central canal and place it in a correct dorso-ventral orientation. A 2D wavelet filter (decomposition of level six with Haar wavelet coupled to a 2-D denoising scheme based on level-dependent thresholds penalizing low signal pixels) was applied to increase the local signal-to-noise ratio of the fluorescent beads. Stacked images were then re-sliced and particle trajectories were automatically tracked and analyzed for every successive dorso-ventral position. The average signal on three consecutive positions was calculated with a moving-window average. Each column was then divided by its average value to account for spatially or temporally inhomogeneous light illumination and photobleaching if there was any. Only pixels above the average signal were kept, and an internal Matlab segmentation algorithm (regionprops) was applied to extract regions of interest in each kymograph. In order to filter only the straight lines corresponding to particle motion, a few additional empirical filtering conditions were applied: Only regions with at least 15 pixels (to avoid random fluctuations), and corresponding to lines with eccentricity above 0.9 (to prevent from selecting big aggregates or any large objects) were kept. Additionally, to prevent selecting particles stuck in the canal (vertical line in the kymograph), or camera/illumination fluctuations (global change of intensity in the entire image, leading to a horizontal line in the kymograph), we filtered any region whose orientation was close to 0 or 90°± 180° (excluding regions whose cosine or sine of orientation were below 0.1). Note that the intensity threshold was not strong, but most of the particle filtering was made through these four last conditions. A few regions corresponding to individual particle rostro-caudal trajectory were thus obtained. The tangent of these line orientations was calculated, and corrected by pixel size and frame time, to obtain the individual particle velocity. All particles from one kymograph were aggregated to the average velocity and standard error to infer the local flow. Note that we refer to the number of events, or number of lines tracked, rather than the number of particles, as this algorithm cannot guarantee that one particle is not tracked several times if its dorsoventral position changes over time. The full MATLAB script can be found on Github. (Thouvenin, 2019; https://github.com/wyartlab/eLife_2019_OriginAndRole; copy archived at https://github.com/elifesciences-publications/eLife_2019_OriginAndRole).

We extracted some parameters from the flow profiles to compare the profile from many embryos. The maximal and minimal values, and their respective positions were extracted. The position at the center of the profile where the speed shifts from a positive to a negative value was also extracted. The position in the canal was normalized between 0 and 1, where 0 was the most ventral position where at least five moving particles were found and velocity exceeded 0.4, and one the most dorsal position where five moving particles were found and the velocity was not higher than −0.4.

In Figure 5D, the average flow versus time was plotted in embryos able to contract in agarose. The experimental protocol was similar as for the flow profile analysis, except that alpha-bungarotoxin was removed from the injection mix. The same automatic kymograph approach was used to obtain events corresponding to motion of the beads. For each event, each trace was kept and its speed, its central position in D-V and R-C axis, and its time of occurrence were extracted. If the speed profile was obtained by averaging all pixels over time along the D-V axis, the average flow versus time was obtained by averaging all particles from all D-V positions at a given time. When bidirectional flow was present (before contraction, or in an anaesthetized embryo), the average flow was close to 0. in order for the flow measurement to be meaningful, and for the Brownian motion to be averaged out, several events should be measured at each timepoint. It gives a constraint on the actual timescale we can probe for the average flow rate. In Figure 5D, the fish contraction strength was estimated to be correlated to the temporary increase in volumetric flow. To estimate such contraction, the average signal from all pixels within central was calculated and the absolute value of their time derivative was calculated (magenta curve in Figure 5D). An arbitrary threshold of four was empirically chosen to select points during contraction. The contraction strength used in Figure 5E was calculated as the sum of the displacement signal during one event with displacements values above threshold. This value mixes the instantaneous derivative value and the contraction duration, as they can both play a role in the establishment of different flows. In datasets where embryos did not move, the average flow and displacement were calculated over the entire time series. The average displacement value was calculated in order to avoid getting a large value due to the large time span of the dataset and the use of absolute value.

Besides the PTV and kymograph approaches, we also tried to use a particle imaging velocimetry (PIV) algorithm, which is generally used to measure local flow (Raffel et al., 2007; Santiago et al., 1998). This section aims at comparing all these techniques qualitatively with respect to our imaging conditions. Before describing why we mainly used the automatic kymograph approach, we have to recall that we could not inject the particles directly in the central canal, and exogenous particles above 100 nm in diameter stayed confined inside the brain ventricles. Therefore, we could only image particles with diameter smaller than 100 nm, and mainly used 20 nm particles with fairly low signal, especially at high imaging frequencies. The Brownian motion of 20 nm beads is considerable and can generate instantaneous displacements of the same order of magnitude as the CSF flow itself. With that in mind, PIV should be efficient to measure the local flow even in low SNR conditions, as it uses a cross-correlation between successive frames. Nevertheless, because of the importance of Brownian motion, the cross-correlation fails at producing a peak of high SNR. Besides, because of the geometry of the canal and the bidirectional nature of CSF flow, we could not use a large square region-of-interest to calculate the cross-correlation, as typically used in PIV. As particles can move in opposite directions in a 6 × 6 μm² (32 × 32 pixels) square, the cross-correlation results in several small competing peaks. Asymmetric region of interests should be more adapted to take full advantage of the steady flow along the rostro-caudal direction, but performing a PIV with asymmetric regions is not standard and not trivial. On the other hand, PTV tracks every individual particle and should produce the trajectory of each particle versus time. One can expect that a given trajectory would then have a Brownian component plus a directed component, each of them being separable. Unfortunately, PTV was not very robust with low SNR, and not very efficient when particle density was quite high. The most critical part was to efficiently detect the particles without error. A first solution is to have a high threshold that selects only the brightest particles, as what was done in Figure 1. Nevertheless, only a few particles could be tracked for a given dataset, and the profile could not be inferred from only a few particle positions. If the threshold were lowered, it would increase the calculation time, and the probability of false detection. Since the PTV algorithm tends to minimize the distance traveled by particles, as the number of false detection increases, the false associations increased even more leading to an important underestimation of particle speed. In contrast, the automatic kymograph seemed more robust in this context because the flow is mostly steady in the rostro-caudal direction, allowing us to perform efficient averaging in this direction. Brownian motion was averaged out by averaging on 3–5 consecutive dorso-ventral positions, so that the particles produced traces long enough to be identified. The main strength of the kymograph approach was that many particles can be detected at every position, so that even if several mistakes were made, the standard error stayed small. The longer the total acquisition time, the more particles could be averaged, and the more precise the technique was, which was not true for both PTV and PIV in which all images were independent.

Nonetheless, most of our efforts focused on the automatic kymograph approach, and that probably PTV and PIV algorithms may be adapted to our configuration. We could not find any available PTV and PIV algorithm that could work with our data, and theories behind both techniques become more complicated in non-optimal configurations. In contrast, the kymograph just required to fit lines in an image, for which many solutions were available and just required small adjustments.

A Z-Stack of 30 images around the central canal of 30 hpf zebrafish embryos filled with Texas Red was acquired. A custom MATLAB program (Thouvenin, 2019, Github, https://github.com/wyartlab/eLife_2019_OriginAndRole) was used to open, rotate, and crop the image around central canal, then to calculate the maximum intensity projection. The rotation was crucial to fit the central canal projection by a manually selected rectangle, whose width measured the canal average diameter. Even though the canal diameter was not constant along the rostro-caudal direction, the rectangle was placed in order to minimize the total error. The canal diameter was averaged on large rostro-caudal portions, to provide fairly good estimate of its size. Data from wild type and mutant siblings were randomly mixed to perform blind analysis. For Kurly embryos displayed in Figure 4, Texas Red was not injected, but we used the time average image of fluorescent beads filling the central canal. The results on WT siblings showed remarkably similar measurements as those obtained with Texas Red. Thus, we decided to combine all these datasets.

30 hpf embryos from Tg(β-actin:Arl13b-GFP; scospondin^icm13/+) incrosses were manually dechorionated, laterally mounted in 1.5% low-melting point agarose, and paralyzed by injecting 1–2 nL of 500 µM alpha-bungarotoxin (TOCRIS) in the caudal muscles of the trunk. A Leica spinning disk confocal microscope equipped with a 40X water immersion objective (N.A. = 0.8) was used to measure cilia dynamics at a 100 frames per second for 3 s. Note that our estimation of beating frequency were limited to 50 Hz. Because we found that the average frequency of cilia beating in the central canal was higher than the value of 12–15 Hz previously reported (Kramer-Zucker et al., 2005), we imaged cilia dynamics with three different microscopes (two commercial spinning disk microscopes from 3i and Leica, and one 2-photon scanning microscope using the scanning artifact, as previously described by Supatto and Vermot (2011), in order to validate our results). However, all the results displayed in this report were obtained with the same Leica commercial microscope to ensure data homogeneity. We can only postulate that the difference found in beating frequency reflects a bias for long cilia with possible slower kinematics when using differential interference contrast (DIC) to identify cilia.

Datasets of cilia position over time were analyzed using custom MATLAB (Thouvenin, 2019, Github, https://github.com/wyartlab/eLife_2019_OriginAndRole) scripts based on Fourier transform. First, we opened, rotated, and cropped the image around the central canal. A straight line was manually drawn from the mean intensity image to define the center of the central canal. We also applied a 4 × 4 spatial sliding-average to increase the local signal-to-noise ratio. The temporal Fourier transform of each averaged pixel was calculated, and the position of its norm extracted and transformed into a frequency value. We could reconstruct 2D maps with the main frequency for each pixel (Figure 2B2). It drew distinct regions of interest of constant frequency, which we assumed to be the regions individual cilia could reach within one beating cycle. Within each region, the calculated main frequency can vary a few Hz, but was grouped around a central frequency with a narrow bandwidth. Successive mask images of pixels with main frequency within increasing bandwidth of 5 Hz (5–10 Hz, then 10–15 Hz, and so on until 50 Hz) were created, and analyzed using MATLAB image processing functions (regionprops). We selected as cilia only regions of constant frequency with area larger than 200 pixels (at least 3 μm in diameter). The position of each region was then compared to the position of the central line to label the cilia as dorsal or ventral cilia. First, the central frequency within each region was calculated. Because we extracted the average frequency over a large spatial scale, this measurement was independent of the cilium phase, and potentially free from sampling artifacts. Then, each region was associated to an ellipse, whose long axis was assumed to be the cilia length. After geometrical calculation, the ellipse orientation gave the angle between the cilia and the dorso-ventral axis. The values corresponding to the negative orientations may principally originate from the few beating dorsal cilia, also oriented towards the caudal end, but from dorsal to ventral (leading to a negative angle). To account for this, we calculated the median of the absolute value of the cilia orientation in the results section. Note that it gave the same value as the median of only the positive tilted cilia. Finally, we calculated a parameter called beating height, as the cilia length multiplied by the cosine of orientation (independent of the sign of the orientation), which gave an estimate of which portion of central canal was occupied by beating cilia. Although some reports (Smith et al., 2012) used the power spectral density to quantify the main beating frequency, we did not find a difference between the main frequency found with the Fourier transform versus the power spectral density calculation.

foxj1a mutants were generated with CRISPR/Cas9 mediated mutagenesis as previously described Hwang et al. (2013); Olstad et al. (2018). The gRNA sequence for Foxj1a gene (GGACGTGTGCTGTCCTGTGC) was identified using ZiFiT Targeter website (http://zifit.partners.org/ZiFiT_Cas9) from the genomic DNA sequence obtained from Ensembl. Annealed oligos were cloned into the pDR274 plasmid (Addgene Plasmid #42250). The sgRNA was in vitro transcribed using MAXIscript T7 kit (Invitrogen). Cas9 mRNA was in vitro transcribed from pMLM3613 plasmid (Addgene Plasmid #42251) using mMessage mMachine T7 kit (Invitrogen) and polyA tailed with a polyA tailing reaction kit (Invitrogen). A mixture of 25 pg of sgRNA and 600 pg of Cas9 mRNA were injected in embryos at one-cell stage. CRISPR/Cas9-induced mutations and germ line transmission of the mutations were detected by T7EI assay and characterized by sequencing. Three separate alleles were recovered from one founder. All mutations resulted in a frame shift from amino acid 66 and an early stop codon, which was before the forkhead domain. All homozygous mutants for the allele exhibited a bent body axis. This study used the foxj1a^nw2 allele, which carries a 5 bp deletion. See more information on Figure 4—figure supplement 1. Heterozygous adults and mutant larvae were genotyped using a KASP assay (LGC genomics) on a qPCR machine (ABI StepOne). The generation of the Foxj1a mutant line and fin clipping procedure on adults were approved by the Norwegian Food Safety Authorities (FOTS applications 16425 and 12395).

For immunostaining, straight body axis and curved body axis embryos resulting from an incross of heterozygous adults were selected and processed separately. The genotype of the clutch was confirmed by KASP genotyping of siblings presenting a curved body axis. Dechorionated and euthanized embryos (collected between 26 and 30 hpf) were fixed in a solution containing 4% paraformaldehyde solution and 1% DMSO for at least 1 hr at room temperature. Embryos were washed with 0.3% PBSTx (3 × 5 min), permeabilized with acetone (100% acetone, 10 min incubation at −20°C), washed with 0.3% PBSTx (3 × 10 min) and blocked in 0.1% BSA/0.3% PBSTx for 2 hr. Embryos were incubated with either glutamylated tubulin (GT335, 1:400, Adipogen) or acetylated tubulin (6-11B-1, 1:1000, Sigma) [MOU1] overnight at 4°C. The next day samples were washed (0.3% PBSTx, 3 × 1 hr) and incubated with the secondary antibody (Alexa-labeled GAM555 plus, Thermo Scientific, 1:1,000) overnight at 4°C. Next, the samples were incubated with 0.1% DAPI in 0.3% PBSTx (Life Technology, 2 hr), washed (0.3% PBSTx, 3 × 1 hr) and transferred to a series of increasing glycerol concentrations (25% and 50%). Stained larvae were stored in 50% glycerol at 4°C and imaged using a Zeiss Examiner Z1 confocal microscope with a 20x plan NA 0.8 objective.

We intentionally damaged the central canal using different levels of light excitation in a few experiments. All experiments were performed after Texas Red injection. In experiments shown in Figure 4D1, a continuous medium intensity (25 mW) of green (565 nm) light was applied to the entire field of view 300 μm on the left and later 300 μm on the right of the displayed field of view for 5 min. After a few minutes of illumination, the central canal collapses, similarly to what was found in ciliary mutants (Figure 4B2). Debris were also seen, as shown on the sides in Figure 4D1. We never observed a breach in the central canal or fluorophores leaking out of it. A possible hypothesis is that TexasRed photobleaching due to the illumination may create damage to the central canal, cells surrounding it, or their motile cilia leading to an obstruction of the canal. In Figure 8A1 and A2, a light pulse of 200 mW (at a wavelength of 810 nm on a two-photon microscope) was scanned 200 times on 50 pixels lined along the dorso-ventral axis, right after the rostral extremity of the CC, or at the beginning of the yolk sac extension for control siblings. This reproducibly created a hole, and hence a breach in the central canal or in the yolk sac extension tissue, around which one can observe fluorescence leaking out in the surrounding tissue.

To measure CSF transport, four groups of eight 30 hpf embryos (16 contracting and 16 paralyzed embryos) were injected with beads and TexasRed within a five minutes time frame. We recorded the exact injection minute for each embryo. Widefield fluorescence images of each embryo were acquired every 10 min for 60 min using an epifluorescence microscope (Nikon AZ100) equipped with a 2X air objective imaged the entire embryo in one image. The image creation time recorded by the computer was automatically extracted and subtracted from the injection time. The front propagation was obtained by applying a low-pass filter on all images and by manually drawing the contour from the rostral extremity of the CC to the propagation front with a MATLAB script. The contour was drawn by successive lines changing direction every time the canal changed direction. The total length was calculated as the sum of all individual length. As for other manual arbitrary measurements, a random permutation was applied to the different images to reduce bias.

We imaged CSF-filled cavitites with a 2-photon microscope (Vivo Multiphoton, 3i, Intelligent Imaging). Live embryos were injected with Texas Red and alpha-bungarotoxin in the diencephalic/mesencephalic ventricle. Other embryos were fixed in a solution containing 4% paraformaldehyde solution and 1% DMSO for at least 1 hr at room temperature. ZO-1 and DAPI immunostaining was then performed following a similar protocol as previously described. Z-stacks of between 30 and 300 μm (depending on the location) around either the brain ventricles or the central canal were obtained and displayed as a maximal projection to visualize the different lumens. Video 9 was obtained by calculating the mean intensity every eight planes. Similar protocol was followed for 3 dpf and five dpf zebrafish larvae.

30 hpf embryos were injected with Texas Red and alpha-bungarotoxin in the DV/MV and imaged with a confocal microscope (SP8, Leica). A Z-stack of 50 μm around the position of the ventral opening of the RV pL was obtained and displayed as a maximal projection to visualize its entire length as well as the beginning of central canal. A MATLAB script was used to perform measurement of this geometry. The contour of the ventral opening of the RV pL was manually drawn by drawing successive lines up to the most caudal part of the image for both the dorsal and ventral walls of this cavity. We could therefore measure the local diameter versus the curvilinear abscissa along the line. Its length was calculated between the initial point up to the point of maximal diameter derivative. Note that this method probably overestimated the lumen diameter, as proper alignment between the dorsal and ventral walls was not perfect. Nevertheless, the ratio of about two between the diameter of this ventral opening of the pL of the RV and the central canal should be in the right order of magnitude.

The full-field optical coherence tomography (FF-OCT) setup was previously described (Sternberg et al., 2018; Thouvenin et al., 2017). The FF-OCT path is based on a Linnik interference microscope configuration illuminated by a temporally and spatially incoherent light source. A high power 660 nm LED (Thorlabs M660L3, spectral bandwidth 20 nm) provided illumination in a pseudo Köhler configuration. A 90:10 beam splitter separates the light into sample and reference arms. Each arm contains a 40x water immersion objective (Nikon CFI APO 40x water NIR objective, 0.8 NA, 3.5 mm working distance). In the reference arm, the light is focused onto a flat silicon wafer with a reflection coefficient of about 23.5% at the interface with water. FF-OCT detects and amplifies any structure that reflects or backscatters light within the sample arm, showing a contrast based on differences in refractive index, including small extracellular vesicles with rich lipid content. The beams from sample and reference interfere only if the optical path length difference between both arms remains within the coherence length of the system, ensuring efficient optical sectioning. A 25 cm focal length achromatic doublet focuses the light to a high speed and high full well capacity CMOS camera (Adimec, Q-2HFWC-CXP). The overall magnification of the FF-OCT path is 50x. The measured transverse and axial resolutions were 0.24 μm and 4 μm, respectively. Camera exposure was 9.8 ms, and images were acquired at 100 Hz. We acquired sequences of consecutive direct images and computed the standard deviation on groups of images to cancel the incoherent light that does not produce interference. Here, we also take advantage of phase fluctuations caused by the vesicle displacements occurring in the 10–100 Hz range to specifically amplify the signal while backscattered signal from the surrounding tissue is more static.

To induce expression of CD63-pHluorin, specifically labeling exosomes, embryos were injected at the 1000 cell stage. The following day, embryos were anesthetized with tricaine and embedded in 1.5% low melting point agarose. Live imaging recordings were performed at 28°C using a Nikon TSi spinning-wide (Yokagawa CSU-W1) microscope (Verweij et al., 2019). Other embryos were used for ultrathin cryosectioning and immunogold labeling, in which case they were fixed in 2% PFA, 0.2% glutaraldehyde in 0.1M phosphate buffer pH 7.4 at 3 dpf. Zebrafish were processed for ultracryomicrotomy and immunogold labeled against GFP using PAG 10. All samples were examined with a FEI Tecnai Spirit electron microscope (FEI Company), and digital acquisitions were made with a numeric camera (Quemesa; Soft Imaging System) (Verweij et al., 2019).

We used MATLAB to obtain the numerical resolution of the system of three coupled equations detailed in the Results subsection ‘A bidirectional flow generated by the asymmetric distribution of beating cilia’. Our simulations were performed by the finite element method solver COMSOL, using the Creeping Flow Module. The flow simulations were carried out by solving the Stokes equation in a rectangular domain composed of 120088 domain elements and 15046 boundary elements (corresponding to 630516 degrees of freedom). The simulations accounting for the transport of particles in the presence of a flow were performed by coupling the diffusion equation to the Stokes equation in the same rectangular domain, using the Transport of Diluted Species Module. Doubling the number of mesh elements did not affect the results obtained throughout this study, and all the simulations converged.

Key datasets and custom MATLAB scripts used in this study can be found on Dryad server (14-04-2019-RA-eLife-47699) and on Github (Thouvenin, 2019, Github, https://github.com/wyartlab/eLife_2019_OriginAndRole).

We here provide a more detailed derivation of the flow velocity profile in a model 2D central canal, with a diameter d and a ciliary region of height h. In the ciliary region (also called ventral region), the Stokes equation can be written:

while in the non-ciliary region (also called dorsal region), the Stokes equation does not contain the bulk force f_v, but only the common pressure gradient:

Integration of Equation A1 between y_A and y leads to:

where y = y_A corresponds to the y-position for which the shear dv/dy is zero. y_A is a constant that will be derived in the following. Similar integration in the dorsal region [B] yields:

where y_B is an analogous constant representing the y-position where a null shear is obtained.

In the ventral region [A], an additional integration of Equation A3, with the no-slip boundary condition (v = 0) at y = 0 (ventral wall) provides:

Similarly in the region [B], with the no-slip b. c. (v = 0) at y = d, we get from Equation A4:

In both Equations (A5) and (A6), three unknowns remain to be determined in order to get the velocity profile: dP/dx, y_A and y_B. At the border between the regions [A] and [B] (y = h), the continuity of velocity (v_A(y = h)=v_B(y = h)) and shear (dv_A/dy(y = h)=dv_B/dy(y = h)) can be written respectively:

and

The last required equation to solve the system is provided by the condition of null global flux across a section of the channel, which is expressed in this bidimensional geometry:

that is:

which finally provides:

The system of three Equations A7, A8 and A11 contains three unknowns dP/dx, y_A and y_B, and may be solved numerically. In the specific case of a ciliary layer occupying half of the central canal diameter (l = d/2), the solution is simple and symmetrical: dP/dx = f_v/2, y_A = d/4 and y_B = 3d/4 which leads to this symmetrical velocity profile:

We here provide further details concerning the derivation of the effective diffusivity D_eff of particles in the central canal in presence of a bidirectional flow.

In the CC, the motile cilia occupy the ventral half (Figure 2). From the flow model derived in Appendix 1 (Equations (A3) and (A4) of the main manuscript), one can rewrite the velocity profile changing the y-coordinate origin (y = 0 at the center of the channel) and using the maximal velocity V_max = 6 µm/s:

in the dorsal region (positive values of y), and:

in the ventral region (negative values of y). The full transport equation of particles of local concentration c writes:

where v(y) corresponds to the two velocity profiles written previously (Equations (A14) and (A15)), as the local velocity field is uniquely oriented in the rostro-caudal x-direction (see also Figure 3 of the main manuscript). D is the molecular diffusion coefficient, here taken as D = 10⁻¹¹ m²/s, corresponding to particles of diameter 40 nm using the Stokes-Einstein law (Equation A7 in the main manuscript). Several terms of Equation (A16) can be neglected considering their order of magnitude. With L the length of the channel, and d its diameter, the Equation (A16) may be written in dimensional analysis for times of order T ' L²/D = τ_diff (characteristic spanwise diffusion time):

or

Both V_maxL/D and L²/d² terms are much larger than unity, due to the high aspect ratio of the channel L/d. This allows to simplify Equation (A16) in the two dorsal (y > 0) and ventral (y < 0) regions, respectively, as follows:

and

Adapting the derivation of Bruus (2008), that was done for a Poiseuille flow in a cylinder, we first assume that the axial derivative of the concentration ∂_xc is independent on y. This assumption will be verified a posteriori and confronted to numerical simulations. Consequently, Equations (A20) and (A21) are ordinary differential equations for c(y), that can be solved. A first integration between y and y = d/2 for Equation (A20) for y > 0 leads to

Similarly, integrating Equation (A21) between y = −d/2 and y for y < 0 leads to

By antisymmetry, ∂_yc |_y=d/2= −∂_yc |_y=−d/2. The equality of the concentration y-derivative in y = 0 then leads to ∂_yc |_y=d/2= −∂_yc |_y=−d/2= 0, which simplifies Equations (A22) and (A23). A second integration between y = 0 and y of Equation (A22) for y > 0 leads to

Similarly, a second integration between y = 0 and y of the Equation (A22) for y < 0 leads to

We now aim at expressing the concentration profile with respect to the average concentration $\stackrel{-}{c}$ over a cross section

which, using Equations (A24) and (A25), yields

and finally

This theoretical prediction is confirmed by numerical simulations solving the complete transport equation without the assumptions made for the theoretical model (Appendix 2—figure 1), for different times, positions in the canal and flow velocities.

Appendix 2—figure 1

Download asset Open asset

Going back to our theoretical model, one can now derive the condition for having a y-independent axial gradient of concentration (condition required for recovering a diffusive-like concentration profile) by differentiating the Equations (A24) and (A25) with respect to x. The latter condition is consequently fulfilled only if ∂_xc = ∂_xc(0), which requires the third term to be negligible. In terms of order of magnitude, this implies that

where L is the length of the central canal. This condition can be rewritten as a condition on the Péclet number Pe:

This condition on the Péclet number is always fulfilled in the CC of zebrafishes for nanobeads or lipidic vesicles (D ∼ 10⁻¹¹ m²/s), where Pe does not exceed 10.

We now aim at deriving the effective diffusivity of particles in the CC. We calculate the average current density J(^¯x) through the cross-section, using the velocity profiles (Equations (A14) and (A15)) and the concentration profiles (Equations (A24) and (A25)):

The zero average flux over the cross section leads to the vanishing of the constant component $\stackrel{-}{c}$ of Equations (A24) and (A25). The previous expression may then be simplified:

which may be approximated to

The effective diffusivity D_eff may then be directly derived from the effective Fick’s law

A more detailed calculation in a cylindrical channel was done by Aris (1956). It generalizes this result at intermediary Péclet numbers Pe ~ 1. Adapted to our geometry and flow, it yields

For particles of radius 20 nm (corresponding to the size of injected nanobeads), we have D ~ 5·10⁻¹² m²/s, which leads for V_max = 6 µm/s to an effective diffusivity D_eff ~2.5D. Similarly, the effective diffusivity of the same particles for a faster flow of amplitude V_max = 20 µm/s is D_eff ' 18D. These two theoretical predictions are in good agreement with the numerical simulations (Figure 6C of the main manuscript).

One strong assumption of the model that enabled to readily integrate Equations (A7) and (A8) considers that the axial derivative of the concentration ∂_xc is only weakly varying over a cross section, such that one can assume it is constant. The Appendix 2—figure 1 plots the relative variation of the normalized concentration gradient ∂_xc/∂_xc (0) as a function of y. It shows that this quantity exhibits minor variations of about 10% over a cross section, demonstrating the relevance of the assumption in the derivation.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank Prof. Rebecca Burdine for the kurly (cfap^tm304 allele) mutant, Prof. Jarema Malicki for the elipsa mutant (ift54/traf3ip) and Prof. Brian Ciruna for the Tg(beta-actin:Arl13b-GFP) line. We would like to acknowledge the kind and expert assistance of Aymeric Millecamps, Tudor Manoliu, and Basile Gurchenkov from the ICM. Quant imaging facility for instrument use and scientific assistance on optical imaging, as well as Sophie Nunes-Figueiredo, Bogdan Buzurin, and Monica Dicu from PHENOZFish platform for fish care. We are grateful to Martin Catala for advising us on zebrafish histology. We greatly thank Jenna Sternberg for article proofreading. We thank the Yaksi laboratory members and the fish facility support team for scientific and technical assistance to NJK. This work was funded by Human Frontier Science Program (HFSP) Research Grant (grant n° RGP063-2018), and the New York Stem Cell Foundation (NYSCF) (grant n° NYSCF-R-NI39) and an ICM postdoctoral fellowship kindly attributed to OT. The research leading to these results has also received funding from the program ‘Investissements d’avenir’ ANR-10-IAIHU-06 (Big Brain Theory ICM Program), ANR-11-INBS-0011 (NeurATRIS: Translational Research Infrastructure for Biotherapies in Neurosciences), and a Helse Midt-Norge Samarbeisorganet grant to NJY and the Kavli Institute for Systems Neuroscience at NTNU.

Animal experimentation: All procedures were performed on zebrafish embryos before 2 days post fertilization in accordance with the European Communities Council Directive (2010/63/EU) and French law (87/848) and approved by the Brain and Spine Institute (Institut du Cerveau et de la Moelle épinière).

1. Received: April 15, 2019
2. Accepted: January 7, 2020
3. Accepted Manuscript published: January 9, 2020 (version 1)
4. Version of Record published: January 29, 2020 (version 2)

© 2020, Thouvenin et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.