Spontaneous circulation of active microtubules confined by optical traps

We propose an experiment to demonstrate spontaneous ordering and symmetry breaking of kinesin-driven microtubules confined to an optical trap. Calculations involving the feasibility of such an experiment are first performed which analyze the power needed to confine microtubules and address heating concerns. We then present the results of first-principles simulations of active microtubules confined in such a trap and analyze the types of motion observed by the microtubules as well as the velocity of the surrounding fluid, both near the trap and in the far-field. We find three distinct phases characterized by breaking of distinct symmetries and also analyze the power spectrum of the angular momenta of polymers to further quantify the differences between these phases. Under the correct conditions, microtubules were found to spontaneously align with one another and circle the trap in one direction.


Introduction
The study of active fluids which have a source of internal propulsion, as is common in biology, has been of interest for several decades. In more recent years, several experiments have demonstrated remarkable emergent phenomena using mixtures of microtubules and kinesin fueled by ATP. Kinesin is a motor protein that binds to a microtubule in such a way that, when powered by ATP, it tends to "walk" along the microtubule in a particular direction governed by the microtubule polarity. As the kinesin moves, it feels a drag force by the surrounding fluid. This has the effect of pushing the fluid in the direction that the kinesin is moving, as well as (by Newton's third law) exerting a force on the microtubule in the opposite direction. As such, flows can be observed in the fluid with no external impetus, from vortex lattices [1] to 2-dimensional active nematics [2][3][4][5][6] to macroscopic effects like coherent flow [7] and shear-induced gelation [8].
In what follows, we use first-principles hydrodynamic simulations to show the types of behavior (e.g., spontaneous circulation) of confined microtubules. We propose an experiment in which microtubules are held in place by an optical trap. To aid in this, we show the feasibility of such an experiment by first calculating the laser power required to contain a bending microtubule, and then we derive the change in temperature one might expect for such an experiment.
We then discuss the methods of simulation. While the specifics of the simulation are somewhat complicated, it uses no phenomenology. In fact, it has (in slightly adapted forms) been used to successfully reproduce and offer important insights into the phenomena of cytoplasmic streaming in Drosophila oocytes [9] and metachronal wave formation in microtubule bundles [10].
The results of the simulations themselves are then presented, in which we identify qualitatively different kinds of observed motion. We examine the dependence of these outcomes on parameters, and provide some interpretation. In addition to examining the motion of the polymers themselves, we also calculate and present analysis of the fluid motion in the near-and far-field. We then perform a power spectrum analysis of the angular momenta of polymers to give some concrete metrics for determining phase.

Preliminary estimations
Before going in-depth into predictions as to what will happen if microtubules are contained in an optical trap, we should first address whether it is feasible to do so, given typical issues common in optical trapping.

Laser frequency
For biological applications, lasers with wavelength >1000 nm are typical, as this reduces damage to cells [11,12]. However, even from the beginning of its usage, optical traps with wavelengths in the visible regime have been successfully used for sub-cellular structures [13].
In what follows, we consider a range of wavelengths from 500 nm to slightly over 1000 nm. These cover the most widely used lasers in optical traps, ranging from approximately 1064 nm for Yag lasers to shorter wavelength lasers frequently employed using laser diodes [14].

Forces
Microtubules confined in an optical trap would bend, and it is useful to get a sense of whether an optical trap would exert a sufficient force to keep the microtubules contained. In what follows, we use typical values to argue that an optical trap with reasonable intensity would be enough to overcome the microtubule rigidity.
The most common and straightforward approach to optical trapping uses a tightly focused Gaussian beam in the TEM 00 mode. For such a beam, the intensity of the radiation as a function of radial position can be written as where P is the power of the laser and w is the radius of the beam waist. Because the diameter of a microtubule is 24 nm and the smallest the beam diameter can be is on the order of the wavelength ∼ 500 nm, we use the electric dipole approximation where is the induced dipole of the trapped particle. For a sphere of radius a, where m ≡ n 1 ∕n 0 ; n 0 and n 1 being the refractive indices of the surroundings and the sphere, respectively [15]. Because, for a monochromatic wave, I = c 0 n 0 2 E 2 , we rewrite 2 as In terms of the Gaussian beam described by Eq. 1, For a bending microtubule, the elastic restoring force per unit length will be If the microtubule is circling the optical trap, then this becomes where r is the radius of circulation. If we now model the microtubule as a chain of beads (each bead having radius a), then the elastic force on a single bead would be Equating Eqs. 8 and 5 and solving for P gives the power required to contain a microtubule circling at radius r, This is minimized at r = w , meaning the minimum power required to contain a microtubule is The microtubule stiffness C has be measured to be approximately 1 × 10 −23 Nm 2 [16,17]. The refractive index of water is n 0 = 1.33 , and the refractive index of tubulin has been measured to be n 1 ≈ 2.5 [18,19]. Modeling the microtubule as a string of beads of radius a means a ≈ 12 nm. Using these quantities, we find that P min ≈ 16 W. This is quite large compared to standard lasers used for optical traps, but it should be emphasized that when a group of microtubules are circling an optical trap, each individual microtubule often has a radius of curvature larger than the radius of the trap if the length of the microtubule is on the order of, or less than, the trap radius (see images in Section 4). The minimum power input P min ∝ b −3 where b is the polymer radius of curvature, so even a factor of 2 increase (i.e., let b = 2r ) reduces P min to only 2 W. Also we note that laser power levels in optical traps experiments range from a few milli-Watts to a Watt or more [14], so the current proposed power is not at levels too different from previous work.

Heating
The relatively large laser power required to perform this proposed experiment raises some concerns about heating in the system, and it is worthwhile to address the degree of heating one might expect so that an experiment may be properly designed. It has been shown [20] that for a Gaussian beam, the change in temperature at the center of the optical trap of uniform absorbance is approximately where is the absorption coefficient, P is the laser power, is the thermal conductivity, is the laser wavelength, and R is a characteristic distance ( R ≫ ) to a boundary at which temperature is held constant, often taken to be the distance to the glass slide in experiments. If the experiment is performed on water, it should be noted that is highly dependent on , but this dependence has been thoroughly studied [21]. Figure 1 shows ΔT(r = 0)∕P as a function of for water, letting R = 10 m. This was computed using Eq. 11, where data for as a function of was obtained from Ref. [21] Fig. 16 (Media 1). From this, we can see that, for a 20 W laser, wavelengths longer than  Fig. 16 (Media 1). The trap is assumed to only contain water a distance of 10 m from a thermally conducting plate. The lower bound of the plot terminates at 600 nm because at this point the absorption coefficient of water becomes similar to that of glass, and the constant temperature boundary condition approximation becomes questionable ∼700 nm quickly become unfeasible, as this would lead to temperature increases of 20 K or more. However, wavelengths shorter than 700 nm would likely only heat by 5-10 K.
Adding microtubules to the optical trap would increase the temperature further, as the absorbance of proteins tends to be several orders of magnitude larger than that of water. However, an important implication of the calculation in the previous section is that the trap radius has no effect on the power required to confine the microtubules. This means that, if needed, the radius could be increased (and microtubules made longer) in order to reduce the relative area of the microtubules. For example, a 1 m long, 24 nm diameter microtubule takes nearly 1% of the area of a trap of radius 1 m. If the length of the microtubule and the radius of the trap are both increased to 10 m, the area fraction per microtubule reduces by a factor of 10 to under 0.1%.
Furthermore, the assumption that the optical trap is far from the plates will not necessarily be true, and therefore, the calculations above should be seen as a upper bound on heating. Two plates separated by a very thin gap (on the order of or even less than the radius of the trap) may be used, as the hydrodynamics of this are accounted for in the simulations that follow. This would further increase the ability of the slides to dissipate heat, especially if the slide material is chosen to have high thermal conductivity and low absorbance at the desired wavelength. Sapphire substrates have been used for this purpose as it has ∼20 times the thermal conductivity than borosilicate glass [22], but quartz (∼3 times the thermal conductivity of water) would likely also be a viable choice.
We also calculated the Rayleigh number [23] for this system with a plate separation of 10 m and conservatively taking the temperature difference to be 100 K. We find that the Rayleigh number is 0.0014, many orders of magnitude below the onset of the Rayleigh-Bénard instability. Therefore convection does not take place.
It should also be noted that when plate separation is at the nanometer scale, hydrophobic dewetting can have a substantial effect on the interface between glass and water. For example, nanoscale pores can transition to a dry state [24]. However, when the plate separation is in the micron range, as we are considering here, these effects should not be noticeable, and indeed similar geometries are employed for experimental work on active matter [6] closely related to the ones proposed here.

Method of simulation
This simulation was adapted from previous work used to successfully model cytoplasmic streaming and metachronal wave formation [9,10], and more detailed explanations of what follows can be found in these papers. This is a first-principles simulation and at a high level utilizes a straightforward approach. Microtubules are expressed as polymers, each made of a chain of N monomers bound together by a spring force which separates them by a distance = 1 . A group of M polymers is initialized, and various forces act on them, see These constants are the same as was used in previous work [10]: k kin has arbitrary units and sets the scale for other forces in the simulation, and dt was chosen to optimize simulation efficiency (verified to produce the same qualitative results as smaller time steps).
The polymers are initialized in the trap in a zigzag pattern, alternating polarity, with random noise given to each monomer's initial position. This is more computationally efficient than true Monte Carlo initialization, and we observed no noticeable difference in simulation outcome.

Kinesin drag force
The drag force propelling the microtubules depends on the concentration of kinesin. Kinesin binds preferentially to microtubules, meaning that a relatively low concentration of kinesin in the solution will result in a high concentration of kinesin on the microtubules. Suppose that these kinesins are modeled as a linear train of spheres (radius a) separated by distance d and moving at speed v 0 . It has been shown [9] that, far from the kinesin, the fluid velocity is the same as that due to thin cylinder moving at speed (a∕d)v 0 . In terms of Eq. 12, k kin ∝ a∕d.
The above demonstrates that there are two main ways of tuning k kin experimentally: one can change the kinesin concentration (effectively changing d), or one can add cargo to the kinesin (changing a). Explicitly adding cargo for kinesin to transport is not strictly necessary-many studies have shown active matter phenomena using no added cargo (although it is possible the kinesin are transporting segments of microtubule). This being said, the size of the dragged cargo is not particularly important in itself: a well-known result of slender-body theory is that the drag force F on a cylinder moving parallel to its axis is where is the dynamic viscosity, u is the cylinder speed (relative to the far-field fluid), and ℓ the cylinder length. This is onlyis the cylinder length. This is only logarithmically dependent on the cylinder radius a. Much more important is the ratio a/d, as u ∝ a∕d for the kinesin train. Fig. 2 Illustrations of some of the important forces applied to simulated polymers. (a) shows stiffness and kinesin drag forces that a polymer experiences regardless of the existence of other polymers (part of the F terms in Eq. 14), and (b) gives an example of how polymers exert hydrodynamic forces on one another via the interaction tensor

Description of other forces
The forces that go into computing ( i ) in Eq. 12 are: is the drag experienced by a small sphere, is the hydrodynamic interaction tensor (described in the following section), and F is the sum of all forces on a monomer not related to kinesin drag or hydrodynamic interactions: where: is the spring force keeping monomer separation approximately constant. For our simulations, k spr = 100 and ℓ = 1.
j = k stiff 2 j − j+2 − j−2 is the stiffness force which resists polymer bending. This is equivalent to C in Eq. 6. k stiff is varied in our simulations, with 0.02 ≤ k stiff ≤ 10.
j is the trapping force which pushes all monomers radially inward. The r 7 dependence of the trap was used for computational efficiency as well as to let polymers travel freely within the trap while providing a firm boundary at the trap radius r = R. A spring force of the form j = −k trap j was also attempted, and similar (albeit somewhat less stable) types of behavior were observed. This was not used for analysis, however, because the trap radius is less well defined and the time required for polymers to exhibit collective behavior is significantly longer. k trap = 1.0 for all simulations, and the trap radius R was varied between 2.5 and 10.0.

Interaction tensor
The hydrodynamic interaction tensor is the same as that which was used previously to simulate metachronal wave formation [10], so we will only describe its significance at a high level here. As this is Stokes flow (Re = 0), the sum of non-hydrodynamic forces exerted on each monomer is transferred perfectly to the surrounding fluid. If the monomer is sufficiently small relative to the interaction distances, such a force can be modeled as a point force (known as a stokeslet). The exact solution for the velocity field v(r) due to a free stokeslet as derived by Oseen has been known for nearly a century: is known as the Oseen tensor. The interaction tensor used in this study is a simplified version of the solution derived by Liron and Mochon for a stokeslet between two infinite flat parallel plates [10,25]. We assume that microtubules in the proposed optical trapping experiment would be confined approximately halfway between two glass slides (separated by distance H = 1.0). For this reason and for computational efficiency, we confine all monomers to this plane in all simulations.

Results
In what follows, we describe the sorts of behaviors that emerge depending on input parameters. A link to videos of simulated behavior is provided in Appendix 1. We then analyze the velocity in the surrounding fluid to give a sense of the mixing ability of each type of behavior.

Types of microtubule motion
While the motion of simulated polymers was often complex, we classify behavior into three categories: uncorrelated circulation, correlated circulation, and stasis. These are broad classifications: the descriptions provided are qualitative, and many parameter choices exhibit intermediate behavior. Below is a description of these categories, and Appendix 2 gives further information regarding how certain parameters affect the behavior exhibited by the system.

Uncorrelated circulation
If polymers do not interact sufficiently (i.e., if k Oseen is small) or polymer density is low, then each polymer tends to act independently of other polymers. As such, there is no symmetry breaking: approximately the same number of polymers circulate clockwise as counter-clockwise. Correspondingly, the polarity of the microtubules is mixed.

Correlated circulation
If k Oseen is increased and the density of polymers is sufficient, then correlated circulation is often seen: that is, polymers interact strongly enough and the stiffness is low enough such that some polymers reverse direction, breaking symmetry to exhibit organized circulation in a single direction. The polarity of the microtubules are all in the same direction, consistent with the circulation of the system.

Stasis
If polymers interact strongly but the polymer length is short compared to the trap diameter, the stiffness is too low, or the polymer density is too high, then no circulation occurs. Polymers interact but are unable to change direction, resulting in a cluster of polymers that remains largely stationary, with occasional and irregular changes in direction. Here static rotational symmetry is broken, with the microtubules clustered in one region of the trap away from the middle. It should be noted that the "Stasis" designation refers only to overall motion of the microtubule cluster. There are occasional microtubules that break away from this cluster, and the far-field fluid speeds in this phase are indeed much larger than those of the circulation phases (see Fig. 4). Figure 3 shows the direction of the velocity field for the three types of motion discussed in the previous section. The polymers are shown as black curves, with black circles indicating the negative (-) poles of the microtubules (i.e., the direction of microtubule propulsion, as motor proteins move toward the positive end). In the videos, the polymers are yellow with the kinesin motion propelling the polymers in the opposite direction to their motion. The initial conditions are discussed at the end of Sect. 3). We immediately notice that the velocity field for circulation behaviors is irregular, while that of stasis resembles the flow field for a stokeslet. This suggests that the far-field behavior for polymers in stasis will be stronger than those in circulation. Fig. 4 shows the average fluid speed far from the trap for the same simulations as Fig. 3. Indeed, this is what we see: in the far field, fluid speeds due to polymers in stasis are an order of magnitude higher than polymers in uncorrelated circulation, and nearly two orders of magnitude higher than polymers in correlated circulation. Furthermore, we verify that fluid speeds for all cases drop off as 1/r 2 . This is not entirely unexpected, as this is the far-field behavior of the Liron-Mochon interaction tensor [10,25]

