Unifying autocatalytic and zeroth order
branching models for growing actin networks

In many situations of high biological relevance, including the migration of animal cells and the propulsion of specific intracellular pathogens, motility results from the directed polymerization of a dendritic actin filament network Carlier (2010). The organization of the growing network is determined mainly at the leading edge, where a small number of proteins regulates the interplay between three fundamental processes. The driving force for propulsion is polymerization of actin filaments from globular actin monomers. This is limited by capping proteins, which bind to the filament ends and prevent further polymerization. New filaments nucleate by branching off from mother filaments Pollard (2007). Although the biochemical details of this process are not yet completely understood, it is widely accepted that the branching complex Arp2/3 is activated by nucleation promoting factors (NPFs) like WASP and SCAR/WAVE proteins Beltzner and Pollard (2008); Xu et al. (2012). When an activated Arp2/3-complex is bound to the side of an existing actin filament, a daughter filament starts to grow at a characteristic angle around $70^{\circ}$ relative to the mother filament (compare Fig. 1a). At the same time, the branch point moves away from the leading edge because of the on-going polymerization of actin filaments.

Due to the high biological relevance and universal nature of the underlying processes, many theoretical models have been suggested to describe the characteristic features of growing actin networks Mogilner (2009). However, in many cases contradictory predictions have been obtained, in particular regarding experimentally observed force-velocity relations Wiesner et al. (2003); McGrath et al. (2003); Marcy et al. (2004); Parekh et al. (2005); Prass et al. (2006); Heinemann et al. (2011); Zimmermann et al. (2012) and the filament orientation distribution of the network Maly and Borisy (2001); Verkhovsky et al. (2003); Schaub et al. (2007); Koestler et al. (2008); Weichsel et al. (2012). Interestingly, many of these contradictions are a direct consequence of two different choices for the order of the branching reaction. In autocatalytic models, the branching rate is assumed to be proportional to the number of existing filaments in the network, i.e. it is modeled as a first order reaction in filament density, implicitly assuming an unlimited reservoir of activated Arp2/3 Maly and Borisy (2001); Carlsson (2003); Schaus et al. (2007). This yields growing actin networks for which a constant filament density is maintained only at a unique steady state growth velocity. Increasing forces acting against the network reduce the speed of growth only transiently, as an increasing filament density subsequently lowers the force per filament back to the stationary level.

In marked contrast to the autocatalytic scenario, another class of models assumes that branching occurs with a constant rate, i.e. it is taken to be a zeroth order reaction in filament density, corresponding to a limited supply of activated Arp2/3 Carlsson (2003); Schaub et al. (2007); Weichsel and Schwarz (2010). Under these conditions, it has been shown that a continuum of steady state velocities exists. Moreover, two competing steady state filament orientation patterns are stable, namely the $\pm 35$ and $+70/0/\!\!-\!\!70$ patterns shown schematically in Fig. 1b and c, respectively. Transitions between these two fundamentally different network architectures can be triggered by changes in network growth velocity Maly and Borisy (2001); Weichsel and Schwarz (2010). Indeed similar structural transitions have been demonstrated recently in electron microscopy data of the lamellipodium of keratocytes, indicating their physiological relevance Koestler et al. (2008); Weichsel et al. (2012). In this Letter, we will show that the two contradictory model scenarios of autocatalytic and zeroth order branching can be unified within a general theoretical framework that reconciles some of the seemingly contradictory observations and predictions.

Arp2/3 activation model. We first introduce a kinetic model for filament branching, based on a likely scenario for Arp2/3 activation Beltzner and Pollard (2008); Ti et al. (2011); Xu et al. (2012). Motivated by the dimensions of the lamellipodium for cells migrating on a flat substrate, we consider a two-dimensional situation in which the network moves away from the leading edge with a well defined retrograde velocity $v_{\rm nw}$. All reactions are assumed to occur in a small reaction zone extending from the leading edge over a nanometer-scale distance $d_{\rm br}$. We consider a system of two variables: $A$ is the concentration of Arp2/3 that is bound to the filaments, but did not lead to a daughter branch yet. $P$ is the concentration of NPFs which is available to activate bound Arp2/3 complexes to nucleate a daughter branch. The kinetic equations are

