Wrinkling Instability in 3D Active Nematics

In nature, interactions between biopolymers and motor proteins give rise to biologically essential emergent behaviors. Besides cytoskeleton mechanics, active nematics arise from such interactions. Here we present a study on 3D active nematics made of microtubules, kinesin motors, and depleting agent. It shows a rich behavior evolving from a nematically ordered space-filling distribution of microtubule bundles toward a flattened and contracted 2D ribbon that undergoes a wrinkling instability and subsequently transitions into a 3D active turbulent state. The wrinkle wavelength is independent of the ATP concentration and our theoretical model describes its relation with the appearance time. We compare the experimental results with a numerical simulation that confirms the key role of kinesin motors in cross-linking and sliding the microtubules. Our results on the active contraction of the network and the independence of wrinkle wavelength on ATP concentration are important steps forward for the understanding of these 3D systems.


■ INTRODUCTION
Active matter is a state of matter where large-scale dynamical structures emerge from the interaction of individual active components each driven by their own internal energy source. A well investigated class of such systems comprises active nematics 1,2 where individual elongated components selforganize into a dynamic state with spatiotemporal chaos with topological defectsa behavior known as active turbulence. 3−10 Several experimental and theoretical model systems have been developed to study networks of cytoskeletal filaments and motor proteins fueled by ATP. 11−13 A remarkable example is the emergent behavior in 2D observed in mixtures of microtubules and kinesin motors arranged in bundles by a depletion agent pioneered by Sanchez et al. 12 Internally generated spatiotemporally chaotic flows at the millimeter scale were observed that persisted as long as ATP was available. The 2D networks exhibit a steady state with permanent flow at large scales, while the dynamics is driven by the buckling and elongation of the microtubule bundles due to the activity of the motors on the nanometer scale. 14 At high concentrations these networks show a transition to a locally ordered nematic liquid-crystalline state with topological defects and spatiotemporal dynamics 15−21 that can be varied by confinement such as structured liquid interfaces, 22 hard wall circular boundaries, 23 curved surfaces of deformable spherical vesicles, 24 and toroidal droplets 25 as well as by radial alignment. 26 These systems share the experimental approach to be condensed at an oil−water interface.
Whereas systems investigated in two dimensions were the focus of copious studies, the prospect of developing threedimensional active matter systems was shown only by a few examples. These 3D experiments included microtubule-based active gels confined in toroids, 27 and sedimented cell-sized droplets. 28 In the present study we show experimentally that a threedimensional suspension of microtubules and kinesin-1 motors self-organizes into an active nematic ribbon that for a range of parameters wrinkles periodically in the third-dimension due to its longitudinal extension before transitioning to a state of active turbulence. Similar results of a parallel study have been recently published. 29 However, the contraction in the mentioned system is passive and due to depletion agent, while our network actively contracts due to the motor action. A related folding phenomenon has been observed also in thin acto-myosin networks 30 whose behavior was generated by contracting forces and not by extensile ones as in the present work. Additionally, we understand the behavior of our system with stochastic simulations and a quantitative theory that describes the wrinkling instability.

