format: text

text:
# Modeling oscillatory Microtubule–Polymerization

## I Introduction

Microtubules are cylindric filaments that are used in cells for many different purposes, being vitally involved in cell motility and division, in organelle transport, and in cell morphogenesis and organization . The precise ways in which microtubules achieve their amazing variety of cellular functions is not fully understood yet. Microtubules in cells are generally dynamic, they assemble, disassemble or rearrange on a time scale of minutes. GTP (guanosine triphosphate) hydrolysis is apparently the driving force of microtubule physiology.

The rich non–equilibrium dynamics of microtubules, including the nucleation and polymerization kinetics etc.  is attracting considerable attention, both experimentally and theoretically . Two phenomena in this area, the dynamical instability of microtubules  and the oscillatory polymerization  challenge theoretical modeling already for a while .

Oscillations during microtubule polymerization have been observed either when GTP is regenerated enzymatically from endogenous GDP (guanosine diphosphate)  or when some amount of GTP is provided at the beginning or during an experiment. In the latter case oscillations occur only as a transient, because GTP is either consumed or some reactions steps may be inhibited due to the accumulation of GDP . If both possibilities are combined, the length of a transient regime depends on the initial concentrations of GTP and GDP and on the capacity to regenerate GTP. Present models for microtubule polymerization focus mainly on a description of transiently occurring oscillations and the solutions of the respective models are mostly numerical .

In recent in vitro experiments, however, the capacity to regenerate GTP has been enhanced and extended up to several hours . Compared to a typical oscillation period during microtubule polymerization, which is of the order of a minute, the reaction conditions in these experiments are almost quasi–stationary over a long range of time. Therefore we focus on modeling microtubule polymerization for time–independent regeneration conditions. As starting point we take common reduced models, where several elementary processes of the real biochemical reaction are described by a few effective reaction steps as explained in Sec. II, cf. Refs. . Reductions of complex chemical reaction schemes are quite common and a famous example is the so–called oregonator  that is a reduced model for the legendary Belousov–Zhabotinsky (BZ) reaction . However, since microtubules are long filaments, there are essential differences between the polymerization of microtubules filaments and common chemical reactions. For instance microtubules may undergo an orientational ordering transition beyond a critical filament density , a phenomenon, which doesn’t occur in common chemical reactions. Accordingly, the length distribution of microtubules is explicitly taken into account for all variants of models investigated in this work. Such models are the basis of future work on interesting pattern–formation phenomena related to the interplay between orientational ordering of the filaments and the kinetics involved in the filament growth .

In addition, we focus on model variants that include the possibility of an oscillatory microtubule polymerization and that allow analytical approaches. However, the reaction steps, such as nucleation, growth and decay of microtubules or the rate limiting factors of oligomer decay or tubulin regeneration, which have been identified to be crucial for oscillations , are taken into account. Moreover we address the question whether microtubule oscillations occur transiently or in a persistent manner beyond a Hopf bifurcation. Whether such a Hopf bifurcation takes place super- or subcritically is investigated in terms of the so–called amplitude expansion.

It is not a major goal of this work to achieve quantitative agreement between the results obtained with phenomenological models and experimental measurements. However, since the present understanding of the mechanism leading to oscillatory microtubule polymerization is incomplete, reduced models may be an appropriate tool for working out typical trends that may be testable in experiments. For comparison it is very helpful that for reduced models trends may be worked out analytically and may be presented by simple formulae. A number of spatiotemporal phenomena involving microtubule polymerization call for a better understanding too , but also in this case simple and effective models are indispensable in order to keep the modeling tractable .

At the transition to oscillatory polymerization the stationary state becomes sensitive against small perturbations, which grow or decay exponentially, $`e^{\sigma t}`$. Here the exponential factor $`\sigma =\sigma _r\pm i\omega _c`$ is the sum of the so–called growth rate $`\sigma _r`$ and the oscillation frequency $`\omega _c0`$. Below the bifurcation point the growth rate $`\sigma _r<0`$ is negative and the perturbations are damped. Beyond the bifurcation point $`\sigma _r`$ is positive and the stationary polymerization state is unstable against oscillatory perturbations. Hence the Hopf bifurcation to oscillatory polymerization takes place when the real part $`\sigma _r`$ of both roots passes zero. The investigation of the polymerization dynamics beyond the Hopf bifurcation requires in most cases a numerical analysis of the basic reaction equations. However, close to threshold $`\sigma _r`$ is small and the oscillation of the polymerization, described by the real part of $`e^{i\omega _ct}`$, is much faster than the temporal evolution of the complex valued amplitude $`A(t)`$ of the oscillations. Therefore the oscillation may be written as a product of both factors, i.e. $`A(t)e^{i\omega _ct}`$, and there is a very general approach, the so–called amplitude expansion, for separating the dynamics at these two disparate time scales . The amplitude equation describing the evolution of the amplitude $`A(t)`$ is obtained by a perturbation expansion of the reaction equations with respect to the slowly varying amplitude $`A(t)`$ and it is of the form

$$\tau _0_tA=\epsilon (1+ia)Ag(1+ic)A^2A.$$
(1)
The control parameter $`\epsilon `$ measures the relative distance from the bifurcation point and $`\tau _0`$ is the relaxation time defined by $`\tau _0=\epsilon /\sigma _r`$, that depends on the system. If the coefficient $`g`$ of the nonlinear term is positive, the bifurcation to the oscillatory state is supercritical and if it is negative, the bifurcation is subcritical. The imaginary parts of the prefactors describe the linear and the nonlinear frequency dispersion. Especially about the extension Eq. (1) including spatial degrees of freedom, there exists a rich literature as summarized e.g. in a recent review . Here in this work we calculate the coefficients of the universal equation (1) for microtubule polymerization and we discuss their variation in terms of the reaction rates.

In the following section II we describe the main steps of the reaction cycle for microtubule polymerization and the respective equations for two models are presented. The time–independent solutions for the stationary polymerization are given for both models analytically in section III. Those become unstable against oscillatory perturbations in the range of high of tubulin–dimer density. The respective linear stability analysis and the derivation of the oscillation threshold are given in section IV, including their dependence on the reaction parameters. Readers who are mainly interested in numerical results about the oscillation threshold may proceed directly to section IV A 4. The partial differential equations for growing and shrinking microtubules are of first order in the length microtubules and first order in time. Their straight forward discretization and numerical solution has to be considered with care, therefore a stable numerical scheme, that becomes exact close to the threshold of the Hopf bifurcation, is described in section V. The derivation of the universal equation (1) is outlined in Sec. VI whereas the technical details are given in appendix A. With a summary and an outlook about modeling microtubule polymerization we conclude this work in Sec. VII.

## II Models for microtubule polymerization

Microtubule assembly and disassembly proceeds in several steps . Aggregation of guanosine triphosphate (GTP) liganded tubulin dimers, the so–called tubulin–t, to microtubules is started by heating up tubulin solutions to a temperature of about $`3037^o`$C in the presence of GTP. Then microtubules spontaneously nucleate and polymerize to long rigid polymers made up of $`\alpha \beta `$ tubulin dimers. An increasing number of long microtubules in a solvent causes an increasing turbidity and the amount of polymerized tubulin–t may be monitored by measuring this turbidity  or by X-ray scattering . The nucleation of microtubules is a rather complex process and it is still a matter of debate whether the nucleation rate depends in experiments only on the initial concentration of tubulin-t, $`c_t`$, or during the polymerization on the temporally varying $`c_t`$ . But once microtubules are formed, they grow and the available tubulin–t dimers will be used up. The growth velocity of microtubules, $`v_g`$, is rather sensitive to temperature variations but it is rather independent on $`c_t`$ . Growing microtubules may change their state to rapidly depolymerizing ones by the so–called catastrophe rate $`f_{cat}`$. In previous works for the catastrophe rate mostly an exponential dependence on the tubulin–t concentration was assumed, i.e. $`f_{cat}\mathrm{exp}(c_t/c_f)`$ with some constant $`c_f`$ . Once microtubules have changed from growth to shrinking, they shrink rather quickly with a large velocity $`v_sv_g`$.

During the depolymerization of microtubules they are fragmented into oligomers or directly into guanosine diphosphate (GDP) liganded tubulin dimers, the so–called tubulin-d dimers. The oligomers themselves are believed to fragment further into tubulin–d dimers and the decay rate depends on the free GTP and GDP. Oligomers are stabilized by GDP and destabilized by GTP . If an excess of GTP is available, then tubulin–d in solution will exchange its unit of GDP for GTP and each tubulin–t dimer resulting from such an exchange step is identical to the initial tubulin–t dimer. Such a regeneration step completes the whole microtubule polymerization cycle. If a continuous source of GTP is provided, for instance by a regeneration process, this cycling may be continued over a long time . The variation of the reaction rates of the polymerization cycle with the concentration $`c_t`$ may depend on the specific experiment.

There are rather detailed models available to describe this reaction cycle of microtubule polymerization, see e.g. . As a simplification of this complex biochemical reaction we only take into account as rate limiting factors one of the two intermediate steps of the polymerization cycle, either the dynamics of shrinking microtubules (upper cycle in Fig. 1) or the decay dynamics of oligomers (lower cycle in Fig. 1). Without both rate rate limiting factors there are no microtubule oscillations, but one is sufficient for oscillations. The two simplified reaction schemes, as sketched in Fig. 1, are analyzed in detail in this work.

### A Dynamics of growing microtubules

During microtubule polymerization there are many growing filaments in a unit volume and their length distribution may be described by a length and time–dependent function $`p_g(l,t)`$ whose detailed form varies with the experimental conditions. A simple model for the dynamics of distribution of growing microtubules $`p_g(l,t)`$ is described by the following first order differential equation

$`_tp_g`$ $`=`$ $`f_{cat}p_gv_g{\displaystyle \frac{p_g}{l}}.`$ (2)
$`f_{cat}`$ describes either the transition from the growing to the shrinking state of microtubules (model I) or the decay of growing microtubules into oligomers (model II). $`v_g`$ is the growth velocity of the microtubules.

#### 1 Growth velocity and catastrophe rate

In recent experiments with high tubulin–t concentration the growth velocity $`v_g`$ was rather independent of $`c_t`$ . Since we are mainly interested in the oscillatory behavior of microtubule polymerization, that occurs at high $`c_t`$–concentrations, we assume a constant $`v_g`$ in this work. In most of the present models a $`c_t`$–dependent catastrophe rate $`f_{cat}`$ is crucial for oscillatory polymerization of microtubules. Rather common is an exponential $`c_t`$–dependence

$`f_{cat}(c_t)`$ $`=`$ $`fe^{c_t/c_f},`$ (3)
with the amplitude $`f`$ and the decay constant $`c_f`$. However, also a linear $`c_t`$–dependence

$`f_{cat}(c_t)`$ $`=`$ $`\overline{f}(c_uc_t),`$ (4)
with an appropriate constant $`c_u>c_t`$ leads to an oscillatory microtubule polymerization as we show in Sec. IV A. A hyperbolic $`c_t`$–dependence of $`f_{cat}`$, as discussed in Ref. , also supports oscillating polymerization.

#### 2 Nucleation and boundary conditions

The nucleation process of microtubules is rather complex and it has been investigated in more detail in Refs. , recently. The nucleation rate $`\nu `$ depends on the initial concentration $`c_0`$ of tubulin dimers, but it is rather independent of the temporal variation of $`c_t`$, as observed in recent experiments . Accordingly, for a given initial concentration $`c_0`$ we assume a constant nucleation rate $`\nu `$. The nucleation rate $`\nu `$ itself defines a boundary condition for the length distribution of growing microtubules $`p_g(l,t)`$ at $`l=0`$,

$$p_g(l=0,t)=\frac{\nu }{v_g}.$$
(5)
### B Model I includes the dynamics of shrinking microtubules

