Selection mechanism at the onset of active turbulence

Active turbulence describes a flow regime that is erratic, and yet endowed with a characteristic length scale1. It arises in animate soft-matter systems as diverse as bacterial baths2, cell tissues3 and reconstituted cytoskeletal preparations4. However, the way that these turbulent dynamics emerge in active systems has so far evaded experimental scrutiny. Here, we unveil a direct route to active nematic turbulence by demonstrating that, for radially aligned unconfined textures, the characteristic length scale emerges at the early stages of the instability. We resolve two-dimensional distortions of a microtubule-based extensile system5 in space and time, and show that they can be characterized in terms of a growth rate that exhibits quadratic dependence on a dominant wavenumber. This wavelength selection mechanism is justified on the basis of a continuum model for an active nematic including viscous coupling to the adjacent fluid phase. Our findings are in line with the classical pattern-formation studies in non-active systems6, bettering our understanding of the principles of active self-organization, and providing potential perspectives for the control of biological fluids. Experiments on microtubule-based nematics, together with active gel theory, suggest that the length scale associated with active turbulence is selected at its onset—balancing activity with the stabilizing effects of nematic elasticity and geometry.

work by Duclos et al., which reports spontaneous flows for spindle-shaped cells confined in stripes 18 . However, these flows do not exhibit an intrinsic length scale, but are system size dependent and thus not enforced by a genuine wavelength selection mechanism. Here, we demonstrate from experiments with microtubule-based active nematics and active-gels theory that the length scale at fully developed turbulence is selected right at the onset of the initial reorientational instability of a radially aligned active nematic free of confinement. In our experiments, the selection mechanism is understood in terms of a balance between activity, as a destabilizing force acting over the whole wavenumber range, and a pair of stabilizing effects acting distinctively at both ends of the spectrum: at the largest wavenumbers, stabilization is associated with the nematic elasticity of the material, and in the smallest range, both a geometric effect of the radial symmetry and hydrodynamic coupling to the adjacent viscous fluids suppress the long wavelength instability.
Our chosen material is based on the self-assembly of micrometre-sized stabilized microtubules condensed at an oil-water interface from a bulk active-gel preparation 4,5,19 (see Methods and Supplementary Fig. 1). The microtubules are bundled under the depleting action of polyethylene-glycol (PEG). This facilitates crosslinking and internal shearing by kinesin motors clustered with streptavidin (see Methods). Fuelled with adenosine triphosphate (ATP), the bundles are permanently subjected to extension and buckling. Left to evolve free of constraints, the two-dimensional (2D) active layer adopts a characteristic long-range orientational (nematic) order, continuously permeated by large-scale flows and locally interrupted by regions void of microtubules that configure semi-integer defects. The above scenario, which displays continuously renovating dynamics of defect creation and annihilation [20][21][22] , corresponds to a state of active nematic turbulence, the origin of which we unveil here.
We begin by aligning the active nematic in contact with a glass tube vertically placed at the oil-water interface ( Supplementary  Fig. 1). In-flow currents flowing into the tube align the nematic active sample with bundles lying radially around the area where the tube is positioned (Fig. 1a). After removal of the capillary, this preparation is unstable, being immediately disrupted by spontaneous buckling of the aligned extensile material, leading to the development of a bend-type distortion (Fig. 1b). The instability organizes into a pattern of periodically spaced concentric crimps that evolve into circular walls. Elastic stress accumulated in the walls is released through the nucleation of pairs of ± 1/2 defects 13 that align along dark circular lanes (Fig. 1c). Subsequent defect motion dismantles the circular pattern (Fig. 1d), and reverts the system into a turbulent regime (Supplementary Video 1).