$A$ increases as more complexes bind to the filaments with rate $k_{+}$ and decreases due to dissociation (rate $k_{-}$), outgrowth (rate $v_{\rm nw}/d_{\rm br}$) and branching (rate $\tilde{k}_{\rm b}$). The last step also decreases available $P$, as NPFs that activate Arp2/3 are occupied for additional interactions with other Arp2/3 complexes at the same time until they become available again at rate $k_{\rm act}$. $P_{\rm 0}$ is the total concentration of NPFs and $N_{\rm fil}$ is the number of actin filaments (because $N_{\rm fil}$ will be a central quantity of interest below, for our purpose it is convenient to consider the number of filaments in a reaction volume of finite lateral size rather than their concentration).

In steady state, Eq. (1) defines an effective rate of branching as $\mathrm{BR}^{\rm ss}\equiv\tilde{k}_{\rm b}A_{\rm ss}P_{\rm ss}$. This rate is a function of filament number $N_{\rm fil}$ as plotted in Fig. 2 for a typical set of parameters. At sufficiently small filament number $N_{\rm fil}$, the effective branching rate is approximately first order in $N_{\rm fil}$ and hence the coefficient of its linear expansion defines an autocatalytic branching rate constant $k_{\rm b}^{\rm ac}$:

In the limit of large filament number $N_{\rm fil}$, $\mathrm{BR}^{\rm ss}$ saturates at a constant rate as assumed in zeroth order branching models:

Thus the reaction smoothly changes from first to zeroth order as the filament number $N_{\rm fil}$ increases.

Actin growth model. We next analyze the effect of the order of the branching reaction on the steady growth states of actin networks. To this end, we extend a deterministic rate equation model that has been used before to describe both autocatalytic as well as zeroth order branching actin networks Maly and Borisy (2001); Carlsson (2003); Weichsel and Schwarz (2010). The generic results reported here can be confirmed in computer simulations based on individual filaments and stochastic reactions si (2012). We consider an ensemble of filaments located in the same reaction zone of width $d_{\rm br}$ as introduced above. Our central quantity is the distribution function $N(\theta,t)$ for the number of uncapped filaments orientated at time $t$ at an angle $\theta$ with respect to the normal of the leading edge, which evolves in time as

Here the three terms on the right introduce capping, outgrowth from the reaction zone and branching, respectively. While capping is simply a first order process with constant rate, independent of filament orientation $\theta$, for outgrowth we have to distinguish two cases. For $|\theta|\leq\arccos(v_{\rm nw}/v_{\rm fil})$, single filaments growing with velocity $v_{\rm fil}$ can keep up with the leading edge and thus $k_{\rm gr}^{\theta}(v_{\rm nw})=0$. If the orientation angle exceeds the threshold, filaments grow too slowly and leave the reaction region with rate $k_{\rm gr}^{\theta}(v_{\rm nw})=(v_{\rm nw}-v_{\rm fil}\cos{\theta})/d_{\rm br}$.

In the branching term, $\mathcal{W}(\theta,\theta^{\prime})=\mathcal{W}(|\theta-\theta^{\prime}|)$ is a distribution function of the relative branching angle between mother and daughter filaments. Motivated by experimental observations, we approximate this function by the sum of two Gaussians centered around $\pm 70^{\circ}$ and each with standard deviation $5^{\circ}$. This corresponds well to the experimentally reported range of branching angles between $67^{\circ}$ and $77^{\circ}$, that have been measured both using purified proteins Mullins et al. (1998); Blanchoin et al. (2000) and in different cell lines Svitkina and Borisy (1999); Vinzenz et al. (2012). The exact value of the branching angle, however, is irrelevant for our results.

The normalization of the branching term in Eq. (4) is appropriate to directly implement a specific reaction order $\mu$ in the actin growth model. In order to couple the actin growth model to the Arp2/3 activation model, below we will use numerical calculations with $\mu=0$ and couple Eq. (1) and Eq. (4) via the filament number dependent branching rate $k_{\rm b}(N_{\rm fil})=\mathrm{BR}^{\rm ss}$ derived above, with $N_{\rm fil}=\int N(\theta,t)d\theta$. For analytical progress and deeper insight, however, it is instructive to first analyze the steady states of the actin growth model with constant parameters $k_{\rm b}$ and $\mu$. Indeed, it can be shown within the linear stability analysis employed below that both procedures are equivalent si (2012).