In model I we take into account as an intermediate step between growing microtubules and tubulin–d dimers the dynamics of shrinking microtubules, $`p_s(l,t)`$. Here the catastrophe rate $`f_{cat}`$ describes the transition of microtubules from the growing to the shrinking state. The depolymerization speed $`v_s`$ of shrinking microtubules, $`p_s(l,t)`$, is mostly much larger than the growth velocity $`v_g`$. Having microtubules in two different states, one may also expect a transition from the shrinking back to the growing state, as described by a rate $`f_{resc}`$. So one has two coupled equations for the growing and shrinking microtubules

$`_tp_g`$ $`=`$ $`f_{cat}p_g+f_{resc}p_sv_g_lp_g,`$ (7)
$`_tp_s`$ $`=`$ $`f_{cat}p_gf_{resc}p_s+v_s_lp_s.`$ (8)
The rescue rate $`f_{resc}`$, however, is usually very small in experiments and therefore it is neglected in this work. The boundary condition for shrinking microtubules is

$$p_s(l\mathrm{},t)=0,$$
(9)
because the transition from growing to shrinking microtubules is the only source for the shrinking ones and $`p_g(l\mathrm{},t)`$ vanishes for large values of $`l`$.

The temporal evolution of the concentration of the tubulin–t dimers $`c_t`$ and tubulin–d dimers $`c_d`$ is described by two equations as follows

$`_tc_t`$ $`=`$ $`\gamma v_g{\displaystyle _0^{\mathrm{}}}𝑑lp_g(l,t)+\alpha c_d,`$ (11)
$`_tc_d`$ $`=`$ $`\gamma v_s{\displaystyle _0^{\mathrm{}}}𝑑lp_s(l,t)\alpha c_d.`$ (12)
The first term in Eq. (11) describes the consumption of tubulin–t during the polymerization (growth) of microtubules and $`\gamma `$ is a length factor describing the number of tubulin dimers that are incorporated in a unit length of microtubules. $`c_t`$ is regenerated from $`c_d`$ by exchange the unit GDP for GTP and this regeneration process, described by the rate $`\alpha `$, occurs in Eq. (11) as a source and in Eq. (12) as a sink. Tubulin–d dimers are released during the depolymerization of microtubules $`p_s(l,t)`$ and this source is described by the integral in Eq. (12).

Tubulin dimers may be a constituent of growing or shrinking microtubules or they carry GTP or GDP as single dimers, but altogether they are conserved as expressed by the condition

$$c_t+c_d+\gamma L=c_0.$$
(13)
Here $`c_0`$ describes the overall concentration of tubulin dimers and $`L(t)`$ is the integrated length of all microtubules per unit volume

$$L(t)=_0^{\mathrm{}}𝑑ll\left(p_g(l,t)+p_s(l,t)\right).$$
(14)
The tubulin–d concentration $`c_d`$ may be eliminated from Eq. (11) by using the conservation condition (13). On the other hand Eq. (13) in combination with Eqs. (II B) and Eq. (11) yield an equation that is identical with Eq. (12). Hence Eqs. (7) and (8) together with

$`_tc_t`$ $`=`$ $`\gamma {\displaystyle _0^{\mathrm{}}}𝑑l\left(v_gp_g+\alpha l(p_g+p_s)\right)+\alpha (c_0c_t),`$ (15)
describe the polymerization dynamics of microtubules for model I, whereby a constant growth and shrinking velocity is assumed in this work.

#### 1 Rescaling of model I

After rescaling time $`t`$ and length $`l`$, i.e.

$`t^{}=\alpha t,l^{}={\displaystyle \frac{\alpha }{v_g}}l,`$ (16)
it is easy to see that model I may be characterized by a set of dimensionless parameters

$`\gamma {\displaystyle \frac{v_g}{\alpha }},{\displaystyle \frac{v_s}{v_g}},{\displaystyle \frac{\nu }{\alpha }},{\displaystyle \frac{c_0}{c_f}},{\displaystyle \frac{f_{resc}}{\alpha }}.`$ (17)
Some of these dimensionless quantities may be further combined to other dimensionless parameters as for instance in the threshold condition given in Sec. IV A.

#### 2 Reduced model

Since the depolymerization velocity $`v_s`$ is much larger than the growth velocity $`v_g`$ one may also consider the limit $`v_sv_g`$. In this case the shrinking microtubules decompose nearly instantaneously into tubulin–d dimers and growing microtubules decay effectively, due to the short life time of the shrinking microtubules, into tubulin–d dimers. In order to describe this direct decay the source term in Eq. (12), $`\gamma v_s_0^{\mathrm{}}𝑑lp_s(l,t)`$, must be replaced by $`\gamma f_{cat}_0^{\mathrm{}}𝑑llp_g(l,t)`$. Eliminating again the density $`c_d`$ one ends up with a reduced model for only two densities:

$`_tp_g`$ $`=`$ $`f_{cat}p_gv_g_lp_g,`$ (19)
$`_tc_t`$ $`=`$ $`\gamma {\displaystyle _0^{\mathrm{}}}𝑑l\left(v_g+\alpha l\right)p_g+\alpha (c_0c_t).`$ (20)
This simplified model reproduces essential aspects of stationary polymerization of microtubules as described in Sec. III.

### C Model II includes the dynamics of oligomers

Oligomers occur as an intermediate product during the decay of microtubules and they are made of several tubulin dimers. This intermediate product is ignored in model I. Here in model II, after the so–called catastrophe, we ignore the dynamics of shrinking microtubules as an intermediate step and instead we take into account the (decay) dynamics of oligomers. Therefore, the catastrophe rate $`f_{cat}`$ in Eq. (2) describes for model II a direct transition of growing microtubules into oligomers. Furthermore, it is assumed that oligomers decay with the rate $`\chi `$ into tubulin–d dimers. The concentration of oligomers is denoted by $`c_{oli}`$ and its dynamics as well as that of $`c_d`$ are described by the two equations

$`_tc_{oli}`$ $`=`$ $`\eta f_{cat}{\displaystyle _0^{\mathrm{}}}𝑑llp_g(l,t)\chi c_{oli},`$ (22)
$`_tc_d`$ $`=`$ $`\chi \lambda c_{oli}\alpha c_d.`$ (23)
$`\eta `$ is a measure for the number of oligomers per unit length of the microtubules and $`\lambda `$ is a measure for the number of tubulin dimers per oligomer. Oligomers decaying with the rate $`\chi `$ build a source term in the equation for tubulin–d dimers in Eq. (23).

The conservation law for the concentration of all tubulin dimers takes the form

$`c_t+c_d+\lambda c_{oli}+\eta \lambda {\displaystyle _0^{\mathrm{}}}𝑑llp_g(l,t)=c_0.`$ (24)
The equation for the growing microtubules is the same as for model I, cf. Eq. (2), but in the equation for $`c_t`$, cf. Eq. (11), one has to replace the length factor $`\gamma `$ by the product $`\eta \lambda `$. Eliminating $`c_{oli}`$ model II is described by Eq. (2) and Eq. (11) together with the following dynamical equation for $`c_d`$

$`_tc_d`$ $`=`$ $`\chi \left(c_0c_tc_d\eta \lambda {\displaystyle _0^{\mathrm{}}}𝑑llp_g\right)\alpha c_d.`$ (25)
As boundary condition for the growing microtubules we again use Eq. (5) with a constant nucleation rate $`\nu `$. For model II we only consider the catastrophe rate given in Eq. (3). This again guarantees a nonlinear feedback of the dynamics of the tubulin–t dimers to the dynamics of the growing microtubules.

#### 1 Reduced model

Similar as for model I also model II becomes in the limit $`\chi \mathrm{}`$ identical with the model described by Eqs. (II B 2). If we assume a very fast dissociation of oligomers into tubulin–d dimers, $`\chi 1`$, we can neglect the intermediate state $`c_{oli}`$. In this case the source term in equation (23), $`\chi \lambda c_{oli}`$, can be replaced by the source in equation (22), cf. $`\eta \lambda f_{cat}_0^{\mathrm{}}𝑑llp_g`$, which describes the direct decay of growing microtubules into tubulin–d dimers. After replacing $`c_d`$ and setting $`\gamma =\eta \lambda `$ in Eq. (11), we again obtain with the help of the conservation law (24) the simple reduced model as described by the equations (II B 2).

## III Stationary solutions

A polymerization cycle with a stationary length distribution of microtubules and time–independent dimer concentrations $`c_t`$, $`c_d`$ or oligomer concentration $`c_{oli}`$ are one type of the solutions of the model equations described in the previous section II. For this stationary state the various polymerization steps, such as nucleation, assembly and disassembly of microtubules as well as the regeneration of tubulin-d dimers are in a balanced state. Under certain conditions a stationary polymerization is observed in experiments . However, it may become unstable against oscillatory perturbations if the initial tubulin dimer concentration $`c_0`$ is large enough, as shown in the following section IV.

### A Model I

Eqs. (II B) are first order linear differential equations with respect to the length $`l`$ and in the stationary case these equations have exponentially decaying solutions, which take in the limit of a vanishing rescue rate the form

$`p_{g,s}^{(0)}(l)`$ $`=`$ $`{\displaystyle \frac{\nu }{v_{g,s}}}\mathrm{exp}\left({\displaystyle \frac{f_{cat}^{(0)}}{v_g}}l\right).`$ (26)
The catastrophe rate $`f_{\mathrm{𝑐𝑎𝑡}}^{(0)}`$ may be given either by Eq. (3) or Eq. (4) but in both cases the stationary tubulin–t concentration, denoted by $`c_t^{(0)}`$, is determined self–consistently as described below. The stationary solutions $`p_{g,s}^{(0)}`$ allow an analytical calculation of the integrals in Eq. (15) and a nonlinear equation in $`c_t^{(0)}`$ follows,

$`c_0c_t^{(0)}`$ $`=`$ $`{\displaystyle \frac{\nu \gamma v_g}{f_{cat}^{(0)}}}\left({\displaystyle \frac{1}{\alpha }}+{\displaystyle \frac{1}{f_{cat}^{(0)}}}(1+\beta )\right).`$ (27)
From this equation the stationary tubulin concentration $`c_t^{(0)}`$ can be determined as a function of the overall concentration of tubulin dimers $`c_0`$ and as function of the other parameters. In Eq. (27) the abbreviation for the velocity ratio

$$\beta =\frac{v_g}{v_s}$$
(28)
has been introduced and the respective length distributions $`p_g^{(0)}`$ and $`p_s^{(0)}`$ follow for a given value of $`c_t^{(0)}`$ via Eqs. (26). The stationary value $`c_t^{(0)}`$ for the reduced model, described by Eqs. (II B 2), follows from Eq. (27) in the limit $`\beta 0`$.

In the range of $`\alpha `$ much larger than the catastrophe rate $`f_{cat}^{(0)}`$ the stationary tubulin–t concentration $`c_t^{(0)}`$ becomes independent of it, because all tubulin–d dimers, that are released during the depolymerization, are immediately regenerated to tubulin–t dimers. Both a large nucleation rate $`\nu `$ and a large growth velocity $`v_g`$ lead to a high consumption of tubulin–t and therefore to a lower stationary concentration $`c_t^{(0)}`$. This tendency is illustrated by the difference between the two curves in Figure 2. On the other hand, large values of the amplitude of the respective catastrophe rate, either $`f`$ or $`\overline{f}`$, act against long microtubules which consist of many tubulin dimers and therefore enhance the densities $`c_t^{(0)}`$ and $`c_d^{(0)}`$. The length distribution of the growing respective shrinking microtubules is determined by $`v_g/f_{cat}^{(0)}`$, which also depends via the catastrophe rate on the concentration $`c_t^{(0)}`$.

These tendencies becomes even more obvious if the linear dependence of the catastrophe rate in Eq. (4) is chosen for the special case $`c_u=c_0`$ and Eq. (27) is expanded in the limit of small and large values of $`\alpha `$. In both cases we obtain the simple formulas

$`c_t^{(0)}`$ $`=`$ $`c_0\left({\displaystyle \frac{\nu \gamma v_g(1+\beta )}{\overline{f}}}\right)^{1/3}(\alpha \overline{f}),`$ (30)
$`c_t^{(0)}`$ $`=`$ $`c_0\left({\displaystyle \frac{\nu \gamma v_g}{\alpha \overline{f}}}\right)^{1/2}(\alpha \overline{f}),`$ (31)
which reflect the described tendencies.