■ RESULTS AND DISCUSSION
We experimentally investigate the dynamics of an active nematic, confined in a long rectangular channel of size 30 mm × 1.5 mm × 100 μm. The active nematic consists of an ionic aqueous solution of microtubules, multiheaded kinesin-1 motors, and poly(ethylene glycol) (PEG) as depletion agent. 12 We observe that when the system is initially prepared with nematic order 31 along the long axis of the channel, it exhibits a sequence of transitions before developing into the active turbulent phase with spatiotemporal chaos. This is summarized in Figure 1. In particular, the system makes a transition from a space-filling 3D active nematic to a narrow ribbon extending over the length of the channel and located in the midplane (away from the boundaries). This initial contraction along the z-direction, which occurs over a time scale of τ c , is followed by a wrinkling instability at time scale τ w , which takes the system back to being 3D. The long-lived state of spatially periodic wrinkling that is oriented along the long axis undergoes another transition to an active turbulent state where the active nematic fills the entire channel after a long time τ at . This remarkable sequence demonstrates 3D active nematic dynamics and instabilities similar to those observed in a recent study. 29 Conversely, no comparison can be made with 2D systems investigated so far 26 due to different spatial distribution of that system and its interaction with substrates. Inside the ribbon, nematically ordered filament bundles are clearly visible. To confirm the alignment of bundles, we quantify the nematic order in the sample by using correlation analysis and reconstruct the nematic director field as can be seen in Figure 2. The first mode of evolution of the active nematic is contraction perpendicular to the nematic director. After contraction, extensile stress builds up along the nematic  Nano Letters pubs.acs.org/NanoLett Letter director 1 ( Figure 1A). Here, contractile and extensile stresses are internal stresses that would lead to contraction and extension, respectively, in the absence of constraints. In the channel, long-range expansion along the whole channel length is not possible and the extensile stress can lead to a spatially periodic wrinkling of the ribbon (Figure 1B−E and Movie S1). No in-plane bending (horizontal deformation of the director field, akin to deformations in 2D active nematics) is observed at this stage. We expect that out-of-plane bending is preferred to in-plane deformations because it is opposed mainly by the stiffness of the microtubules whereas the in-plane bending additionally requires shearing of the cross-linked sheet. By acquiring consecutive areas along the x-axis, we showed that the 3D instabilities are not local but rather emerging in the whole network over centimeters ( Figure S3). As the wrinkles grow, they reach the channel walls and start folding over. At this point, the sheetlike structure disappears and the individual bundles continue to extend in different directions in space (Movie S3). After the time τ at , the system reaches a dynamic state of active turbulence ( Figure 1E). The characteristics of the material needed to form a 3D active nematic differ significantly from those in the previously studied 2D situations. 12,26 Those studies typically use short microtubules with lengths around 1 μm, higher depletant concentration, an oil−water interface and a radial arrangement. Nematic ordering and active stress buildup in our 3D systems have a relatively low volume fraction of microtubules (0.001) and require longer filaments. The average microtubule length in our experiment was 19 ± 10 μm ( Figure S2). We verified the importance of length by shearing the filaments to much shorter lengths, in which case no pattern was observed. Contraction Characterization. The evolution of the system from 3D to 2D includes contraction along both yand z-directions. The average lateral (y) contraction was 169 ± 8 μm or 11.2%. Vertically (z-direction), the ribbon thickness reaches 30 ± 2 μm, corresponding to a shrinkage by 70%. The nonuniform contraction likely reflects the dimensions of the channel. Lateral contraction involves larger absolute displacements and is therefore opposed by stronger drag forces. Similar nonuniform contraction is also seen in the simulation ( Figure  5H,I). Whereas Senoussi et al. 29 attribute the contraction to the effect of the depletion agent alone (2% pluronic in that study), we show that in our system it results from the active component due to cross-linking activity of kinesin motors combined with the passive bundling effect of PEG. By carrying out the experiment without including motors to the filament mixture, we could observe the bundling of microtubules due to PEG alone, their nematic orientation along the longitudinal axis of the channel and the y-contraction ( Figure 3A,B). The network contracts by 7.6% and reaches its final, stable state after ∼330 min. No instabilities are present. The final thickness of ca. 40−45 μm corresponds to a contraction of 57%. These contractions can be explained by the entropic mechanism of the depletion interaction. The nonadsorbing polymer (PEG) bundles the microtubules in order to minimize the excluded volume. 32 The dimension of the resulting bundles depends on the molecular weight and radius of gyration of PEG (in our case20 kDa and 7 nm, respectively).
In order to show the active contribution to the contraction, we repeated the experiment with motor proteins alone and without PEG. Interestingly, the system actively contracts as in the case with depletion agent alone ( Figure 3C,D) and reaches a maximum shrinkage of 36%, much stronger then the passive contraction through depletion forces. However, the time scale for shrinkage was 3 times slower than in the passive case. Although in this case the microtubules are ordered in a parallel fashion, neither neither thick bundles nor instabilities are visible ( Figure S4). We explain the absence of instability by the missing bundling effect of PEG, which normally reduces the relative distance between microtubules to nanometric gap comparable to the dimension of motor clusters. This represents an optimal value for the continuous binding of kinesin motors and therefore for their efficacy. The development of the extensile stress that is needed for the wrinkling lso requires the combined action of the multiheaded kinesins and the depletion forces. Again the PEG-induced force that holds the filaments together reduces the distance between the microtubules within the bundles, enhances the binding of the motors to the polymers and gives the system the possibility to contract. As a positive feedback, reduced distance enhances the cross-linking of the motors and leads to further contraction. Under this configuration, the system reaches the ribbon arrangement and the motors extend the polymers randomly oriented within the ribbon. The extensile stress can be explained by an asymmetric distribution of kinesin motors along microtubules. Because kinesins move toward the microtubule plus ends and then exert a minus-directed force on them, there will be an excess of compressive loads on the microtubules. These, in turn, lead to an extensile stress in the sheet (see the section Simulation).
Wavelength Analysis. In our experiments at the onset of the wrinkling instability, the observed wavelength λ ranged from 200 to 500 μm for constant ATP concentration (2 mM) and the time until the onset of instability (τ w ) varied between 10 min and 2 h ( Figure 4D). We find that the wavelength of the pattern and the time scale of its formation are correlated; Nano Letters pubs.acs.org/NanoLett Letter wrinkling instabilities that appear at shorter times tend to have shorter wavelengths than those forming later. To further investigate this correlation, we varied the ATP concentration from 2 mM, at which the motors will be reaction-limited, to 10 μM that is in the diffusion-limited regime. At 1 mM ATP concentration the wavelength of the pattern was 272 ± 36 μm and the wrinkling set in after few minutes; at 10 μM ATP the wavelength was 249 ± 38 μm and the onset of wrinkling occurred at around 13 min. Contrary to what one might expect, no discernible effect on the wavelength of the wrinkling instability could be observed. In other studies a 2D system shows a faster response and a shorter selected wavelength with increasing ATP concentration. 26 However, due to the differences between the two setups (2D system, high depletant concentration, confinement of the system to an oil−water interface), no direct comparison can be made. Remarkably, in our case, the instability does not change with ATP concentration as the time for wrinkling and its wavelength at 10 μM fall into the range observed at 2 mM ATP. This is an interesting result as the variation of the characteristic wrinkling formation times cannot be understood through the biochemical process linking velocity and force to ATP hydrolysis, as is the case in classical 2D active nematics. 33 It was not observed in the other 3D studies as those systems were investigated at constant ATP concentration at reaction-limited regime. 29,30 These observations suggest that the active stress is determined by the stall force of kinesin motors, which is independent of the ATP concentration, and by the crosslinking nature of the motors. To test the importance of the motor concentration, we decreased the kinesin concentration by 10-fold (to 1.7 nM) and indeed did not observe the instability ( Figure 4E). However, the network contracted continuously for around 48 h, resembling the contraction of the stabilized microtubule network in Xenopus oocyte extracts. 34 This result proves that the active force exerted by a lower amount of motors is able to contract the network through the capacity of kinesin to cross-link the filaments but is not strong enough to trigger the instability. Theory. We can develop a theoretical understanding of the observed phenomena based on the combination of the contraction and the mechanical activity inside the resulting ribbon. The director field is mostly aligned in the direction of the channel, n̂≃ ex. For microtubules and kinesin motors, we expect an extensile stress along the director field. This stress is accompanied by contractile stresses across the channel (directions eŷ and eẑ). After the contractile stress in the zdirection leads to the formation of the ribbon, we can describe the system as a two-dimensional active nematic sheet that can explore out-of-plane bending deformations. For this effective description, the nematic order is described by the symmetric with an activity coefficient ζ. In principle, an isotropic contribution to the stress can also exist, but it can be shown that this contribution plays a subdominant role and can be ignored for simplicity (see the Supporting Information). The ribbon can be represented by a quasi 2D sheet, whose shape can be described (in a Monge representation) by a function h(x, y) for the vertical deflection relative to the center of the channel, and a sheet stiffness κ eff that can be derived from the bending modulus of the individual filaments and the nematic order parameter (see the Supporting Information).
Focusing for simplicity on the variations along the x-axis, we can construct a generalized free energy functional as where the coupling between the active stress and the deformation arises from the nonlinear contribution to the strain tensor (see the Supporting Information for a full derivation that includes deformation in both directions). The emerging instability of the sheet is subject to viscous dissipation (either in the fluid or in the remaining microtubule network outside of the sheet). We describe it with a local drag term that resists vertical movement of the sheet, which is written as Γ∂ t h. In the Supporting Information we show a posteriori that it largely exceeds the expected hydrodynamic drag in the fluid layer surrounding the sheet and we therefore attribute it to the residual filament network filling the volume. Within a linear stability analysis, we use the ansatz h(x, t) ∝ e Λt e iqx for a perturbation with the wave vector q. We find a growth rate of For negative values of the active stress corresponding to extensile stress, the growth rate Λ(q) will be positive for an interval of wave numbers, with the fastest growing mode having the wavenumber q S 4 eff * = ζ κ and a growth rate Assuming that the samples do not differ in other parameters, we obtain a relationship between the characteristic time of pattern formation and the wavelength λ = 2π/q* as follows w 4 τ λ ∝ (5) The experimental data shown in Figure 4D confirm the trend predicted by eq 5, both when comparing the outcome in equal samples and in those with different ATP concentrations. With the above equations we can quantitatively estimate the active stress per filament (see the Supporting Information) as 0.004 pN, several orders of magnitude below the force of a single kinesin motor. This indicates that the motor forces almost cancel out in the network and only a small imbalance contributes to the macroscopic stress.
Our results also suggest that the nematic order is the requirement for wrinkling and that the instability is independent of the size of the channel. This explains the coexistence of buckling along perpendicular axes in adjacent areas of the same channel shown in Figure 4F. In this case, their nematic domains are oriented along perpendicular directions.
Simulation. We used a complementary computer simulation to directly verify that sheet formation and wrinkling arise from the interplay at nanometer scale between the action of kinesin motors and attractive depletion forces. We simulated the dynamics of 20 000 filaments in a sufficiently large box to reproduce the salient phenomena (see the Supporting Information). The simulated system consists of elastic filaments (microtubules), tetrameric kinesin motors (pairs of kinesin dimers that can bind to two adjacent filaments), and passive cross-linkers. Filaments are subject to hard-core repulsion and attractive depletion forces. The simulations use a Langevin-dynamics algorithm from the Cytosim package. The results are shown in Figure 5, Figure S6, and Movie S4.
Starting from a nematically ordered state ( Figure 5A), the simulation first shows the condensation in the middle of the channel and shrinking in the lateral direction ( Figure 5B). As in the experiment, formation of microtubule bundles also becomes visible at this stage. Afterward, the extensile stress ( Figure 5G) leads to a wrinkling instability in the sheet ( Figure  5C). The growth of the selected out-of-plane deformation mode continues until the amplitude becomes limited by the simulation box size ( Figure 5D). Afterward, the ribbon structure disintegrates and the system ends in a state of 3D active turbulence ( Figure 5E,F). The simulation shows a very good qualitative agreement with the experimental observations. The lateral and vertical distribution of filaments as a function of time are shown as kymographs in Figure 5H,I, respectively. The maximal contraction is 73% vertically and 19% horizontally. The contraction is not uniform but proceeds from the borders toward the center, explaining why the horizontal contraction is relatively weaker, although the absolute displacement is larger.
The simulations are in agreement with our earlier conclusion that the kinesin motors responsible for stress generation are close to the stall conditions. At the onset of instability, the kinesins connecting antiparallel filaments act with an average force of around 4pN or 80% of the stall force ( Figure S6B). Furthermore, the simulation provides an insight into the mechanism behind the buildup of extensile stress. Because motors move toward the microtubule plus ends, a density gradient is established along the microtubules ( Figure S6C). This, in turn, causes, on average, a compressive force in the microtubule or an extensile stress in the network.
We also repeated the simulations in the absence of motor proteins ( Figure S6E) and without the attractive interaction caused by PEG ( Figure S6F). Without motors, bundles are still formed and the network contracts, albeit slower than with motors. Without attractive forces, the contraction is further slowed down and no bundling is visible.