Steady state analysis. The actin growth model Eq. (4) contains only four relevant parameters, the rates $k_{\rm c}$ and $k_{\rm b}$ for capping and branching, respectively, the network growth velocity $v_{\rm nw}$ and the order of the branching reaction $\mu$. By integrating the reaction model Eq. (4) over $35^{\circ}$ sized angle bins and neglecting contributions from filaments growing in directions $>87.5^{\circ}$, we obtain three simplified coupled equations for the evolution of $N_{0^{\circ}}$, $N_{\pm 35^{\circ}}$ and $N_{\pm 70^{\circ}}$, which can be analyzed analytically. As an alternative which does not require any additional assumptions, we determine the stable regimes of network growth numerically, by propagating a finely discretized version of the equation until a steady state is reached.

In the analytical approach, there exist exactly two physically meaningful steady state solutions, $N^{\rm ss70}$ and $N^{\rm ss35}$, given by

and

where

These two fixed points correspond to the two competing orientation patterns depicted schematically in Fig. 1c and b, respectively. Linear stability analysis shows that for $\mu>1$, both are saddle points and thus no stable solution exists. In contrast, $\mu\leq 1$ leads to mutually exclusive stability of the two solutions si (2012). Fig. 3a shows the regions of stability for each of the two orientation patterns within the two dimensional parameter space spanned by $k_{\rm c}$ and $v_{\rm nw}$. The dashed contour indicates transitions between a $+70/0/\!\!-\!\!70$ pattern outside and a $\pm 35$ pattern inside. Remarkably, the transition is independent of $k_{\rm b}$ and $\mu$ and thus all cases with $\mu\leq 1$ show no difference in the locations of the transitions. The result from the full numerical analysis of Eq. (4) is shown as inset. The main difference between the analytical and the numerical result is that the stability of the $\pm 35$ pattern vanishes for large $k_{\rm c}$ in the analytical model, because it disregards contributions from filament orientations $\gtrsim 90^{\circ}$. This increases stability of the $+70/0/\!\!-\!\!70$ pattern, when outgrowth of filaments is negligible compared to capping.

Until now we have shown, that for all reaction orders of interest ($0\leq\mu\leq 1$), two orientation patterns compete for stability, with phase boundaries being independent of the exact value of $\mu$. Nevertheless the limit $\mu\to 1$ (autocatalytic growth) is special, because in this case, finite steady state solutions only exist if additional constraints are satisfied. Due to Eq. (5), $N^{\rm ss70}$ is finite only when $C_{1}=1$. From Eq. (7) this requires $k_{\rm gr}^{70^{\circ}}(v_{\rm nw})k_{\rm c}+k_{\rm c}^{2}=k_{\rm b}^{2}/2$. Due to Eq. (6), $N^{\rm ss35}$ is finite only when $C_{2}=1$. From Eq. (7) this requires $k_{\rm gr}^{35^{\circ}}(v_{\rm nw})+k_{\rm c}=k_{\rm b}/2$. Therefore, if a stable steady state solution exists for given values of $k_{\rm c}$ and $k_{\rm b}$, then for $\mu=1$ it corresponds to a unique network growth velocity $v_{\rm nw}$. In this way, the most prominent feature of autocatalytic growth Carlsson (2003) emerges in our unifying model. Due to these additional conditions, stable solutions are restricted in parameter space to a lower dimensional manifold, which in case of the analytical model has a jump discontinuity si (2012). In Fig. 3a we show the projection of this manifold onto the $(k_{\rm c},v_{\rm nw})$-plane with bright and dark gray regions marking the stability regions for the $+70/0/\!\!-\!\!70$ and $\pm 35$ patterns, respectively. As the inset indicates, a jump discontinuity is not observed in the full numerical treatment.

In Fig. 3b, the manifold for $\mu=1$ is projected onto the $(k_{\rm c},k_{\rm b})$-plane (in this projection, the jump discontinuity cannot be seen). Only a subset of ($k_{\rm b},k_{\rm c})$-combinations yields stable autocatalytic growth with finite filament number $N_{\rm fil}$. This important result is predicted both by the analytical and the numerical approach (compare inset).

Limits of autocatalytic network growth. Using the insights obtained in the preceding sections from the actin growth model Eq. (4) with $\mu$ as a model parameter, we now combine the Arp2/3 activation model Eq. (1) and the actin growth model Eq. (4) with $\mu=0$ to arrive at a unifying theoretical framework for actin network growth with a branching reaction that is determined by a regulatory process. Fig. 4 shows our numerical results for network growth velocity $v_{\rm nw}$ as a function of filament number $N_{\rm fil}$ (solid lines) for various values of the capping rate $k_{\rm c}$. They agree very well with the results from stochastic computer simulations shown as inset si (2012). At sufficiently low filament number, we observe an autocatalytic regime where a whole range of values for $N_{\rm fil}$ corresponds to the same velocity. For larger $N_{\rm fil}$, however, the steady state network velocity starts to decrease similarly to a pure zeroth order description (dashed lines). These changes include transitions between the two dominant filament orientation patterns as predicted in Fig. 3.