### B Model II

Stationary solutions for model II can be calculated in a similar manner as discussed in the previous section for model I. The length distribution of growing microtubules is again given by Eq. (26) and the integral in Eq. (25) can be calculated analytically. $`c_d`$ may be eliminated from Eq. (25) by using Eq. (11) with $`\dot{c}_t=0`$ and by setting $`\gamma =\eta \lambda `$. Then the nonlinear equation for tubulin–t dimers $`c_t^{(0)}`$ takes the form

$`c_0c_t^{(0)}={\displaystyle \frac{\nu \eta \lambda v_g}{f_{cat}^{(0)}}}\left({\displaystyle \frac{1}{\alpha }}+{\displaystyle \frac{1}{\chi }}+{\displaystyle \frac{1}{f_{cat}^{(0)}}}\right).`$ (32)
Eq. (32) is invariant under permutation $`\alpha \chi `$. If $`\alpha `$ and $`\chi `$ become much larger than the catastrophe rate, the stationary concentration $`c_t^{(0)}`$ becomes rather independent of both. In the limits $`v_s\mathrm{}`$ in Eq. (27) and $`\chi \mathrm{}`$ in Eq. (32) we again obtain the concentration $`c_t^{(0)}`$ for the reduced model. No stationary solution is possible in the limit $`\chi 0`$ because in this limit all tubulin–d dimers are stored in oligomers and the polymerization cycle becomes interrupted.

## IV Threshold for oscillatory polymerization

Stationary microtubule polymerization becomes unstable against oscillating modes in the range of high tubulin dimer concentrations $`c_0`$ and the parameter range where this happens is calculated by a linear stability analysis. Starting from the model equations given in Sec. II we derive linear equations for small perturbations with respect to the stationary state and such perturbations exhibit an exponential time dependence, $`e^{\sigma t}`$. For the exponential factor $`\sigma `$ we derive a nonlinear equation from which both the critical dimer concentration $`c_{0c}`$ and the critical frequency $`\omega _c`$ for the Hopf bifurcation is calculated numerically for various parameter combinations and in limiting cases also analytically.

### A Model I

We introduce small perturbations $`p_{g,s}^{(1)}`$ and $`c_t^{(1)}`$ with respect to the stationary solutions $`p_{g,s}^{(0)}`$ and $`c_t^{(0)}`$ as determined in the last section. With the ansatz

$`p_{g,s}`$ $`=`$ $`p_{g,s}^{(0)}+p_{g,s}^{(1)},`$ (34)
$`c_t`$ $`=`$ $`c_t^{(0)}+c_t^{(1)},`$ (35)
one obtains, after linearization of Eqs. (II B) and Eq. (15), the following set of linear equations describing the dynamics of the perturbations

$`_tp_g^{(1)}`$ $`=`$ $`\left(f_{cat}^{(0)}+v_g_l\right)p_g^{(1)}p_g^{(0)}f_{cat}^{(1)},`$ (37)
$`_tp_s^{(1)}`$ $`=`$ $`v_s_lp_s^{(1)}+f_{cat}^{(0)}p_g^{(1)}+p_g^{(0)}f_{cat}^{(1)},`$ (38)
$`_tc_t^{(1)}`$ $`=`$ $`\alpha c_t^{(1)}`$ (40)
$`\gamma {\displaystyle _0^{\mathrm{}}}𝑑l\left(v_gp_g^{(1)}+\alpha l(p_g^{(1)}+p_s^{(1)})\right).`$
Here $`f_{cat}^{(1)}`$ is the first order contribution of an expansion of the catastrophe rate $`f_{cat}=f_{cat}^{(0)}+f_{cat}^{(1)}+\mathrm{}`$ with respect to the perturbation $`c_t^{(1)}`$:

$$f_{cat}^{(1)}=f_{cat}^{(0)}\frac{c_t^{(1)}}{c_f}.$$
(41)
Since the first order linear equations (IV A) have constant coefficients, their solutions depend exponentially in time and $`c_t^{(1)}`$ may be written as

$$c_t^{(1)}=Ae^{\sigma t}+c.c.$$
(42)
(c.c.=conjugate complex). With this ansatz the three equations in (IV A) can easily be integrated and the solutions for the growing and shrinking microtubules are given by

$`p_g^{(1)}`$ $`=`$ $`{\displaystyle \frac{\nu f_{cat}^{(0)}}{v_gc_f\sigma }}\mathrm{exp}\left(\sigma t{\displaystyle \frac{f_{cat}^{(0)}}{v_g}}l\right)`$ (45)
$`\left[\mathrm{exp}\left({\displaystyle \frac{\sigma }{v_g}}l\right)1\right]A+c.c.,`$
$`p_s^{(1)}`$ $`=`$ $`{\displaystyle \frac{\nu f_{cat}^{(0)}}{v_sc_f\sigma }}\mathrm{exp}(\sigma t{\displaystyle \frac{f_{cat}^{(0)}}{v_g}}l)[k_1\mathrm{exp}({\displaystyle \frac{\sigma }{v_g}}l)`$ (46)
$`+`$ $`k_2+K\mathrm{exp}\left({\displaystyle \frac{\sigma }{v_s}}l\right)]A+c.c..`$ (47)
Herein we have introduced the abbreviations

$`k_1`$ $`=`$ $`{\displaystyle \frac{f_{cat}^{(0)}}{\sigma (1+\beta )+f_{cat}^{(0)}}},`$ (48)
$`k_2`$ $`=`$ $`{\displaystyle \frac{\sigma f_{cat}^{(0)}}{f_{cat}^{(0)}+\sigma \beta }},`$ (49)
and the boundary condition in Eq. (9) requires a vanishing integration constant $`K=0`$. The boundary condition for the time–dependent part of the growing microtubules, $`p_g^{(1)}(l=0,t)=0`$, is also fulfilled. According to the analytic expressions for $`p_g^{(1)}`$ and $`p_s^{(1)}`$ given in Eqs. (42) both may be eliminated in Eq. (40). The remaining integral in Eq. (40) can be calculated analytically and the nonlinear dispersion relation for $`\sigma `$ follows

$`1+\sigma (\sigma +\alpha )G+{\displaystyle \frac{\alpha }{f_{cat}^{(0)}}}\left(1+\beta {\displaystyle \frac{f_{cat}^{(0)}\sigma }{f_{cat}^{(0)}+\sigma \beta }}\right)`$ (50)
$`{\displaystyle \frac{f_{cat}^{(0)}}{f_{cat}^{(0)}+\sigma }}\left[1+{\displaystyle \frac{\alpha }{f_{cat}^{(0)}+\sigma }}\left(1+\beta {\displaystyle \frac{f_{cat}^{(0)}}{f_{cat}^{(0)}+\sigma (1+\beta )}}\right)\right]=0,`$ (51)
with a reduced parameter

$$G=\frac{c_f}{\gamma \nu v_g}$$
(52)
for the catastrophe rate given in Eq. (3) and with

$$G=\frac{f_{cat}^{(0)}}{\overline{f}\gamma \nu v_g}$$
(53)
for the rate given in Eq. (4). After a few rearrangements the dispersion relation in Eq. (50) can be written as a fourth order polynomial in $`\sigma `$

$`\sigma ^4G\beta \left(1+\beta \right)+\sigma ^3G\left[\alpha \beta (1+\beta )+f_{cat}^{(0)}\left(1+3\beta +\beta ^2\right)\right]`$ (54)
$`+\sigma ^2\left[G\alpha f_{cat}^{(0)}\left(1+3\beta +\beta ^2\right)+\left(\beta +2Gf_{cat}^{(0)}{}_{}{}^{2}\right)\left(1+\beta \right)\right]`$ (55)
$`+\sigma [\alpha (1+\beta )(1+\beta +2Gf_{cat}^{(0)^2})+Gf_{cat}^{(0)^3}`$ (56)
$`+f_{cat}^{(0)}(1+2\beta )]+\alpha f_{cat}^{(0)}[2+2\beta +Gf_{cat}^{(0)^2}]+f_{cat}^{(0)^2}=0.`$ (57)
This polynomial describes the linear stability of the stationary solutions given by Eq. (26) and Eq. (27) completely and they are unstable in the parameter range where the growth rate becomes positive, $`Re(\sigma )>0`$. Keeping for instance all parameters besides the dimer concentration $`c_0`$ fixed, then the neutral stability condition $`Re(\sigma )=0`$ provides an equation for the critical dimer concentration $`c_{0c}`$. For concentrations larger than this critical value, $`c_0>c_{0c}`$, the stationary solutions are unstable.

The smallest critical dimer concentrations $`c_{0c}`$ for an oscillatory polymerization are required if the parameters $`\alpha `$,$`\beta `$ and $`G`$ take intermediate values, as discussed in more detail below. At the threshold the real part $`Re(\sigma )=0`$ vanishes and the imaginary part of $`\sigma `$ is the so-called Hopf–frequency $`\omega _c=Im(\sigma )`$. In this special case with a purely imaginary $`\sigma `$ the polynomial in Eq. (54) can be decomposed into its real and imaginary parts giving two coupled equations for the determination of the two unknowns $`f_{cat}^{(0)}`$ and $`\omega _c`$. Having determined $`f_{cat}^{(0)}`$ numerically, $`c_{0c}`$ may be calculated via Eq. (27).

#### 1 Limiting cases with the rate in Eq. (3)

For the limiting cases $`\beta 0`$, $`\beta \mathrm{}`$, $`\alpha 0`$ and $`\alpha \mathrm{}`$ analytical expressions can be given for both the threshold concentration $`c_{0_c}`$ and the Hopf frequency $`\omega _c`$. This is explained at first for the catastrophe rate given by Eq. (3) and for the parameter $`G`$ given in Eq. (52). At threshold one has $`\sigma =i\omega _c`$ and two equations follow from the nonlinear dispersion relation in Eq. (54) which determine the two unknowns $`\omega _c`$ and $`f_{cat}^{(0)}`$. The critical initial concentration $`c_{0c}`$ follows via $`f_{cat}^{(0)}`$ from Eq. (27).

##### a  $`\alpha \mathrm{}`$:

In this limit one obtains from (54)

$`f_{cat}^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{\alpha G}}{\displaystyle \frac{1+\beta }{1+\beta +\beta ^2}},`$ (59)
$`\omega _c`$ $`=`$ $`\sqrt{{\displaystyle \frac{1+\beta }{G\beta }}}.`$ (60)
Accordingly the critical tubulin concentration diverges as $`c_{0c}\alpha ^2`$, that agrees with the full numerical results shown in Fig. 3, besides small logarithmic corrections. In this limit the Hopf frequency $`\omega _c`$ becomes independent of $`\alpha `$ and with increasing values of $`\beta `$ it decreases slightly to a constant value $`\omega _c\sqrt{1/G}`$.

##### b $`\alpha 0`$:

In this case $`\omega _c\sqrt{1/G}`$ becomes also independent of $`\alpha `$ and the catastrophe rate vanishes as $`f_{cat}^{(0)}\alpha `$. Therefore the critical tubulin concentration diverges according to Eq. (27) as $`c_{0c}\alpha ^2`$.

##### c $`\beta 0`$:

In this limit one obtains

$$f_{cat}^{(0)}=\sqrt{\frac{\beta }{G}}\mathrm{and}\omega _c=\left(\frac{\alpha ^2}{G}\right)^{1/4}\left(\frac{1}{\beta }\right)^{1/4}.$$
(61)
The Hopf frequency diverges with increasing values of the shrinking velocity as $`\omega _cv_s^{1/4}`$ in agreement with the numerical results shown in Fig. 4. With this expression for $`f_{cat}^{(0)}`$ one obtains via Eq. (27) for the critical initial concentration $`c_{0c}G/\beta +c_f\mathrm{ln}(f(G/\beta )^{1/2})`$. For medium parameter values this is essentially $`c_{0c}1/\beta `$ as indicated in Fig. 4.

In experiments the shrinking velocity was always larger then the growth velocity, therefore the limit $`\beta \mathrm{}`$ is discarded.