Selection mechanism at the onset of active turbulence
Berta Martínez-Prat 1,2 , Jordi Ignés-Mullol 1,2 *, Jaume Casademunt 3,4 and Francesc Sagués 1,2 The development of the fully turbulent regime is better visualized in situations of low motor concentration (Supplementary Video 2). In these cases, the instability repeats periodically following orthogonal directions until the material arrives at the characteristic turbulent state. First, self-propelled defects remain confined within the circular lanes that originated in the original bend instability. Defect motion exerts a shear flow that forces the progressive tangential alignment of the bundles between neighbouring lanes. In this aligned configuration, the system is again prone to bend instability, which this time organizes radially directed walls, along which new motile defects unbundle and propel ( Fig. 2a-d). As the wavelength of the instability remains practically unaltered, the signature of the whole process is illustrated by lines of moving defects that form a square grid after time averaging the recorded fluorescence micrographs (Fig. 2e). Alternation between radial and azimuthal modes is evidenced by the trace of the corresponding fast Fourier transform (FFT) amplitudes (Fig. 2f).
This pattern-forming instability can be quantitatively described in terms of the characteristic wavenumber k* and growth rate Ω* corresponding to the leading (most unstable) mode of distortion (see Methods and Supplementary Figs. 2 and 3). Four different control parameters were considered: ATP, motor units, microtubules and PEG. ATP concentration is considered a direct measure of the chemical energy put into the system, and the other three parameters are presumably related to the power effectively transduced on the microtubule network and also its mechanical properties.
The first set of results are taken on varying the ATP concentration. Plots of both k* and Ω* indicate that the instability is 'stronger' (that is, a smaller imprinted wavelength or, equivalently, larger wavenumber, and a faster response) with increasing ATP concentration (Fig. 3a). Experiments with microtubule-based active nematics have shown that activity α is related to a logarithmic measure of the concentration of ATP 19 . We will thus use α ≃ ln[ATP] for the experimental activity in our analysis, assuming all other control parameters are kept constant. For very low α, the extracted k* values collapse into a plateau, as the corresponding wavelengths become of the order of the sample size (a few hundreds of micrometres). Past a threshold activity, we find that our data are compatible with both (k*) 2 and Ω* increasing linearly with α. The role of motor units (Fig. 3b) is evidenced by the direct linear dependence of both k* and (Ω *) 1/2 on increasing concentration of clustered kinesins. The dependence on microtubule content shows that k* and Ω* decrease past a range of minimal microtubule concentration (Fig. 3c). Finally, non-monotonous behaviour is observed when varying the concentration of the PEG depleting agent (Fig. 3d). At small concentrations, an increase in [PEG] renders the system unstable towards smaller wavelengths. Conversely, high concentrations of PEG render the material more rigid, and it is only able to accommodate distortions with wavelength comparable to the system size.
Interestingly, the separated data sets collapse into a single master curve, Ω* ≈ (k*) 2 (Fig. 3e), that fully characterizes the involved length and timescales at the onset of the bending instability of the active nematic. Wavelength selection is rationalized within the theory of active gels 23,24 . We build on the analysis in ref. 25 for parallel alignment. In an unconfined system and for extensile stresses, this gives a linear growth rate of the form Ω(k) = a′ − b′ k 2 for longitudinal perturbations (see Supplementary Information), which yields no intrinsic wavelength selection. In the case of radial alignment, the growth rate spectrum is expected to be discrete and hence it is necessary to select a finite scale with maximal growth rate. However, the effect of the unavoidable viscous coupling to liquids in contact with the active nematic fluid 26 has been shown to qualitatively modify the growth rate for parallel alignment into the form Ω(k) = ak − bk 2 (refs. 27,28 ), which indeed provides a genuine selection mechanism at the maximal growth rate. Therefore, for the benefit of discussion and physical insight, we will assume the latter form of growth rate as an appropriate qualitative description valid also for the radial case.
Parameter a of the linear part in the dispersion relation should depend on activity, but in a way that is not known precisely, while b, characterizing the equilibrium elastic relaxation, is explicitly known in terms of the material parameters (see Supplementary  Information). A critical wavenumber k c = a/b, above which random perturbations decay, sets the upper bound on the instability range, while the maximum growth rate is set at k* = k c /2, independent of the material parameters. Accordingly, we may assume that, for the real system at hand, we generically have k* = θk c , where θ < 1 is a numeric factor that is to be determined, but is essentially independent of the parameters. From the experimental measurements of k * we may estimate a from the relationship a = bk c ; however, because the value of k c should approach that in the absence of viscous coupling in the appropriate limit we may expect a ∝ (α) 1/2 (see Supplementary Information).
Interestingly, the growth rate of the selected mode Ω* = Ω(k*) can be written in the form Supplementary Information). Notice that this is the well-known scaling reported in the literature for the active length scale. Using the explicit form of b (see Supplementary Information), this results in Here, K, γ, η and ν are, respectively, the elastic constant, the rotational and shear viscosities, and the flow-alignment coefficient. From the above considerations, we can readily interpret the experimental results when varying the ATP content (that is k* ∝ α 1/2 ). For kinesin concentration, the experimental scaling can be understood by referring to earlier studies on the velocity of the active microtubule-based material. First, experiments 29  has been predicted theoretically and confirmed experimentally 1,19 . Consequently, the linear scaling of k* with motor concentration (Fig. 3b) is fully consistent. Concerning microtubule content, the behaviour (Fig. 3c) can be interpreted by assuming, plausibly, that when densifying the active nematic by increasing the microtubule concentration, the system becomes stiffer (larger K), making k* and Ω* smaller. Finally, we can now interpret the two regimes observed when changing the concentration of PEG. At low concentrations, PEG makes the system more efficient by optimizing motor translocation, thus resulting in an enhanced effective α. At high PEG concentrations, our results highlight that, beyond a threshold, PEG has a similar effect to the one observed for the case of increasing microtubule concentration.
One of the main results of this analysis is equation (1), which justifies the observed scaling with a constant of proportionality that depends solely on material parameters (K, γ, η, ν) but not on the activity. In this way, we conclude that there is also a characteristic timescale of the emerging pattern, which is given as a simple ratio  Fig. 3e).
of the activity and shear viscosity parameters (see Supplementary  Information), Ω ∝ α η * . We identify this timescale as characteristic of the turbulent regime, as previously suggested 30 .
The presented results demonstrate a route to active turbulence in nematics endowed with characteristic length and timescales that are imprinted right at the onset of an extensional instability of the material. These results should motivate a renewed interest in better understanding the fundamentals and realizations of active flows, both in the laboratory and in natural systems.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, statements of data availability and associated accession codes are available at https://doi.org/10.1038/ s41567-018-0411-6.
Kinesin expression. The kinesin used in the experiments was the heavy-chain kinesin-1 K401-BCCP-6His from Drosophila melanogaster (truncated at residue 401, fused to biotin carboxyl carrier protein (BCCP), and labelled with six histidine tags). This was expressed by Escherichia coli using the plasmid WC2 from the Gelles Laboratory (Brandeis University) and then purified with a nickel column. Finally, the suspension was dialysed against 500 mM imidazole buffer and stored in a 40% wt/vol sucrose solution at − 80 °C. The final concentration of kinesin was estimated by means of absorption spectroscopy.
Kinesin/streptavidin clusters. Kinesin/streptavidin clusters were assembled by mixing freshly thawed K401 kinesin with streptavidin (Invitrogen 43-4301) at a ratio of 2:1 with a streptavidin concentration of 30.6 μ g ml −1 . DTT was added to a concentration of 2 μ M. The mixture was incubated on ice for 30 min.
Assembly of the active nematic. The standard preparation consisted of a M2B aqueous mixture with ATP (1.5 mM), kinesin/streptavidin clusters (8.2 μ g ml −1 of streptavidin), microtubules (1.3 mg ml −1 ) and the depleting agent PEG (1.7% wt/wt). To maintain constant the concentration of ATP during experiments, we also added an ATP-regenerating system consisting of pyruvate kinase/lactate dehydrogenase (Sigma, P-0294) (6% vol/vol) and phosphoenol pyruvate (27.8 mM). The mixture also contained antioxidants to avoid photobleaching during the fluorescence imaging: DTT (5.8 mM), glucose oxidase (0.2 mg ml −1 ), catalase (0.04 mg ml −1 ), trolox (6-hydroxy-2,5,7,8-tetramethylchroman-2-carboxylic acid) (2.1 mM) and glucose (3.5 mg ml −1 ). To maintain a biocompatible oil-water interface in subsequent steps, we also included in the mixture the surfactant pluronic (Sigma P2443) (0.4% wt/vol). In some cases we varied the concentrations of ATP, motor clusters, PEG or microtubules to study the influence of their concentrations. To prepare the samples, 1 μ l of the final microtubule suspension was introduced into a polydimethylsiloxane (PDMS) pool (5 mm diameter and 2 mm depth) already filled with an isotropic oil (PDMS) with a kinematic viscosity of 1.0 cm 2 s −1 (Rhodorsil Oils 4747V 100 Bluestar Silicones) (Supplementary Fig. 1). The pool was not covered, leaving the sample open. Spontaneous adsorption of the microtubules onto the oil-water interface led to the formation of the active nematic. The pool was made by polymerizing PDMS inside a custom-made mould and finally gluing it on a polyacrylamide-functionalized slide with ultraviolet glue (Norland NOA-81).
Observation of the active nematic. Observation of the active nematic layer was carried out by fluorescence microscopy. We used a custom-built inverted optical microscope equipped with a white LED source (Thorlabs MWWHLP1) and a Cy5 filter cube (Edmund Optics). Images were captured using a Zyla 4.2 PLUS sCMOS camera (Andor) operated with the open-source software ImageJ μ-Manager.
Alignment of the active nematic. A capillary tube (Sigma, Z611239), with outer and inner diameters of 1.2 mm and 0.3 mm, respectively, was held with a linear motion stage (KDC101, Thorlabs) above the pool containing the sample, and it was slowly moved until it touched the oil-water interface, creating a radial flow that led to alignment of the active nematic layer (Supplementary Fig. 1). Afterwards, the tube was carefully removed, allowing observation of the interface and enabling the active nematic to evolve freely.
Characterization of the bend instability wavenumber k*. A time-averaged micrograph of different frames displaying the pattern was used to measure the fluorescence intensity profile along a direction orthogonal to the crimps, using ImageJ ( Supplementary Fig. 2). Finally, the wavenumber was extracted by performing a FFT of the intensity profile using MatLab. Error bars were calculated as the s.d. of ten measurements along different radial directions in a given experimental realization.
Characterization of the bend instability timescale Ω*. The instability timescale was determined by monitoring the time evolution of the FFT power spectra with ImageJ. Ω* was finally extracted by fitting the FFT intensity of the known k* at short times to an exponential trend ( Supplementary Fig. 3). Error bars correspond to uncertainty in the fitted parameter.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data that support the plots within this paper and other findings of this study are available from the corresponding author on request.

Corresponding author(s): Jordi Ignés-Mullol
Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see Authors & Referees and the Editorial Policy Checklist.

Statistical parameters
When statistical analyses are reported, confirm that the following items are present in the relevant location (e.g. figure legend, table legend, main text, or Methods section).

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement An indication of whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistics including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers upon request. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability The data that support the plots within this paper and other findings of this study are available from the corresponding author on request