■ CONCLUSION
We present a 3D active system without substrate interactions that collapses into a sheet due to activity and further evolves to exhibit a 3D wrinkling instability before transitioning to fully developed active turbulence. The active stress is determined by the stall force of the molecular motors and simultaneous crosslinking function as the simulation confirmed. We expect our observation and explanation of these novel emergent properties to broaden the perspective of active nematic systems in 3D.

■ MATERIALS AND METHODS
Motile Bundle Solution. The motile bundle solution was prepared as described before. 36 Briefly, the kinesin-streptavidin complexes were mixed with 2.7 mg/mL tubulin and 0.6% PEG. (See the Supporting Information for a more detailed description.) By adding PEG, attractive interactions between microtubules are induced through the depletion force and lead to bundle formation. 37,38 After the polymerization of microtubules, the active solution was injected into a PDMS microfluidic channel with a rectangular cross-section. To reduce unspecific protein adsorption, the channel was functionalized with PLL-g-PEG (see the Supporting Information). The peak-to-peak amplitude of the wrinkling can be determined by vertical scans (Movie S2) as the difference between the two z-positions at which the structure is in focus ( Figure 1F,G).
Imaging and tracking. Image acquisition was performed using an inverted fluorescence microscope Olympus IX-71 with a 4× objective (Olympus, Japan) and a DeltaVision imaging system (GE Healthcare). The images were acquired with a variable frame rate according to the experiment with an exposure time of 500 ms for a variable time according to the experiment (see the Supporting Information for more details).
Simulation. We simulated the evolution of a 3D active nematic system using the open source Cytosim package 39 (www.cytosim.org) (see the Supporting Information for more details).

* sı Supporting Information
The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.nanolett.0c01546. S1: Experimental results of the 3D wrinkling instabilities in microfluidic channel. The transition from a 3D system Nano Letters pubs.acs.org/NanoLett Letter to a 2D ribbon wrinkling in the direction perpendicular to the director field is clearly visible (MP4) S2: Experimental results of the 3D wrinkles scanned along their peak-to-peak amplitude (z-direction). Their height is about 70 μm whereas the channel is 100 μm high (MP4) S3: Experimental results of the transition from a ordered nematic to a chaotic state (MP4) S4: Computer simulation of a 3D active nematic, consisting of filaments (shown in green), motors and cross-linkers. Starting from the initial nematic state, the filaments first form a flat ribbon, which then undergoes a wrinkling instability.