#### 2 Limiting cases for the rate in Eq. (4)

The tendencies for the parameter dependence of the threshold for the Hopf bifurcation as discussed in Sec. IV A 1 are by far not a special property of the choice of the catastrophe rate in Eq. (3). Therefore we consider the same limiting cases as before for the catastrophe rate given in Eq. (4) and with $`G`$ as defined in Eq. (53).

a. $`\alpha \mathrm{}`$: In this limit one obtains from Eq. (54)

$`f_{cat}^{(0)}`$ $`=`$ $`\left({\displaystyle \frac{g(1+\beta )}{\alpha (1+\beta +\beta ^2)}}\right)^{1/2},`$ (63)
$`\omega _c`$ $`=`$ $`{\displaystyle \frac{\left(\alpha g(1+\beta +\beta ^2)(1+\beta )\right)^{1/4}}{\beta ^{1/2}}},`$ (64)
with $`g=\overline{f}\gamma \nu v_g`$. Hence the critical tubulin concentration required for a Hopf bifurcation diverges as $`c_{0c}\alpha `$.

b. $`\alpha 0`$: In this limit one has $`f_{cat}^{(0)}\alpha `$ and $`c_{0c}\alpha ^2`$ diverges too.

c. $`\beta 0`$: For this limit we obtain

$`f_{cat}^{(0)}`$ $`=`$ $`(g\beta )^{1/3},`$ (66)
$`\omega _c`$ $`=`$ $`g^{1/6}\alpha ^{1/2}\left({\displaystyle \frac{1}{\beta }}\right)^{1/3}.`$ (67)
This confirms the importance of a finite ratio of $`\beta =v_g/v_s`$, because the threshold diverges for $`\beta 0`$, similar as for the catastrophe rate given in Eq. (3).

#### 3 Traveling waves solutions

At the threshold of the Hopf bifurcation the rate $`\sigma `$ is purely imaginary, $`\sigma =i\omega _c`$, and the expressions given in Eq. (42) and Eqs. (42) are oscillatory in time

$`c_t^{(1)}`$ $`=`$ $`2A\mathrm{cos}(\omega _ct),`$ (69)
$`p_g^{(1)}`$ $`=`$ $`{\displaystyle \frac{S_1}{v_g}}\mathrm{exp}({\displaystyle \frac{f_{cat}^{(0)}}{v_g}}l)\left[\mathrm{sin}(\omega _ct)\mathrm{sin}(\omega _c(tl/v_g))\right],`$ (70)
$`p_s^{(1)}`$ $`=`$ $`{\displaystyle \frac{S_1}{v_s}}\mathrm{exp}({\displaystyle \frac{f_{cat}^{(0)}}{v_g}}l)[k_2\mathrm{sin}(\omega _ct+\phi _2)`$ (71)
$`+`$ $`k_1\mathrm{sin}(\omega _c(tl/v_g)+\phi _1)],`$ (72)
whereby the following abbreviations for the amplitudes

$`S_1`$ $`=`$ $`{\displaystyle \frac{2A\nu f_{cat}^{(0)}}{\omega _cc_f}},`$ (74)
$`k_1`$ $`=`$ $`{\displaystyle \frac{\sqrt{f_{cat}^{(0)}{}_{}{}^{4}+\omega _c^2f_{cat}^{(0)}{}_{}{}^{2}(1+\beta )^2}}{f_{cat}^{(0)}{}_{}{}^{2}+\omega _c^2(1+\beta )^2}},`$ (75)
$`k_2`$ $`=`$ $`{\displaystyle \frac{\sqrt{\left(\omega _c^2\beta f_{cat}^{(0)}{}_{}{}^{2}\right)^2+f_{cat}^{(0)}{}_{}{}^{2}\omega _c^2\left(1+\beta \right)^2}}{f_{cat}^{(0)}{}_{}{}^{2}+\left(\omega _c\beta \right)^2}},`$ (76)
and phases

$`\phi _1`$ $`=`$ $`\mathrm{arctan}\left({\displaystyle \frac{\omega _c(1+\beta )}{f_{cat}^{(0)}}}\right),`$ (78)
$`\phi _2`$ $`=`$ $`\mathrm{arctan}\left({\displaystyle \frac{\omega _cf_{cat}^{(0)}(1+\beta )}{\beta \omega _c^2f_{cat}^{(0)^2}}}\right),`$ (79)
have been introduced. The analytical expressions for $`p_g^{(1)}`$ and $`p_s^{(1)}`$ indicate that the time–dependent contribution to the length distribution of the microtubules includes homogeneous amplitude oscillations and waves with a wavelength $`\omega _c/v_g`$ that travel to larger values of the length $`l`$. Hence, the length distributions $`p_{g,s}^{(1)}`$ depend on two different length scales, the decay length $`v_g/f_{cat}^{(0)}`$ and the wave length $`v_g/\omega _c`$ of the traveling waves. With the explicit solutions for $`c_t^{(1)}`$ and $`p_g^{(1)}`$ the phase of the oscillating part of the tubulin–d concentration relative to the phase of $`c_t^{(1)}`$ as well as its oscillation amplitude is calculated via Eq. (11).

#### 4 Numerical results for the threshold of model I

Since $`\sigma `$ is purely imaginary at the threshold of a Hopf bifurcation, $`\sigma =i\omega _c`$, the dispersion relation in Eq. (54) can be decomposed in its real and imaginary part and from these two equations the critical concentration $`c_{0c}`$ and the frequency $`\omega _c`$ may be calculated numerically as a function of the parameters. The numerical calculations in this section are restricted to the exponential tubulin dependence of the catastrophe rate as given by Eq. (3).

The critical tubulin dimer concentration $`c_{0c}`$ and the critical frequency $`\omega _c`$ at the Hopf bifurcation are shown in Fig. 3 as function of the regeneration rate $`\alpha `$ and in Fig. 4 as function of the velocity ratio $`1/\beta =v_s/v_g`$, whereby the reduced parameter $`G`$ has been chosen at the values $`G=3000`$ and $`G=300`$, respectively. Since $`G`$ includes a number of parameters the curves in both figures represent a larger parameter set. In the limit of a vanishing regeneration and in the limit of very large values of $`\alpha `$, where the regeneration process is much faster than any other process, the critical tubulin concentration $`c_{0c}`$ diverges and therefore the Hopf bifurcation is suppressed. In addition, both figures indicate that the smallest values of the critical tubulin concentration $`c_{0c}`$ are obtained at intermediate values of the parameters $`\alpha `$, $`\beta `$ and $`G`$. The location of the threshold minima, however, depends on the actual values of the rest of the parameters. The frequency $`\omega _c`$ becomes rather small in the limit $`\alpha 0`$ and for large values of $`\alpha `$ this frequency becomes independent of it, cf. section IV A 1.

In Fig. 4 the threshold minimum is less pronounced than in Fig. 3 and in the limit $`\beta =v_g/v_s0`$ the threshold $`c_{0c}`$ increases linearly with $`v_s`$ and in agreement with limits given in section IV A 1. Accordingly, there is no Hopf bifurcation for the reduced model that follows in the limit $`\beta 0`$ as described in section II B 2. Hence, the dynamics of shrinking microtubules is one essential degree of freedom favoring oscillating microtubule polymerization. The dynamics of oligomers, as discussed in Sec. IV B, is an alternative degree of freedom that favors oscillations.

For large values of $`v_s`$ the frequency $`\omega _c`$ becomes large too and the oscillation period becomes much shorter than any relaxational dynamics of $`p_g`$ and $`c_t`$. According to the quick shrinking, the life time of a depolymerizing microtubules vanishes and therefore the amplitude of the density of shrinking microtubules is small too, $`p_s1/v_s`$. In other words, in the limit of large values of $`v_s`$ the intermediate step of shrinking microtubules may be neglected and the transition from $`p_g`$ to tubulin-d dimers is effectively a direct process as explicitly assumed for the reduced model. If either the regeneration or the shrinking dynamics becomes too fast, the Hopf bifurcation is suppressed. The two intermediate steps, the depolymerization and the regeneration, act obviously as antagonistic steps or jam processes that favors oscillations.

Since one has at threshold $`\sigma =i\omega _c`$, $`p_g^{(1)}`$ and $`p_s^{(1)}`$ in Eq. (70) and Eq. (71) include both traveling wave contributions $`\mathrm{exp}[i(\omega _ctkl)]`$ with a wave number $`k=\omega _c/v_g`$ and that always travel towards larger lengths of the microtubules. The length distribution is exponentially decaying on the length scale $`v_g/f_{cat}^{(0)}`$. If this is large compared to the wave length $`\lambda =2\pi v_g/\omega _c`$, as for instance in the limit $`v_sv_g`$, then one has a kind of self averaging in the respective integrals and the Hopf bifurcation is suppressed.

The phase difference between the oscillations of tubulin–t and the oscillations of the total amount of polymerized tubulin, described by $`L(t)`$ in Eq. (14), is another experimentally accessible quantity . The difference between the phases of the oscillatory contributions of $`c_d`$ and $`c_t`$ as well as the difference between the phases of $`L(t)`$ and $`c_t`$ are given in Fig. 5. These phase differences as well as the ratios between the amplitudes of the fields, cf. lower part of Fig. 5, are calculated at the threshold of the Hopf bifurcation by using the analytical solutions calculated in Sec. IV A 3.

For large values of $`\alpha `$, tubulin–d is quickly regenerated into tubulin–t and therefore the density $`c_d`$ becomes smaller as shown by the lower part in Fig. 5. In the opposite limit of small values of $`\alpha `$ tubulin–t is consumed by nucleation and growth of microtubules, but the source, which is supplied by the regeneration of tubulin–d, decays and therefore one obtains large values for the ratio between the amplitude of $`c_d^{(1)}`$ and $`c_t^{(1)}`$ as well as between the amplitudes of $`L^{(1)}`$ and $`c_t^{(1)}`$.

The decay of the ratio between the amplitudes of $`L^{(1)}`$ and $`c_t^{(1)}`$ is less obvious. $`L^{(1)}=L_g^{(1)}+L_s^{(1)}`$ has the two contributions $`L_g^{(1)}=\gamma _0^{\mathrm{}}𝑑llp_g^{(1)}=\overline{A}\mathrm{cos}(\omega _ct+\phi _{\overline{A}})`$ and $`L_s^{(1)}=\gamma _0^{\mathrm{}}𝑑llp_s^{(1)}=\overline{B}\mathrm{cos}(\omega _ct+\phi _{\overline{B}})`$. The two amplitudes $`\overline{A}`$ and $`\overline{B}`$ increase with the regeneration rate $`\alpha `$. However, with increasing values of $`\alpha `$ the phase difference $`\phi _{\overline{A}}\phi _{\overline{B}}`$ increases as well up to unity leading to an effective decay of the sum $`L^{(1)}`$ as shown by the lower part of Fig. 5.

The phase shifts of $`c_d`$ and $`L(t)`$ with respect to the phase of $`c_t`$ are rather independent of the regeneration rate $`\alpha `$ as shown in Fig. 5. The absolute values of these shifts are in qualitative agreement with the expectation as described in the following. At the maximum of $`c_t`$ the catastrophe rate takes its minimum and therefore, since the nucleation and the growth velocity are constant, $`L(t)`$ is increasing for a while and up to the moment when enough $`c_t`$ is consumed and the catastrophe rate increases again. Due to an increasing decay of microtubules the maximum of the latter will also lead to a delayed maximum for $`c_d`$. For large values of $`\alpha `$ the amount of polymerized tubulin is nearly in anti phase with respect to $`c_t`$, which is also mentioned in Ref. . The slightly stronger $`\alpha `$–dependence of the phase difference between $`L(t)`$ and $`c_t`$ is mainly due to the $`\alpha `$-dependence of shrinking microtubules $`p_s`$, because the relative phase of $`p_g^{(1)}`$ is nearly independent of $`\alpha `$ (see also Section IV B).

### B Model II

Here the stability of the stationary polymerization state of model II, described by $`c_t^{(0)}`$,$`c_d^{(0)}`$ and $`p_g^{(0)}`$, is investigated with respect to small perturbations $`c_t^{(1)}`$, $`c_d^{(1)}`$ and $`p_g^{(1)}`$. With the ansatz

