Results are presented here in terms of the normalised attachment rate kA/kD, the normalised motor rate kMτ
b
where τ
b
= Lbγ/4kBT is the approximate time for a filament to freely diffuse one monomer distance; the normalised pressure P/P0 with P0 = ε/b3 with ε = 5 kBT the Lennard-Jones repulsion energy, and, where relevant, the scaled end-detach rate kE/kD.
3.1 Stationary states
For the parameter space considered, we observe four classes of steady-state configuration as shown in Figure 2. For a low density of fast motors, spindles are observed as in Figure 2(a), which crossover to a radially-symmetric aster as the motor density is increased as show in Figure 2(b). For slower motors, we observe a nematic at low motor densities and semi-asters at high motor densities, as shown in Figures 2(c) and 2(d) resp. Semi-asters typically arise for higher pressures than asters and are more compressed, consistent with the emulsion experiments of Pinot et al. [11] and justifying our use of the term. Movies demonstrating the spontaneous emergence of all of these states from the initial conditions for exactly the same parameters are available as Additional File 1 (spindle, corresponding to Figure 2(a)), Additional File 2 (aster, corresponding to Figure 2(b)), Additional File 3 (nematic, corresponding to Figure 2(c)) and Additional File 4 (semi-aster, corresponding to Figure 2(d)).
To quantify to which state a system belongs, each filament's polarity vector is projected onto the x-y plane to give a two-dimensional unit vector aligned towards the [+]-end. This is averaged over all filaments whose centers of mass have azimuthal angle θ with respect to the center of the system, giving rise to the mean orientation . This is then decomposed into angular mode vectors a
m
and b
m
,
(1)
from which can be defined the mode amplitudes Q
m
,
The Q
m
are invariant under global rotations of the whole box.
To determine the corresponding state, the measured Qm up to m = 3 are compared to known values for ideal states, and that with the closest Euclidean distance is taken to be the state. The values for pure asters and nematic phases are easy to derive; for an aster , (Q0, Q1, Q2, Q3) = (0,1,0,0), whereas for the nematic state where , (Q0, Q1, Q2, Q3) = (0,0,0,0). Note that while an isotropic state would give the same Qm as for the nematic, such states only arise for kA and P well below the considered ranges. For spindles and semi-aster states there is a degree of choice in how the target Qm are calculated, so we choose simple forms that permit exact evaluation of the Q
m
. For the spindle, for θ ∈ (-π/4, π/4) or θ ∈ (3π/4, 5π/4) and zero otherwise, for which (Q0, Q1, Q2, Q3) = (0, 2/π2 + 1/2, 0, 2/π2). For the semi-aster, for θ ∈ (-3π/4, 3π/4) and zero otherwise, for which (Q0,Q1,Q2,Q3) = (9/π2, 9/4π2, 333/1715π2, 9/100π2). Variations in these forms have been tested and although the boundaries between states shift slightly, the underlying trends remain the same.
The occurrence of the four steady-states, plus a fifth 'vortex' state to be discussed below, with motor density and speed are presented in Figure 3 for two different pressures. The observed configurations correlate with the density and distribution of motors along the filaments. The mean number of motors per filament nmot/N is plotted in Figure 4 and shows an expected increase with the attachment rate kA as well as the pressure. The increase with pressure can be understood as due to the closer packing of the filaments, increasing the number of potential attachment points for motors and hence nmot. The approximate scaling is faster than the linear relationship measured for constant volume, two-dimensional simulations [19], presumably due to similar reasons: As kA increases so does the motor density which, in this constant-pressure ensemble, allows the system to contract, presenting more potential attachment points between monomers and hence further increasing the motor density. A derivation of the value 3/2 of the exponent is not available so far.
Motors move to the [+]-end and dwell there until detaching, thus a greater fraction are expected to occupy filament [+]-ends, potentially resulting in tight binding mediated via many motors. Plotted in Figure 5 is the fraction of motors with at least one head at a filament's [+]-end, , for the same runs as in Figure 4. By comparing to the configuration plots in Figure 3 it is possible to infer signatures of crossovers between states in the inflection points in these curves. For the slower motors, there is an increase in as the nematic state changes to a state with a greater degree of polar ordering (spindle or semi-aster depending on the pressure). For the faster motors there is a marked increase in with kA, which corresponds to the crossover to the aster state with a high degree of [+]-end binding. Further indication of the importance of [+]-end binding is presented in section 3.3 where enhanced end-unbinding rates kE > k
D
are considered.
3.2 Dynamics and vortices
The stationary states described above admit no spontaneous non-equilibrium flows, despite the motor motion generating a positive energy flux: The increase in the stored motor elastic energy due to motor motion and thermal drift of the connected filaments is balanced by the loss due to detachment, with no observed net translocation or rotation of the filaments in steady-state. Collective rotation of all filaments about a fixed center arises for one region of the considered parameter ranges, but appears to be a transient flow that irreversibly contracts to a non-rotating semi-aster configuration. These states are referred to here as vortices due to their superficial similarity with the rotational modes observed in microtubule-kinesin assays [15, 16] and are described in detail in this section. A snapshot is given in Figure 6 and movies are provided as Additional File 5 (all filaments shown) and Additional File 6 (one filament highlighted for the same run as Additional File 5).
Collective rotation of the whole system can be quantified by the mean angular velocity of filament centre-of-mass vectors r relative to the system center, or alternatively by the net transverse velocity of each filament's centre-of-mass relative to its polarity, . Here we employ the latter as it is available for all of our runs, but we have confirmed that it closely tracks the angular velocity in those runs for which both were measured. Examples of for 3 independent runs are given in Figure 7, and show finite rotation of either sign until the system irreversibly contracts to a semi-aster state and rotation ceases. This contraction time can be confirmed by visual observation of system states, and can be precisely located by fitting the system radius as a function of time, R(t), to the four-parameter hyperbolic tangent R(t) = Rmin + ΔR tanh[(t - tcont)/Δt]. The mean of is presented in Figure 7 for each run as a horizontal line segment, that extends from t = 0 to the contraction time tcont found from this fit. In all cases, tcont coincides with the rapid decay of to zero, providing independent confirmation that rotation ceases when the system contracts to a semi-aster.
It is now possible to define a vorticity order parameter for each point in parameter space. For each run α, the mean μ
a
and standard deviation σ
α
of is calculated starting from t = 0 up to the time that the system contracts. This is regarded as significant if the mean is comparable to or larger than the standard deviation, but since the sign is arbitrary we also take the absolute value to give the vorticity for a single run,
This is then averaged over all runs with the same parameters to give the mean vorticity . A given point in parameter space is then regarded as exhibiting a (transient) vortex if exceeds some arbitrarily-chosen value of order unity. The corresponding region of parameter space for is plotted in Figure 3 and arises for higher densities of motors that are not so fast that they aggregate at [+]-ends, which would stabilize an aster relative to a vortex. Independently varying the fraction of [+]-ended motors by increasing the end-detachment rate kE supports the existence of a critical fraction for vortex formation, as discussed in Sec. 3.3.
The reciprocal relationship between vorticity and contraction time is clearly evident when both quantities are plotted together; see Figure 8. Stronger vortices have a shorter lifetime than weaker vortices. The distribution of contraction times is presented in Figure 9(a) for a single pressure. As the data is noisy we do not attempt to extract an arbitrary distribution, but instead make two statistical tests for the most basic possibilities, i.e. that contraction happens at a fixed time, which would give a Gaussian distribution, or that it happens at a fixed rate corresponding to an exponential distribution. To give some measure of the goodness-of-fit, the Anderson-Darling statistics for an exponential distribution with an unknown mean has a significance level of P ≈ 0.2, whereas that for a normal distribution of unknown mean and variance has a significance level of P ≈ 0.025 [31, 32]. This clearly favours the exponential over the Gaussian distribution. Attaining even this noisy data consumed considerable computing resources and we were unable to repeat this procedure for other parameter values.
Assuming the true distribution is exponential, this would suggest that contraction is triggered by spontaneous fluctuations that occur at a constant rate in time. From observation of movies of filament arrangements, a likely candidate is the transient void formation frequently observed near the outer wall, where nearby filaments are attached purely by motors at their [+]-ends and not along their length. Such voids, when large enough, lead to a 'hinge'-like mechanism in which the void expands and one section of the polarity field inverts, leading to the semi-aster.
The onset of vorticity is also evident in the histogram of the signed vorticity, i.e. the before taking the modulus in eqn. (3), which can be positive or negative depending on the direction of rotation. For low pressures with this distribution is unimodal around the origin, but becomes bimodal when vorticity is more evident as demonstrated in Figure 9(b). Of the 20 runs presented here, 12 rotated in one direction and 8 in the other, which has a significance interval of P ≈ 0.5 as determined from a Binomial test with equal probabilities for both directions. This is to be expected given our use of stochastic initial condition that does not predispose the system to any preferred rotational direction.
Independent confirmation of vorticity can be inferred from the mean-squared angular deviation 〈(Δθ)2〉 already defined in Sec. 2. This is plotted in Figure 10 for the same parameters as above as a function of the lag time Δt, averaged over all waiting times t up until tcont. There is a crossover from linear behavior (Δθ)2 ~ Δt for low pressures with low vorticity, to a more rapid scaling (Δθ)2 ~ (Δt)2 for pressures well into the vortex regime. Since this quantity is the angular analogue of the mean squared displacement for translation degrees of freedom, these two limits can be regarded as diffusive and ballistic, respectively. Microscopically the diffusive limit corresponds to fluctuations with no net drift, whereas the ballistic limit arises when all filaments are rotating around the system with a constant angular velocity in the same direction. Therefore the vortex state should correlate with ballistic motion, and comparison of Figures 8 and 10 confirms this. This figure also demonstrates that the integrated angular rotation of the vortex before contraction is typically larger than π/2, much larger than the diffusive drift ≈ π/10 over the same time frame.
3.3 Enhanced detachment from ends
In the simulations that accompanied the microtubule experiments, it was claimed that the residence time at the microtubule [+]-ends played a crucial role in determining the vortex stability, with an enhanced end-detachment rate required to form vortices [16]. By contrast, for our model it is the fraction of motors at filament [+]-ends that determines vortex stability relative to an aster. The vorticity, motor density and fraction of motors at [+]-ends are plotted in Figure 11 against end-detachment rates kE ≥ kD for two sets of kA, kM and P. It is clear that increasing kE can both destroy a vortex that existed when kE = kD, and create a vortex when kE = kD gave an aster. In order of increasing kE, the sequence aster → vortex → semi-aster (where any vortex is either absent or too short lived to be discerned) is typically observed, although we do not claim this sequence is followed by all points in parameter space. There is a slight increase in motor density in the semi-aster state as evident from the figure, resulting from an increase in potential attachment points due to the increased density.
Thus residence time at [+]-ends, which is ∝ , is not the determining factor with regards vortex stability here. Rather, vortices coincide with around 25% of motors at [+]-ends as highlighted in the figure. There is no critical dependency on the motor density, although we speculate that below some minimum value spindles or nematic states would be observed instead. The critical fraction 25% will likely depend on parameters that were not varied in this work, such as filament length M and the motor spring stiffness. A systematic survey of these parameters, or of kE ≠ kD for all kA, kM and P, is however beyond the scope of this work.
3.4 Controlled volume
One message from the previous sections is that the observed steady-state is predominately determined by the density of motors and the fraction at [+]-ends. It may appear that the primary role of motor motion, which would be the source of any non-equilibrium effects in this model, is merely to select the fraction of [+]-ended motors, faster motors giving a higher fraction. It might even be speculated that even the transient vortex state is driven, not by motor motion, but rather as a protracted buckling event powered by the pressurised walls.
It is straightforward to show that motor motion can drive vortex motion, however. Plotted in Figure 12 is the rotational velocity for two independent runs in a box with fixed radius, where v is a characteristic filament velocity. The system is initially in an aster configuration, but when the radius is suddenly reduced by b/2 at a time t/τ
b
≈ 1.6 × 104, the system switches to a rotating vortex state that appears to be long-lived; the total time window in this figure is an order of magnitude longer than the longest vortex described in Sec. 3.2 (which has the same parameters). Since there is no energy input from the walls, the only possible cause for this rotation is the motor motion. Thus the pressure ensemble is important to let the system adjust its density to the vortex state; however, the same pressure also destabilizes the vortex state, because it favors further contraction into the semi-aster.
Although the magnitude of the rotational velocity remains fixed (note that the characteristic velocity v is the same for both runs and constant in time), the direction aperiodically reverses as evident in the changes of sign in the figure. The statistics of time intervals between direction switching suggest that the underlying mechanism may be the same as for the contraction to the semi-aster state in the constant pressure case. Specifically, the mean switching time Δtswitch/τ
b
= 21.6 × 104 ± 5.8 × 104 is consistent with the mean contraction time ≈ 1.7 × 104 measured earlier, and again is consistent with an exponential distribution (significance level P ≈ 0.2 from n = 7 values using the Anderson-Darling test [31, 32]). This suggests that the same spontaneous fluctuation that permits contraction under constant pressure instead promotes rotational reversal under constant volume.