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 self-organize 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 three-dimensional 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 three-dimensional 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 . 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 . 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 director[1] (Figure A). 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 B–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 E). 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.

Figure 1

3D wrinkling formation. (A) Schematic representation of the wrinkling process as a result of contractile and extensile forces that emerge at different times. The 3D active nematic system evolves to exhibit contraction into a sheet, a wrinkling instability simultaneously as a lateral contraction, and a final transition to active turbulence with spatiotemporal chaos. The red arrow indicates the nematic director before the onset of wrinkling instability. The black arrows indicate the elastic deformations induced by the active stress along different directions. (B)–(E) Time evolution of 3D active nematics through the different stages. The nematic order along the x-direction is visible. Scale bar: 500 μm. (F), (G) Pattern visualization at two different heights. The 3D structure can be seen by the wave crests in focus and the troughs out-of-focus when the focal plane is set to z = 0 μm. At z = 70 μm the situation is reversed (see Movie S2). The experiment was conducted at 2 mM ATP. Scale bar: 250 μm.

Figure 2

Reconstruction of the nematic director field. This confirms the nematic order of the microtubule bundles in the samples. The image is from Figure C. For details see the Supporting Information.

Contraction Characterization

The evolution of the system from 3D to 2D includes contraction along both y- and 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 H,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 A,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).

Figure 5

Simulation of the active nematic using the Cytosim package. It reproduces the dynamics observed in the experiment (see also Movie S4). (A) Initially, the 3D volume is filled with nematically ordered filaments (with orientations within θmax = 30° of the x-axis), mixed with tetrameric active motors and passive cross-linkers. Filaments are subject to repulsive hard-core interactions and attractive depletion forces. (B) The filaments first form a ribbon that contracts laterally. (C) Due to extensile stress, the ribbon forms wrinkles. At the same time, increased bundling of filaments is visible. (D) The wrinkles grow until they get constrained by the top and bottom walls of the channel. (E) Under constraint, the bundles increasingly extend in horizontal direction and the ribbon structure is dissolved. (F) Finally, a state of 3D active turbulence is reached. Upper row: top view. Lower row: side view. (G) Average stress −σ in a cross-section of the channel as a function of simulation time. Positive values represent extensile stress. The letters indicate the times of snapshots in (A)–(E). (H), (I) Kymographs showing the average filament density across the width (H) and height (I) of the channel as a function of time. Contraction in both directions is maximal at the onset of the wrinkling instability.

Figure 3

Contraction with and without depletion agent PEG. (A), (B) Lateral contraction over time of a network made of microtubule and PEG mixture without addition of motor proteins. The system passively contracts due to the depletion effect of PEG. (C), (D) Lateral contraction and longitudinal expansion of a microtubule–motor proteins network over time without addition of PEG. The contracting effect due to the motors proves that the network exhibits an active contraction contribution. The experiment was conducted at 2 mM ATP concentration. Scale bar: 500 μm.

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 C,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 D). We find that the wavelength of the pattern and the time scale of its formation are correlated; 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 cross-linking 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 E). 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.

Figure 4

Different ATP and motors concentration. (A)–(C) Wrinkle formation at different ATP concentrations: 2 mM; 1 mM; limiting concentration 10 μM. The concentration of motor proteins is 17 nM and constant in the three experiments. (D) Wrinkling wavelength, plotted against the appearance time. The continuous line shows the dependence λ ∝ τw1/4, predicted from the model for samples that differ only in the amount of active stress. (E) When the kinesin concentration is reduced by 10-fold, a strong lateral collapse is observed. (F) In some samples, pattern formation can be observed in both x-direction and y-direction. Scale bars: 500 μm

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̂ ≃ ê. 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 ê and ê). After the contractile stress in the z-direction 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 traceless tensor , with the order parameter S, which has the value 1 for perfect nematic alignment and 0 for an isotropic state. The active stress that is caused by the motors can be written as[2,31,35]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 aswhere 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 Γ∂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Λei for a perturbation with the wave vector q. We find a growth rate ofFor 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 and a growth rateAssuming 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

The experimental data shown in Figure D confirm the trend predicted by eq , 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 F. 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 , Figure S6, and Movie S4.

Starting from a nematically ordered state (Figure A), the simulation first shows the condensation in the middle of the channel and shrinking in the lateral direction (Figure B). As in the experiment, formation of microtubule bundles also becomes visible at this stage. Afterward, the extensile stress (Figure G) leads to a wrinkling instability in the sheet (Figure C). The growth of the selected out-of-plane deformation mode continues until the amplitude becomes limited by the simulation box size (Figure D). Afterward, the ribbon structure disintegrates and the system ends in a state of 3D active turbulence (Figure E,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 H,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 cross-linking 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 F,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).

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).