$`p_g`$ $`=`$ $`p_g^{(0)}+p_g^{(1)},`$ (81)
$`c_{t,d}`$ $`=`$ $`c_{t,d}^{(0)}+c_{t,d}^{(1)},`$ (82)
the equations for model II are linearized with respect to these perturbations and one obtains the following set of linear differential equations with constant coefficients

$`_tp_g^{(1)}`$ $`=`$ $`f_{cat}^{(1)}p_g^{(0)}\left(f_{cat}^{(0)}+v_g_l\right)p_g^{(1)},`$ (84)
$`_tc_t^{(1)}`$ $`=`$ $`\eta \lambda v_g{\displaystyle _0^{\mathrm{}}}𝑑lp_g^{(1)}+\alpha c_d^{(1)},`$ (85)
$`_tc_d^{(1)}`$ $`=`$ $`\chi \left(c_t^{(1)}+c_d^{(1)}+\eta \lambda {\displaystyle _0^{\mathrm{}}}𝑑llp_g^{(1)}\right)\alpha c_d^{(1)}.`$ (86)
$`f_{cat}^{(1)}`$ is the first order correction with respect to its value in the stationary state and it is given by equation (41).

The time–dependent contributions to the tubulin–t and tubulin–d dimer densities are described by

$`c_t^{(1)}`$ $`=`$ $`Ae^{\sigma t}+c.c.,`$ (88)
$`c_d^{(1)}`$ $`=`$ $`EAe^{\sigma t}+c.c.,`$ (89)
with the common complex amplitude $`A`$ and relative complex factor $`E`$ that describes via $`E=|E|e^{i\phi _d}`$ the amplitude ratio $`|E|`$ as well as the phase difference $`\phi _d`$ between both fields. With the solution for growing microtubules as given in Eq. (45) we can eliminate $`c_t^{(1)}`$ and $`p_g^{(1)}`$ in Eq. (85) and we again obtain from the resulting solubility condition a nonlinear dispersion relation for the exponential factor $`\sigma `$

$`\left(G\sigma +{\displaystyle \frac{1}{\sigma +f_{cat}^{(0)}}}\right)(\sigma +\alpha +\chi )`$ (90)
$`+`$ $`\alpha \chi \left(G+{\displaystyle \frac{\sigma +2f_{cat}^{(0)}}{f_{cat}^{(0)}\left(\sigma +f_{cat}^{(0)}\right)^2}}\right)=0,`$ (91)
with $`G=c_f/(\eta \lambda v_g\nu )`$. After a few rearrangements of this equation one obtains a fourth order polynomial in $`\sigma `$ for model II as well

$`\sigma ^4f_{cat}^{(0)}G+\sigma ^3f_{cat}^{(0)}G\left(2f_{cat}^{(0)}+\alpha +\chi \right)`$ (92)
$`+`$ $`\sigma ^2f_{cat}^{(0)}\left[1+\alpha \chi G+Gf_{cat}^{(0)}\left(2\chi +f_{cat}^{(0)}+2\alpha \right)\right]`$ (93)
$`+`$ $`\sigma [f_{cat}^{(0)}(f_{cat}^{(0)}+\alpha +\chi )+\alpha \chi `$ (94)
$`+`$ $`f_{cat}^{(0)}{}_{}{}^{2}G(2\alpha \chi +f_{cat}^{(0)}(\alpha +\chi ))]`$ (95)
$`+`$ $`f_{cat}^{(0)}\alpha \chi \left(Gf_{cat}^{(0)}{}_{}{}^{2}+2\right)+f_{cat}^{(0)}{}_{}{}^{2}(\alpha +\chi )=0,`$ (96)
that determines the linear stability of the stationary polymerization for model II. Again we are interested in the neutrally stable case, $`Re(\sigma )=0`$, that separates the stable from the unstable regime. At the neutral stability point of the Hopf bifurcation one has $`\omega _c=Im(\sigma )`$ and Eq. (92) can be decomposed into its real and imaginary part. From these two equations $`f_{cat}^{(0)}`$ and $`\omega _c`$ are determined by standard methods. $`c_{0c}`$ may be calculated via Eq. (32).

#### 1 Traveling waves solutions

At the Hopf bifurcation the non-stationary part of growing microtubules is again described by the distribution given by Eq. (70) and the fields $`c_{oli}`$ and $`c_d`$ are not in phase with $`c_t`$ in general. The two fields may be written in terms of the amplitude ratio $`|E|`$ and the relative phase $`\phi _d`$ in the following form

$`c_t^{(1)}`$ $`=`$ $`2A\mathrm{cos}(\omega _ct),`$ (98)
$`c_d^{(1)}`$ $`=`$ $`2A|E|\mathrm{cos}(\omega _ct+\phi _d).`$ (99)
The amplitude ratio $`|E|`$ and the phase $`\phi _d`$ can be determined from the two coupled equations (85) and (86) and they are given by

$`|E|`$ $`=`$ $`{\displaystyle \frac{\sqrt{f_{cat}^{(0)}{}_{}{}^{2}+\omega _c^2(G(f_{cat}^{(0)}{}_{}{}^{2}+\omega _c^2)1)^2}}{\alpha G(f_{cat}^{(0)}{}_{}{}^{2}+\omega _c^2)}},`$ (100)
$`\phi _d`$ $`=`$ $`\mathrm{arctan}\left({\displaystyle \frac{\omega _c}{f_{cat}^{(0)}}}\left(G(f_{cat}^{(0)}{}_{}{}^{2}+\omega _c^2)1\right)\right).`$ (101)
In a similar manner the oligomer density $`c_{oli}^{(1)}`$ may be written in terms of an amplitude ratio $`|F|`$ and a relative phase $`\phi _{oli}`$ between $`c_{oli}^{(1)}`$ and $`c_t^{(1)}`$

$`c_{oli}^{(1)}`$ $`=`$ $`2A|F|\mathrm{cos}(\omega _ct+\phi _{oli}),`$ (102)
with

$`|F|`$ $`=`$ $`{\displaystyle \frac{\sqrt{\alpha ^2+\omega _c^2}}{\chi }}|E|,\mathrm{and}`$ (103)
$`\phi _{oli}`$ $`=`$ $`\mathrm{arctan}\left({\displaystyle \frac{\alpha \mathrm{tan}(\phi _d)+\omega _c}{\alpha \omega _c\mathrm{tan}(\phi _d)}}\right).`$ (104)
The oscillatory contribution to the polymerized tubulin $`L^{(1)}`$ can also be written as a harmonic function, $`L^{(1)}=2A|H|\mathrm{cos}(\omega _ct+\phi _L)`$. For both, the amplitude ratio $`|H|`$ and the relative phase $`\phi _L`$, one obtains long expressions that are not presented here. The phase shifts and the amplitude ratios between the oscillating fields are shown in Figure 7 as function of the regeneration rate $`\alpha `$. As discussed in Section IV A, the phase shift of the polymerized tubulin $`L^{(1)}(t)`$ with respect to $`c_t^{(1)}`$ is rather independent of $`\alpha `$, whereas the phase of oligomer oscillations change slightly with $`\alpha `$. A phase shift $`\pi `$ between the polymerized tubulin and oligomers is measured in experiments, cf. Refs.  and . In this model this is only possible in the limit of a dissociation rate $`\chi `$ much smaller than the regeneration rate $`\alpha `$.

#### 2 Numerical results for the threshold of model II

At the threshold one has again $`\sigma =i\omega _c`$ and from the imaginary together with the real part of the dispersion relation in Eq. (92) the critical concentration $`c_{0c}`$ and the Hopf frequency $`\omega _c`$ may be calculated as function of the parameters. Also for model II we restrict our numerical calculations to the catastrophe rate with the exponential dependence given in Eq. (3).

The critical tubulin concentration $`c_{0c}`$ and the critical frequency $`\omega _c`$ at the Hopf bifurcation are shown in Fig. 6 as function of the regeneration rate $`\alpha `$ and for two different values of the decay rate of oligomers $`\chi `$, whereby for the reduced parameter $`G`$ the value $`G=3000`$ has been chosen. Since $`G`$ includes a number of parameters the curves in both parts represent a larger parameter set. For a fixed finite value for $`\chi `$ in the limit of a vanishing regeneration $`\alpha 0`$, where the polymerization cycle is interrupted, and in the limit of very large values of $`\alpha `$, where the regeneration process is much faster than any other process, the critical tubulin concentration $`c_{0c}`$ diverges similar as for model I and therefore the Hopf bifurcation is suppressed. If $`\alpha `$ is kept fixed at a medium value the threshold curve $`c_{0c}(\chi )`$ as function of the decay rate $`\chi `$ for oligomers has a similar shape as shown as function of $`\alpha `$ in Fig. 6. The critical tubulin concentration $`c_{0c}`$ also takes its smallest values at intermediate values of $`\alpha `$, $`\chi `$ and $`G`$, whereby the location of the threshold minima depends on the actual values of the rest of parameters. With a decreasing rate $`\alpha 0`$ of tubulin regeneration also the frequency $`\omega _c`$ becomes small. On the other hand for large values of $`\alpha `$ the tubulin regeneration is not anymore a rate limiting factor and the critical frequency $`\omega _c`$ becomes rather independent of $`\alpha `$ , cf. section IV A 1.

The stability of oligomers and therefore the decay rate $`\chi `$ depend very much on the available GTP: Increasing GTP concentrations destabilize oligomers and increase the decay rate $`\chi `$ . With increasing GTP concentrations also the rate $`\alpha `$ of the transition from $`c_d`$ to $`c_t`$ is enhanced. However, if the tubulin regeneration and the oligomer decay become to quick, an oscillatory polymerization is suppressed. In other words, if one increases $`\alpha `$ and $`\chi `$ beyond some minimum values, the threshold concentration for tubulin increases too. Such a tendency for the GTP dependence of the oscillation is in agreement with the results reported from experiments .

With increasing values of $`\alpha `$ tubulin–d is again quickly transfered by the regeneration process into tubulin–t, leading to a small amplitude ratio $`c_d^{(1)}/c_t^{(1)}`$. Accordingly more tubulin is left to be stored in $`L^{(1)}`$ and $`c_{oli}^{(1)}`$. Therefore both increase with larger values of $`\alpha `$ as indicated in Fig. 7. This has to be compared with $`L^{(1)}`$ for model I, where it decays as function of $`\alpha `$ because the phase shift between the contributions of the growing and shrinking microtubules changes too. For model II the relative phases are also rather independent of $`\alpha `$, whereby due to the quick regeneration of $`c_d`$ the relative phase between $`c_t^{(1)}`$ and $`c_d^{(1)}`$ is slightly decreasing.

#### 3 Reduced models

The dispersion relation for model I and II, considered in the previous section, became equivalent in the limits $`\beta 0`$ and $`\chi \mathrm{}`$ and in both cases one obtains the same dispersion relation

$`\sigma ^3Gf_{cat}^{(0)}+\sigma ^2\left[G\alpha f_{cat}^{(0)}+2Gf_{cat}^{(0)}{}_{}{}^{2}\right]`$ (105)
$`+`$ $`\sigma \left[\alpha \left(1+2Gf_{cat}^{(0)^2}\right)+Gf_{cat}^{(0)^3}+f_{cat}^{(0)}\right]`$ (106)
$`+`$ $`\alpha f_{cat}^{(0)}\left[2+Gf_{cat}^{(0)^2}\right]+f_{cat}^{(0)^2}=0,`$ (107)
with the reduced parameter $`G`$ as given in Eq. (52) for model I and with $`G=c_f/(\eta \lambda v_g\nu )`$ for model II. This polynomial in $`\sigma `$ has always negative growth rates, $`Re(\sigma )<0`$, and therefore stationary solutions are always stable.

## V Numerical method for model I

The two differential equations for growing and shrinking microtubules in (II B) are of first order with respect to the length $`l`$ of the microtubules and first order in time. A straight forward spatial discretization of such first order equations often leads to numerical instabilities. Especially the equation for shrinking microtubules, cf. Eq. (8), has problematic stability properties. For this reason we approximate the solutions of Eqs. (II B) by a two mode ansatz