Angular momentum power spectra
One method to quantify the differences between these types of behavior is to consider the power spectrum of the angular momentum. We first calculate the angular momentum of the jth monomer of the ith polymer, Note that, due to our 2D geometry, the cross-product is always in the z-direction and can as such be treated as a scalar. We then sum these to find the angular momentum L i of the ith polymer, We then find the power spectrum for each polymer over time and sum these, i.e., where L i ( ) = F[L i (t)] is the Fourier transform of Eq. 19. Figure 5(a) and (b) shows this power spectrum for the uncorrelated circulation and stasis phases. The correlated circulation phase is not shown as the power spectrum is very strongly peaked at frequency = 0 : that is, all polymers are more or less locked into stable circulation, and there is very little change in the angular momentum.
We immediately notice that the uncorrelated circulation state tends to have a single peak for > 0 . This makes sense, as any two given polymers circulating in opposite directions tend to collide twice per rotation, i.e., ≈ 2 T , where T is the period of circulation. The power spectra in Fig. 5 were over a span of 1000 time steps (giving the x-axis units of Δ = 1 1000 ), and the period of circulation for a single polymer was observed to be T ≈ 200 time steps. Indeed, we see that peak ≈ 2 T = 1 100 = 10Δ .  Figure 5(c) shows these same spectra on a log-log scale, but each averaged over its four initializations. The intent of this plot is to show that the power spectrum for the uncorrelated circulation phase has a much larger high-frequency tail compared to that of the stasis phase. Indeed, if we assume the power spectrum has the approximate form then this figure shows that ≈ 3∕2 and 1/2 for the stasis phase and the uncorrelated circulation phase, respectively.  Fig. 3(a) and (c), respectively. Each plot shows the power spectrum for four system initializations. Units on both axes are arbitrary. The power spectrum for correlated circulation is not shown, as it is so heavily peaked at zero frequency. Plot (c) averages the initializations from plot the other plot on a log/log scale. The dashed lines are a guide to the eye for |L| 2 ∼ −3∕2 (blue) and |L| 2 ∼ −1∕2 (red) In summary, this power spectrum is a very useful tool for determining type of motion: -For the correlated circulation phase, |L| 2 is essentially a delta function at = 0.
-For the uncorrelated circulation phase, |L| 2 has a distinctive peak at ≈ 2 T , where T is the period of circulation. Additionally, for high frequencies, |L| 2 ∼ −1∕2 .