Interestingly, the details of the crossover from first to zeroth order branching strongly depend on the capping rate. This can be understood in the analytical model analyzed above. At low filament density, branching is effectively a first order reaction and thus the conditions, $C_{1}=1$ and $C_{2}=1$, previously derived from Eq. (5)–Eq. (7) for $\mu=1$, need to be satisfied here for stable growth as well. By inserting the autocatalytic branching rate $k_{\rm b}^{\rm ac}$ defined in Eq. (2) into Eq. (7) and applying the relevant first order condition, we are able to derive estimates for minimum and maximum capping rates $k_{\rm c}^{\rm min}$ and $k_{\rm c}^{\rm max}$ corresponding to the largest and smallest possible network velocities, $v_{\rm nw}/v_{\rm fil}=1$ and $v_{\rm nw}/v_{\rm fil}=0$, respectively:

where $c=\cos 70^{\circ}-1$. These two threshold values are shown in Fig. 3b as the intersection of $k_{\rm b}^{\rm ac}$ with the boundaries of the stable autocatalytic parameter subset. Comparison of Eq. (8) with the numerical results from the full model presented in Fig. 4 shows that our analytical approach captures the location of this crossover very well and thus accurately explains the observed behavior. For increasing $k_{\rm c}$, the network velocity in the autocatalytic region decreases until at around $k_{\rm c}\simeq k_{\rm c}^{\rm max}$ the filament number decays to zero for all accessible network velocities. For decreasing $k_{\rm c}$, the network growth velocity reaches its maximal value at $k_{\rm c}\lesssim k_{\rm c}^{\rm min}$ (thick solid line), when the network is not able to balance filament branching by capping and outgrowth anymore. In a purely autocatalytic model, this would lead to a diverging and therefore unphysical filament number. Within our unifying framework, the number of filaments increases only to the point where zeroth order growth behavior starts to dominate and stabilizes a steady state at finite filament number. In this regime, the results from the full model (solid) agree with a zeroth order branching model (dashed).

Relation to experiments. In this Letter, we have developed a theoretical framework that reconciles conflicting results from two classes of actin growth models and explains many experimental observation: an autocatalytic growth regime at low filament density Wiesner et al. (2003); Parekh et al. (2005), zeroth order characteristics at high density McGrath et al. (2003); Marcy et al. (2004), network velocity-dependent transitions in filament orientation patterns Koestler et al. (2008); Weichsel et al. (2012) and bistability and hysteresis at these transitions Parekh et al. (2005); Weichsel and Schwarz (2010). Strikingly our model naturally avoids the instability which occurs at low capping rate in the autocatalytic model.

Our model also makes testable predictions that can guide future experiments. Using single molecule microscopy either in migrating cells Millius et al. (2012) or in reconstituted assays, the number of branching events can be directly correlated to filament density, which can be compared to the effective branching rate as predicted in Fig. 2. From electron microscopy data, filament orientations can be extracted and correlated with the growth velocity as demonstrated in Koestler et al. (2008); Weichsel et al. (2012). This can be compared to the unified phase diagram in Fig. 3. Force-velocity relations could be calculated for our model along the lines of Refs. Carlsson (2003); Weichsel and Schwarz (2010), but would require additional assumptions regarding for example network mechanics, load sharing and filament-membrane interactions. However, some general conclusions can already be drawn at this point and are explained best for the case of reconstituted actin networks growing against a functionalized AFM cantilever or bead Marcy et al. (2004); Parekh et al. (2005); Chaudhuri et al. (2009). In this context, our model predicts an autocatalytic (i.e. force-insensitive) growth velocity for sufficiently low load and high concentration of capping protein. In this regime the filament density near the obstacle is thus expected to grow proportional to the applied force. When either the concentration of capping protein is reduced below the threshold $k_{\rm c}^{\rm min}$ (Eq. (8)) or the load on the network is sufficiently increased, zeroth order behavior is predicted to take over as is illustrated in Fig. 4. If combined with mechanical models like Ref. Zimmermann et al. (2012), in the future these kinetic considerations might lead to a complete understanding of the intriguing physics of growing actin networks.