$`p_g(l,t)`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{f_{cat}^{(0)}}{v_g}}l\right)\left({\displaystyle \frac{\nu }{v_g}}+F_g(l,t)\right),`$ (109)
$`p_s(l,t)`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{f_{cat}^{(0)}}{v_g}}l\right)\left({\displaystyle \frac{\nu }{v_s}}+F_s(l,t)\right),`$ (110)
where the first mode describes just the stationary solution and the second one the oscillatory contribution. This approximation becomes exact close to the threshold and this ansatz leads with Eqs. (II B) to two differential equations for $`F_{g,s}`$

$`_tF_g`$ $`=`$ $`\left(f_{cat}^{(0)}f_{cat}\right)\left({\displaystyle \frac{\nu }{v_g}}+F_g\right)v_g_lF_g,`$ (112)
$`_tF_s`$ $`=`$ $`f_{cat}\left({\displaystyle \frac{\nu }{v_g}}+F_g\right){\displaystyle \frac{v_s}{v_g}}f_{cat}^{(0)}\left({\displaystyle \frac{\nu }{v_s}}+F_s\right)`$ (114)
$`+v_s_lF_s.`$
Both fields may be expanded with respect to the first two spatial Fourier modes $`e^{inlk}`$ ($`n=0,1`$)

$`F_g(l,t)`$ $`=`$ $`B(t)+{\displaystyle \frac{1}{2}}\left(C(t)e^{ikl}+C^{}(t)e^{ikl}\right),`$ (116)
$`F_s(l,t)`$ $`=`$ $`D(t)+{\displaystyle \frac{1}{2}}\left(H(t)e^{ikl}+H^{}(t)e^{ikl}\right),`$ (117)
in order to remove the spatial dependence from Eqs. (V). Herein the wave number is chosen at its value at the threshold of the Hopf bifurcation, $`k=\omega _c/v_g`$. This ansatz leads to a set of ordinary differential equations for the time dependent amplitudes $`B(t),C(t),D(t),H(t)`$ that are described in the following.

Due to Eq. (5) one has the boundary condition $`F_g(l=0,t)=0`$ that gives the relation $`C_R=B`$, with $`Re(C)=C_R`$, between these two functions. Ansatz (116) together with equation (112) leads to the relation

$$Im(C)=C_I=\frac{\nu }{kv_g^2}\left(f_{cat}f_{cat}^{(0)}\right),$$
(118)
and to the first order differential equation in $`B`$

$$_tB=\left(f_{cat}^{(0)}f_{cat}\right)\left(\frac{\nu }{v_g}+B\right).$$
(119)
Ansatz (117) in equation (114) gives the set of coupled differential equations

$`_tD`$ $`=`$ $`f_{cat}\left({\displaystyle \frac{\nu }{v_g}}+B\right){\displaystyle \frac{\nu f_{cat}^{(0)}}{v_g}}{\displaystyle \frac{v_s}{v_g}}f_{cat}^{(0)}D,`$ (121)
$`_tH_R`$ $`=`$ $`f_{cat}B{\displaystyle \frac{v_s}{v_g}}f_{cat}^{(0)}H_Rkv_sH_I,`$ (122)
$`_tH_I`$ $`=`$ $`f_{cat}C_I{\displaystyle \frac{v_s}{v_g}}f_{cat}^{(0)}H_I+kv_sH_R.`$ (123)
With the periodic $`l`$–dependence given in Eqs. (V) the integrals in Eq. (40) can be evaluated and one obtains the following differential equation for the tubulin–t dimer density

$`_tc_t=\gamma v_gK(t)\alpha \gamma L(t)+\alpha c_{0c}\left(1+\epsilon \right)\alpha c_t,`$ (124)
where the coefficients are given by

$`K(t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑lp_g(l,t)={\displaystyle \frac{\nu }{v_g\delta }}+{\displaystyle \frac{k(kB\delta C_I)}{\delta (\delta ^2+k^2)}},`$ (125)
$`L(t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑ll\left[p_g(l,t)+p_s(l,t)\right]`$ (126)
$`=`$ $`{\displaystyle \frac{\nu }{\delta ^2}}({\displaystyle \frac{1}{v_g}}+{\displaystyle \frac{1}{v_s}})+\mathrm{\Delta }(3\delta ^2k^2B+k^4B2\delta ^3kC_I`$ (127)
$`+`$ $`\delta ^2(\delta ^2k^2)H_R2\delta ^3kH_I+(\delta ^2+k^2)^2D),`$ (128)
and where the abbreviations $`\mathrm{\Delta }=[\delta ^2(\delta ^2+k^2)^2]^1`$ and $`\delta =f_{cat}^{(0)}/v_g`$ have been introduced. The reduced control parameter $`\epsilon =(c_0c_{0c})/c_{0c}`$ measures the difference between the tubulin dimer concentration $`c_0`$ and the critical one, $`c_{0c}`$. For $`\epsilon >0`$ sustained oscillations occur but they are damped below threshold $`\epsilon <0`$. For the numerical solution of model II we use either the same approximation scheme, where only the factors of the scheme take a different form, or in the absence of $`p_s`$ a direct spatial discretization provides also a stable algorithm. The five differential equations for model I in Eq. (119), Eq. (119) and Eq. (124) and the corresponding three differential equations for model II are integrated numerically by a second-order Runge-Kutta method with a time step $`\mathrm{\Delta }t=0.01`$.

The time–dependence of the fields of model I are shown in Figure 8 for one parameter set and these fields obviously have different phases. The extrema (maxima) of the polymerized tubulin $`L(t)`$ and the tubulin–d concentration $`c_d(t)`$ are delayed with respect to the extrema of the tubulin–t concentration $`c_t(t)`$, a behavior that is already indicated by the reaction cycle shown in Fig. 1. At the threshold the parameter dependence of these phase differences may be calculated from the formulas given in Sec. IV A 3 and Sec. IV B 1. In Fig. 5 and in Fig. 7 these phases are shown as function of the regeneration rate $`\alpha `$. In Fig. 8 the initial concentration of tubulin $`c_0`$ was chosen very close to the threshold $`c_{0c}`$ with $`\epsilon =0.01`$. At this value of the control parameter the oscillations behave harmonically and the agreement between the numerical solution and the amplitude approximation is rather good, as described in the following section. For larger values of the reduced control parameter the oscillations become anharmonic.

As already been mentioned in the introduction, the length distribution of filaments is a crucial difference between the biochemical reaction discussed in this work and the common oscillatory chemical reactions. For model I we show in Fig. 8d) at two different times and at $`\epsilon =0.01`$ the superposition of the length distribution of growing and shrinking microtubules, cf. $`P(l)=p_g(l)+p_s(l)`$. The exponential decay of the envelope of the length distribution is described by Eq. (70) and Eq. (71), with a decay rate $`v_g/f_{cat}^{(0)}`$, and the modulation is due to the traveling waves in the time–dependent contribution. The amplitude of the shrinking microtubules is rather small for $`\beta =0.1`$, cf. Eq. (71), and therefore the contribution to $`P(l)`$ comes mainly from growing microtubules.

For model II a discretization of the length coordinate in the equation for growing microtubule, cf. Eq. (2), also provides a stable numerical algorithm. Hence, the nonlinear oscillatory solution can be obtained numerically without the approximations as described for model I in Eqs. (V) above. The respective results are shown in Fig. 9a)-c), where the densities $`L,c_t`$ and $`c_{oli}`$ are shown as function of time for three different regeneration rates $`\alpha `$. For both values of $`\alpha `$ in part a) and b) the tubulin concentration is smaller than the corresponding threshold value but the absolute distance $`c_0c_{0c}`$ to the threshold is equal. These transient subthreshold–oscillations are remarkable, because the transient oscillations in experiments might be subthreshold ones in contrast to the common interpretation that the oscillations are transient due to the tubulin–t consumption during the microtubule polymerization.

For the simulations shown in Fig. 9 a narrow length distribution $`p_g`$ has been used as initial condition. In such cases the tubulin-t concentration corresponds almost to the total initial concentration $`c_0`$. Starting with such an initial condition, at first tubulin–t dimers are consumed during the growth of microtubules. This leads to a first maximum of the polymerized tubulin $`L`$, but the oligomers, the decay product of the microtubules, are negligible and as consequence the densities of tubulin-d and tubulin-t drop down too. But a small $`c_t`$ increases the catastrophe rate $`f_{cat}`$ and microtubules decay with a higher rate, which increases the density of oligomers etc.. After a few of such oscillations the densities reach their stationary values for subthreshold concentrations. In the case of a large regeneration rate the oscillation frequency is much higher than for small regeneration rates and more oscillations are performed until the stationary values are reached. The stationary value of the polymerized tubulin and of the oligomers is also much larger in part b) than in part a) whereas the stationary value for the tubulin–t dimers remains nearly constant. The reason is the low density $`c_d`$ in the case of large values of $`\alpha `$.

Far beyond the threshold of the Hopf bifurcation, the oscillations become anharmonic as shown in Fig. 9c). Hereby the oscillations of the total polymerization and of their decay product, $`c_{oli}`$, are nearly in antiphase, similar as in experiments described in Ref.. At the threshold a phase difference of $`\pi `$ was only possible in the limit of large regeneration rates $`\alpha `$ and for a much smaller dissociation rate $`\chi `$. If we consider initial concentrations which are much larger than the corresponding critical concentration, the phase shift of $`\pi `$ between $`L(t)`$ and $`c_{oli}(t)`$ is also possible at intermediate values of $`\chi `$ and $`\alpha `$.

In Fig. 9d) the distribution for growing microtubules $`p_g(l)`$ is shown for model II at three different times. The reduced control parameter $`\epsilon =0.21`$ is rather large compared to its value in Fig. 8 and the curves indicate that the length distribution of microtubules isn’t described anymore by harmonic traveling waves as in the vicinity of the bifurcation. As long as the tubulin-t density is large and the catastrophe rate small, microtubules grow with a constant velocity $`v_g`$ and only a few of them experience a catastrophe. During this period a plateau in the length distribution is build. But after a large amount of tubulin-t has been used up the catastrophe rate $`f_{cat}`$ increases very steeply leading to a strong decay of the microtubules at all lengths. If $`f_{cat}`$ drops down again the low density of long microtubules grow further with a constant velocity $`v_g`$ and short microtubules are nucleated at a higher density. This leads to a step in the distribution that travels with the growth velocity $`v_g`$ to larger values of $`l`$. Therefore, far beyond the oscillatory threshold the temporally anharmonic behavior of $`c_t`$ leads in this manner to step–like length distribution of the microtubules.

## VI Amplitude expansion

Here we focus on a semi-analytical treatment of the oscillating polymerization slightly beyond its onset. The interesting question here is whether the bifurcation to these oscillations is continuous (supercritical) or discontinuous (subcritical). In physical systems Hopf bifurcations are mostly subcritical , but for the models discussed in this work we always find a supercritical one. This is advantageous, because for a supercritical Hopf bifurcation a semi–analytical treatment is possible. The appropriate frame work is the universal amplitude equation of oscillatory fields, cf. Eq. (1). This equation will be derived in the present section from the basic reaction equations of microtubule polymerization.

The perturbation analysis employed for the derivation of Eq. (1) is an expansion of the solutions of the basic equations with respect to small amplitudes of the oscillatory contributions . As a small perturbation parameter the relative difference between the actual tubulin concentration $`c_0`$ and the critical tubulin concentration $`c_{0c}`$ is introduced

$$\epsilon =\frac{c_0c_{0c}}{c_{0c}}.$$
(129)
A signature for the $`\pm `$ symmetry of the oscillatory behavior  is the power law for the oscillation amplitude $`A\sqrt{\epsilon }`$. Accordingly the solutions of the basic equations at the threshold are expanded with respect to powers of $`\sqrt{\epsilon }`$