Conclusion
We have shown that, under the right conditions, interesting and unique motions can occur for confined microtubules. While the design of such an optical trap experiment poses some challenges (a powerful visible-frequency laser would likely be required to contain microtubules while mitigating temperature increase), these are not prohibitive. In fact, we have shown that the trap radius should have no effect on the amount of power required to confine the microtubules, giving a fair amount of leniency in experimental design. Three distinct types of polymer motion were identified and analyzed using first-principles simulations, with insights presented as to what parameter regimes might lead to preference of one type of motion over another. We also calculated velocity fields in the vicinity of the traps and in the far-field. These calculations are important for any future applications in mixing, as it shows that fluid motion is far more localized for correlated circular motion than for other phases (although v ∝ 1∕r 2 for all types of motion).
From an experimental perspective, it may seem as though many of the parameters varied in the simulations are not tunable. For instance, one cannot substantially vary polymer stiffness when dealing with real microtubules. However, it is possible to vary other experimental parameters to achieve the same effects. For instance, the stasis phase would more likely be encouraged with: (1) a very thin system, as this both increases the strength of hydrodynamic interactions and makes the microtubules less prone to sliding over one another, when they have crossed; (2) with more kinesin (or cargo for kinesin added) to increase viscous drag; and/ or (3) increased fluid viscosity. Any increase in hydrodynamic interactions would, in effect, be equivalent to lowering the microtubule stiffness. In fact, this is one of the reasons why it is believed that Drosophila oocytes transition from the slow to fast streaming phase. Experiments [26] show that loosening the actin network, and hence lowering the viscosity, causes premature fast streaming. The slow streaming phase shows some similarity to the stasis phase described above (characterized by slow uncorrelated motions where the microtubules appear disordered), and fast streaming resembles correlated circulation.
Although the hydrodynamics of the boundary would be different, a related system would be to place microtubules inside a hard-walled cylinder instead of an optical trap. Because of the strong hydrodynamic screening induced by the plates, the types of behavior seen in a hardwalled experiment may be similar to that of an optical trap if the height of the cylinder is much less than its radius. It might therefore be of interest to try to confine microtubules this way as well, perhaps by forming them in situ inside of a thin circular boundary (with radius on the order of the microtubule length) that is then confined between two plates. In addition to having fewer design obstacles, such an experiment would also perhaps be more representative of intracellular systems.
The experimental confirmation of this effect would have important implications. For instance, the core forces and geometry of this work are very similar to those present in active nematics [3,6], so understanding this behavior would be highly relevant to such systems. Such experiments could also provide a logical next step to applications in localized mixing in microfluidics.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.