$`𝐮=𝐮^{(0)}+\epsilon ^{1/2}𝐮^{(1)}+\epsilon 𝐮^{(2)}+\epsilon ^{3/2}𝐮^{(3)}+𝒪(\epsilon ^2),`$ (130)
where the vector notation $`𝐮^{(j)}=(\stackrel{~}{p}_g^{(j)},\stackrel{~}{p}_s^{(j)},\stackrel{~}{c}_t^{(j)})`$ is used with $`j=0,1,2,3`$. The components of $`𝐮^{(j)}`$ differ by a factor of $`\sqrt{\epsilon }`$, such as $`\sqrt{\epsilon }\stackrel{~}{c}_t^{(1)}=c_t^{(1)}`$ etc. . The components of $`𝐮^{(0)}`$ describe the stationary microtubule polymerization as given in Sec. III and the components of $`𝐮^{(1)}`$ describe the linear oscillatory solutions that may be written at the threshold in the following form

$`𝐮^{(1)}=B𝐮_ee^{i\omega _ct}+c.c..`$ (131)
Here $`𝐮_𝐞`$ includes the amplitude ratios between the fields $`c_t^{(1)}`$, $`p_g^{(1)}`$, $`p_s^{(1)}`$ at the threshold and $`\sqrt{\epsilon }B=A`$ as explained below.

Close to the threshold one has $`Re(\sigma )\epsilon 1`$ and the linear solution $`𝐮^{(1)}e^{\sigma t}`$ grows or decays only by a very small amount during one oscillation period $`2\pi /\omega _c`$. These two disparate time scales near the threshold, that of the oscillation period ($`2\pi /\omega _c`$) and that of growth and decay ($`1/\epsilon `$), may be separated within a perturbation expansion by introducing a slow time scale $`T=\epsilon t`$ . The fast time scale is included in the exponential function $`e^{i\omega _ct}`$ and the other one will be described by a slowly varying amplitude $`B(T)`$. Accordingly the linear solution near threshold may be written as

$`𝐮^{(1)}(t,T)=B(T)𝐮_ee^{i\omega _ct}+c.c..`$ (132)
In order to differentiate this product of time dependent functions, instead of applying the chain rule of differentiation, one may replace this operation by the following sum $`_t_t+\epsilon _T`$. Here $`_t`$ acts only on the fast time–dependence occurring in the exponential function and $`_T`$ acts only on the amplitude $`B(T)`$.

Using this replacement and the $`\epsilon `$-expansion of $`𝐮`$ the basic equations given in Sec. II can be ordered with respect to powers of $`\sqrt{\epsilon }`$. In this way we obtain a hierarchy of partial differential equations, which we need up to $`𝒪(\epsilon ^{3/2})`$. The whole procedure is described in more detail in Appendix A. The amplitude equation follows from a solubility condition for the equation at order $`𝒪(\epsilon ^{3/2})`$ and it has the following form

$`\tau _0_TB=\left(1+ia\right)Bg\left(1+ic\right)|B|^2B.`$ (133)
$`\tau _0`$ is the relaxation time, $`a`$ is the linear and $`c`$ is the nonlinear frequency shift. The nonlinear coefficient $`g`$ determines the bifurcation structure. For $`g>0`$ the bifurcation is supercritical (steady) and for $`g<0`$ the bifurcation is subcritical (unsteady). For the coefficients $`\tau _0`$, $`a`$, $`g`$ and $`c`$ one obtains long expressions in terms of the reaction constants of the basic equations, that have been calculated by using computer algebra. The respective formulas are not presented but the parameter dependence of the coefficients is shown in Fig. 10 for model I and in Fig. 11 for model II.

Rescaling the time $`T=\epsilon t`$ and amplitude $`A=\sqrt{\epsilon }B`$ back to the original units yields the amplitude equation

$`\tau _0_tA=\epsilon \left(1+ia\right)Ag\left(1+ic\right)|A|^2A,`$ (134)
as introduced in Sec. I. This equation has simple nonlinear oscillatory solutions of the form $`A=Fe^{i\mathrm{\Omega }t}`$, with an amplitude $`F`$ and a frequency $`\mathrm{\Omega }`$ as follows

$`F`$ $`=`$ $`\sqrt{{\displaystyle \frac{\epsilon }{g}}},`$ (135)
$`\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \frac{1}{\tau _0}}\left(\epsilon agcF^2\right)={\displaystyle \frac{ac}{\tau _0}}\epsilon .`$ (136)
$`\mathrm{\Omega }`$ describes the deviation of the oscillation frequency from the critical one, $`\omega _c`$.

The linear coefficients $`\tau _0`$ and $`a`$ of Eq. (134) may directly be calculated from the dispersion relation in Eq. (50) or Eq. (92) in the following way. The solution $`A=0`$ of Eq. (134) corresponds to the stationary polymerization described in Sec. III, which is in the range $`\epsilon <0`$ stable against small perturbations $`A\stackrel{~}{F}e^{\sigma t}`$ ( with $`\stackrel{~}{F}\sqrt{|\epsilon |}`$) and unstable for $`\epsilon >0`$. Neglecting in Eq. (134) the contributions due to the cubic nonlinearity one obtains from its linear part the dispersion relation $`\sigma =\epsilon (1+ia)/\tau _0`$. This formula gives the relaxation time $`\tau _0`$ and the linear frequency dispersion $`a`$ in terms of derivatives of the growth rate with respect to the control parameter $`\epsilon `$: $`\tau _0=(Re(\sigma )/\epsilon )^1`$ and $`a=\tau _0Re(\sigma )/\epsilon `$. If $`\epsilon `$ is expressed in terms of the dimer density $`c_0`$, cf. Eq. (129), then both quantities may also be written in terms of the derivatives with respect to $`c_0`$

$`\tau _0={\displaystyle \frac{1}{c_{0c}Re(\sigma )/c_0}},a=c_{0c}\tau _0{\displaystyle \frac{Im(\sigma )}{c_0}},`$ (137)
whereby both derivatives are taken at the threshold concentration $`c_{0c}`$. The dispersion relation $`\sigma (\epsilon )`$ as obtained on the one hand by the amplitude equation and on the other hand by solving Eq. (50) or Eq. (92), both have to reproduce the growth or decay dynamics of small perturbations with respect to the stationary polymerization. Therefore, the coefficients $`\tau _0`$ and $`a`$ of the amplitude equation can directly be calculated via Eq. (137) from the numerical solutions of Eq. (50) or Eq. (92).

One aim of the amplitude expansion is the determination of the type of the Hopf bifurcation. For two different nucleation rates $`\nu =0.01`$ and $`\nu =0.04`$ the variation of the nonlinear coefficient $`g`$ as function of the regeneration rate $`\alpha `$ is shown for model I in Fig. 10 (top) and for model II with oligomer dynamics in Fig. 11. In both cases $`g`$ behaves rather similar and $`g`$ is positive for the models investigated in this work. Therefore the Hopf bifurcation is supercritical. In Fig. 10 and Fig. 11 the nonlinear coefficient $`g`$ increases at first with the regeneration rate $`\alpha `$ and reaches a maximum in a range where the threshold concentration $`c_{0c}(\alpha )`$ takes its minimum. The corresponding threshold curves $`c_{0c}(\alpha )`$ for two different nucleation rates $`\nu =0.01`$ and $`\nu =0.04`$ also cross each other. Beyond this maximum of $`g`$ ( where $`c_{0c}(\alpha )`$ takes its minimum) decreases again.

It should be mentioned that for a given value of the control parameter $`\epsilon `$ a large value of $`g`$ corresponds to a small value of the oscillation amplitude. Since the threshold $`c_{0c}(\alpha )`$ varies too, the variation of the oscillation amplitude with $`\alpha `$ is much less when $`c_{0c}\epsilon `$ is kept constant. According to the sign of the nonlinear coefficient $`c`$ the oscillation frequency $`\omega _c+\mathrm{\Omega }`$ decreases with increasing values of $`\epsilon `$. The results in terms of the amplitude equation are in fairly good agreement with the behavior of the full numerical solution of the basic reaction equations.

A determination of the bifurcation structure by numerical simulations of the basic equations is error prone compared to results of perturbation calculation described here. Besides the lower accuracy, parameter studies such as in Fig. 10 and Fig. 11 are with numerical simulations much more time consumptive.

Close to threshold the advantages of the perturbation calculation are obvious. However, it is a priori unknown in which $`\epsilon `$-range the amplitude equation approach (134) applies quantitatively. For some systems the amplitude equation is valid in a rather large range of the control parameter $`\epsilon `$, but for other systems its validity is restricted to very small values of it, cf. Ref. . In order to check this for our models of microtubule polymerization, we compare in Fig. 12 the variation of the oscillation amplitude of $`c_t^{(1)}`$ with the control parameter $`\epsilon `$ as obtained by the numerical solution described in Sec. V and by the solution $`A=\sqrt{\epsilon /g}`$ of the amplitude equation for two different values of the nucleation rate $`\nu `$. At larger values of the control parameter, $`\epsilon =0.1`$, the difference between the results for the numerical solution and the amplitude equation is still less than $`8\%`$. For model II the deviations are larger between the amplitude determined by the amplitude equations approach and the numerical solutions with ansatz in Eq. (V). However, when we solve the equation for growing microtubule, cf. Eq. (2), numerically by discretization of the length coordinate, the deviations becomes smaller.

For both models, the linear coefficient $`\tau _0`$ and the nonlinear coefficients $`g,c`$ don’t differ very much from each other. However, the linear frequency shift $`a`$ increases in model II for large regeneration rates $`\alpha `$ whereas for model I it decreases. Nevertheless in both models the nonlinear frequency correction due to the values of $`c`$ is much larger than the linear correction due to $`a`$, $`c>>a`$.

Whether the bifurcation to oscillatory polymerization is also supercritical in experiments under stationary regeneration conditions, is an open question. Therefore an experimental determination of the bifurcation–type would be an important test for the reduced models investigated in this work.

The variation of the other coefficients with the regenerate $`\alpha `$ provides further contact between the model parameters and experimentally measurable quantities. The parameters $`c`$, $`a`$ and $`\tau _0`$ may be determined as follows. Studying the growth of small perturbations, $`c_t^{(1)}e^{Re(\sigma )t}=e^{\epsilon t/\tau _0}`$ by plotting the logarithm $`\epsilon t/\tau _0\mathrm{log}(c_t^{(1)})`$ as function of time and for different values of $`\epsilon `$, the relaxation time $`\tau _0`$ may be determined. In a similar manner $`a`$ may be determined by studying the frequency of a perturbation far below its nonlinear saturation amplitude. If the perturbation saturates finally, the oscillation frequency of the nonlinear solution changes with $`\epsilon `$ as indicated by Eq. (135). From this $`\epsilon `$–dependence the nonlinear coefficient $`c`$ may be extracted. An experimental determination of these coefficients, as described, would be a further test of the basic model equations.

## VII Summary and Conclusion

Two reduced models, that capture oscillating microtubule polymerization and the length distribution of the microtubule filaments, as described in Sec. II, have been analyzed. In both models the complex biochemical reaction steps of microtubule polymerization are described by a few ones, which have been identified in experiments to be important. The focus on a few essential degrees of freedom leads to some simplicity of the models that allows for instance a derivation of analytical expressions for the threshold and the oscillation frequency at the Hopf bifurcation. Such analytical results make trends as function of the reaction rates more easily visible. Some of these trends may be tested in experiments and some of the reaction constants may be measured.

At threshold also analytical expressions could be derived for the temporal evolution of the concentrations and the length distribution of microtubules. Those provide a detailed picture about the temporal variation of the fields, their relative phases and the amplitude ratios between them. The formula for the length distribution is especially instructive, cf. Eq. (70), it describes a superposition of amplitude oscillations of the distribution and traveling waves, where the waves always travel towards larger lengths. This qualitative behavior of the length distribution during oscillatory polymerization is rather independent of the respective model and it is a rather general feature. A few snap shots of the numerically generated distribution far beyond threshold are shown in Fig. 9b). The distribution includes still traveling waves, but far beyond threshold these behave rather anharmonic.

For stationary reaction conditions, as assumed in this work, Fig. 9 shows a remarkable result. In part a) and b) of this figure a subthreshold concentration for tubulin was assumed, i.e. $`c_0<c_{0c}`$, and in both cases the final state is stationary polymerization. However, on the route to the stationary state transient oscillations occur. Therefore, the transient character of the microtubule oscillations observed in an experiment with an enzymatic regeneration process for GTP , could have, according to the results described in this work, its origin in a to low tubulin concentration. A higher tubulin concentration or an appropriate regeneration rate of GTP and a different live time of oligomers could lead in a similar experiment to persistent microtubule oscillations.

For transient oscillations observed in other experiments the common interpretation is as follows. During the microtubule polymerization in experiments the available GTP is used up and the oscillations last only for a few periods. If GTP is continuously supplied, the simultaneously increasing amount of GTD inhibits various reactions steps and slows down the reaction cycle. The results shown in Fig. 9 a) and Fig. 9 b) indicate that oscillations may occur as a transient because either GTP is used up or the tubulin concentration has only a subthreshold value. Accordingly, there may be several reasons for transient oscillations in experiments. Either the initial tubulin concentration has a subthreshold value, the decay rate of oligomers and the regeneration rate for GTP do not have their optimal values which would explain transient oscillations in experiments with a regenerative enzyme system, the reaction conditions are not constant, because GTP is used up. In the latter case the available tubulin–t decreases with time and the microtubule polymerization decays.

In an in vitro experiment with constant reaction conditions the threshold of oscillations might be measured by increasing the tubulin dimer concentration by appropriate steps. Immediately after each step transient oscillations might occur, but the threshold is only crossed when the oscillations persist over a long time.

In order to avoid numerical instabilities during long time simulations of the reaction equations, including Eq. (8), we use analytical approximations for the length dependence of the microtubule distributions as described in Sec. V. This stable numerical scheme can be generalized in future work to an effective algorithm for dealing with microtubule polymerization in one and two spatial dimensions  in order to investigate spatial patterns occurring during polymerization of microtubules . The respective extension of the amplitude equation may also lead to new interesting insights.

Microtubule filaments at a high density show a isotropic–nematic phase transition , similar as it has been observed for F-actin filaments . Since the early theory of Onsager  this transition for rods in a solvent is a well understood phenomenon . For a monodisperse distribution of filaments of fixed length, i.e. without polymerization kinetics, many aspects of the isotropic–nematic transition have been understood. For polydisperse rod mixtures some aspects of the isotropic-nematic transition can be considered to be understood too . The effects of nucleation, growth of filaments and the decay of filaments on the isotropic–nematic transition are not known at present and one may expect interesting phenomena related to this kinetics . The effect of oscillating microtubule polymerization on the isotropic–nematic transition is also completely unexplored at present and will be investigated in forthcoming works .

It is a great pleasure to thank M. Breidenich, H. Flyvbjerg, E. Mandelkow, H. Müller–Krumbhaar, J. Prost and J. Tabony for interesting discussions.

## A Amplitude expansion for model I

For model I and the catastrophe rate given in Eq. (3) the major steps of the derivation of the amplitude equation (1) are described in this appendix. Since we neglect the rescue of shrinking microtubules, $`f_{resc}=0`$, the only nonlinear term in the basic equations of model I is the product $`f_{cat}(c_t)p_g`$ in Eq. (II B). At first we expand the concentrations $`c_t,c_0`$ and length distributions $`p_{g,s}`$ with respect to deviations from their stationary value at the threshold for oscillation, such as for instance for the tubulin-t concentration, i.e. $`c_tc_t^{(0)}=\sqrt{\epsilon }\stackrel{~}{c}_t^{(1)}+\epsilon \stackrel{~}{c}_t^{(2)}+\mathrm{}`$. Note that the tilded fields differ just by a power of $`\sqrt{\epsilon }`$ from the untilded fields as introduced in Sec. IV, cf. $`(\sqrt{\epsilon })^j\stackrel{~}{c}_t^{(j)}=c_t^{(j)}`$ etc. . In order to simplify the notation of this appendix we drop the “tilde” and the expansion of the catastrophe rate takes the form

$`f_{cat}=f_{cat}^{(0)}+\epsilon ^{1/2}f_{cat}^{(1)}+\epsilon f_{cat}^{(2)}+\epsilon ^{3/2}f_{cat}^{(3)}+\mathrm{},`$ (A1)
whereby the coefficients of this expansion are

$`f_{cat}^{(1)}`$ $`=`$ $`f_{cat}^{(0)}{\displaystyle \frac{c_t^{(1)}}{c_f}},`$ (A3)
$`f_{cat}^{(2)}`$ $`=`$ $`f_{cat}^{(0)}\left({\displaystyle \frac{1}{2}}\left({\displaystyle \frac{c_t^{(1)}}{c_f}}\right)^2{\displaystyle \frac{c_t^{(2)}}{c_f}}\right),`$ (A4)
$`f_{cat}^{(3)}`$ $`=`$ $`f_{cat}^{(0)}\left({\displaystyle \frac{c_t^{(1)}c_t^{(2)}}{c_f^2}}{\displaystyle \frac{c_t^{(3)}}{c_f}}{\displaystyle \frac{1}{6}}\left({\displaystyle \frac{c_t^{(1)}}{c_f}}\right)^3\right).`$ (A5)
Collecting in Eq. (II B) and Eq. (15) the contributions to the order $`\epsilon ^{1/2}`$ we recover the linear equations given in Sec. IV

$`_tp_g^{(1)}`$ $`=`$ $`f_{cat}^{(0)}{\displaystyle \frac{c_t^{(1)}}{c_f}}p_g^{(0)}f_{cat}^{(0)}p_g^{(1)}v_g_lp_g^{(1)},`$ (A7)
$`_tp_s^{(1)}`$ $`=`$ $`f_{cat}^{(0)}{\displaystyle \frac{c_t^{(1)}}{c_f}}p_g^{(0)}+f_{cat}^{(0)}p_g^{(1)}+v_s_lp_s^{(1)},`$ (A8)
$`_tc_t^{(1)}`$ $`=`$ $`\gamma {\displaystyle _0^{\mathrm{}}}𝑑l\left(v_gp_g^{(1)}+\alpha l(p_g^{(1)}+p_s^{(1)})\right)`$ (A10)
$`\alpha c_t^{(1)}.`$
At order $`\epsilon `$ we obtain the three equations

$`_tp_g^{(2)}`$ $`=`$ $`\left(f_{cat}^{(0)}{\displaystyle \frac{c_t^{(2)}}{c_f}}p_g^{(0)}f_{cat}^{(0)}p_g^{(2)}v_g_lp_g^{(2)}\right)`$ (A13)
$`+f_{cat}^{(0)}{\displaystyle \frac{c_t^{(1)}}{c_f}}p_g^{(1)}{\displaystyle \frac{f_{cat}^{(0)}}{2}}\left({\displaystyle \frac{c_t^{(1)}}{c_f}}\right)^2p_g^{(0)},`$
$`_tp_s^{(2)}`$ $`=`$ $`\left(f_{cat}^{(0)}{\displaystyle \frac{c_t^{(2)}}{c_f}}p_g^{(0)}+f_{cat}^{(0)}p_g^{(2)}+v_s_lp_s^{(2)}\right)`$ (A15)
$`f_{cat}^{(0)}{\displaystyle \frac{c_t^{(1)}}{c_f}}p_g^{(1)}+{\displaystyle \frac{f_{cat}^{(0)}}{2}}\left({\displaystyle \frac{c_t^{(1)}}{c_f}}\right)^2p_g^{(0)},`$
$`_tc_t^{(2)}`$ $`=`$ $`\gamma {\displaystyle _0^{\mathrm{}}}𝑑l\left(v_gp_g^{(2)}+\alpha l(p_g^{(2)}+p_s^{(2)})\right)`$ (A17)
$`+\alpha c_{0c}\alpha c_t^{(2)}.`$
With the solutions of the equations at the previous order $`\epsilon ^{1/2}`$, that are already given in Sec. IV A 3, the equations at order $`\epsilon `$ have to be solved. These solutions have the following form

$`c_t^{(2)}`$ $`=`$ $`A_0+A_2\mathrm{exp}(2i\omega _ct)+c.c.,`$ (A19)
$`p_g^{(2)}`$ $`=`$ $`e^{f_{cat}^{(0)}l/v_g}(B_0(l)+[B_2(l)e^{2i\omega _ct}+c.c.]),`$ (A20)
$`p_s^{(2)}`$ $`=`$ $`e^{f_{cat}^{(0)}l/v_g}(D_0(l)+[D_2(l)e^{2i\omega _ct}+c.c.]),`$ (A21)
whereby the expressions for the coefficients $`A_0`$, $`A_2`$, $`B_i`$ and $`D_i`$ are rather lengthy in terms of the coefficients of the solutions at order $`\epsilon ^{1/2}`$ and are not given here.

The equations at the next higher order $`\epsilon ^{3/2}`$ are

$`_Tp_g^{(1)}+_tp_g^{(3)}`$ $`=`$ $`\left(f_{cat}^{(0)}{\displaystyle \frac{c_t^{(3)}}{c_f}}p_g^{(0)}f_{cat}^{(0)}p_g^{(3)}v_g_lp_g^{(3)}\right)`$ (A26)
$`\left(f_{cat}^{(0)}{\displaystyle \frac{c_t^{(1)}c_t^{(2)}}{c_f^2}}{\displaystyle \frac{f_{cat}^{(0)}}{6}}\left({\displaystyle \frac{c_t^{(1)}}{c_f}}\right)^3\right)p_g^{(0)}`$
$`+\left(f_{cat}^{(0)}{\displaystyle \frac{c_t^{(2)}}{c_f}}{\displaystyle \frac{f_{cat}^{(0)}}{2}}\left({\displaystyle \frac{c_t^{(1)}}{c_f}}\right)^2\right)p_g^{(1)}`$
$`+f_{cat}^{(0)}{\displaystyle \frac{c_t^{(1)}}{c_f}}p_g^{(2)},`$
$`_Tp_s^{(1)}+_tp_s^{(3)}`$ $`=`$ $`\left(f_{cat}^{(0)}{\displaystyle \frac{c_t^{(3)}}{c_f}}p_g^{(0)}+f_{cat}^{(0)}p_g^{(3)}+v_s_lp_s^{(3)}\right)`$ (A30)
$`+\left(f_{cat}^{(0)}{\displaystyle \frac{c_t^{(1)}c_t^{(2)}}{c_f^2}}{\displaystyle \frac{f_{cat}^{(0)}}{6}}\left({\displaystyle \frac{c_t^{(1)}}{c_f}}\right)^3\right)p_g^{(0)}`$
$`\left(f_{cat}^{(0)}{\displaystyle \frac{c_t^{(2)}}{c_f}}{\displaystyle \frac{f_{cat}^{(0)}}{2}}\left({\displaystyle \frac{c_t^{(1)}}{c_f}}\right)^2\right)p_g^{(1)}`$
$`f_{cat}^{(0)}{\displaystyle \frac{c_t^{(1)}}{c_f}}p_g^{(2)},`$
$`_Tc_t^{(1)}+_tc_t^{(3)}`$ $`=`$ $`\gamma {\displaystyle _0^{\mathrm{}}}𝑑l\left(v_gp_g^{(3)}+\alpha l(p_g^{(3)}+p_s^{(3)})\right)`$ (A32)
$`\alpha c_t^{(3)}.`$
The two fields $`p_g^{(3)}`$ and $`p_s^{(3)}`$ have to be calculated explicitly at this order from Eq. (A26) and Eq. (A30) as well. With both solutions the integral on the right hand side of Eq. (A32) can be calculated. Eq. (A26) and Eq. (A30) include both contributions proportional to $`e^{i\omega _ct}`$ and $`e^{3i\omega _ct}`$, but only the single harmonic terms are relevant in Eq. (A32). The coefficient of $`e^{i\omega _ct}`$ in Eq. (A32) must vanish. Part of it vanishes automatically, because it reproduces the threshold condition and the rest provides the amplitude equation with all the coefficients now given in terms of the reaction rates of the basic equations.