Low-dimensional manifold of actin polymerization dynamics

The BC model of actin polymerization is an 11-dimensional system of ODEs tracking the concentrations of non-tip actin subunits in different states as well as the concentrations of filament tip subunits in different states [20]. It is assumed that the number concentration of filaments N remains constant, implying an absence of filament nucleation, splitting, or joining. Additionally, the total concentration of actin subunits M is assumed to remain constant, such that actin subunits are not created or destroyed in any reaction. Since there are typically many actin subunits belonging to a given actin filament, we have $N\ll M$. All species are assumed to be well-mixed and in large enough quantities to be treated effectively via a mean-field description. In other words, the size of the stochastic fluctuations is negligible compared to the concentrations of the species. Unpolymerized (globular) actin subunits are referred to as G-actin, while polymerized (filamentous) actin subunits are referred to as F-actin. Actin filaments are helical, but they are more easily modeled as linear chains, which is a realistic approximation if one assumes that the reaction propensities of a given F-actin subunit are determined only by the state of the nucleotide bound by that subunit and not by the subunit's neighbors. Such a chain is displayed in figure 1, along with some of the reactions allowing interconversion between subunit types. The variables representing the concentrations of these actin species are superscripted by the hydrolysis state of the bound nucleotide (for what follows we refer more simply to a subunit being in a certain hydrolysis state as opposed to the nucleotide attached to a subunit as being in that state). The hydrolysis states are ATP, ADP-Pi, and ADP, denoted T, Pi, and D, respectively. The tip subunits are denoted T and are further subscripted according to which tip they are on. Thus, for example, the concentration of tip subunits at the plus end bound to ADP-Pi is denoted ${{T}}_{+}^{{\rm{Pi}}}$. Because inorganic phosphate rapidly dissociates from G-actin, ${G}^{{Pi}}$ is taken to be 0 and is not tracked. With the 3 hydrolysis states of each of the 2 tips, the 3 states of the F-actin and the 2 states of G-actin, the number of tracked variables is 11.

Zoom In Zoom Out Reset image size

The different subunit types interconvert through chemical reactions. These reactions can be classified as polymerization/depolymerization reactions which change G actin to F actin and vice versa, or as reactions in which the hydrolysis state of the subunit changes. The evolution of the concentrations of a given tip subunit hydrolysis states in principle depends on the hydrolysis state of the subunit adjacent to the tip, which itself depends on the hydrolysis state of the next subunit in the filament, and so on. Every subunit in the filament then requires explicit tracking, causing the dimensionality of the model to be roughly equal to the degree of polymerization of a filament, which typically contains hundreds of subunits. The major accomplishment of the BC model is to truncate this set of recursion equations by assuming that the hydrolysis state of the subunit adjacent to the tip depends only on the hydrolysis state of the tip subunit. Using the results of stochastic simulations of a more complete model, they write empirical equations to capture these relationships, and in doing so they close off an 11-dimensional subset of the original hundreds of equations. They find close agreement between their truncated model and the full stochastic simulation over a wide and realistic range of parameters.

The equations of motion for the 11 variables in the BC model can be written as a nonlinear system of ODEs:

$$\begin{eqnarray}&&\dot{{\bf{x}}}={\bf{f}}({\bf{x}}),\end{eqnarray} \tag{ 1 }$$

where ${\bf{x}}$ is a vector of the concentrations of the 11 species, and ${\bf{f}}$ is a nonlinear vector-valued function of ${\bf{x}}$. This system of equations must be solved numerically, but the results match well with simulations that accurately model experimental data [21]. The steady-state vector ${{\bf{x}}}^{{\rm{ss}}}$ satisfying ${\bf{f}}({{\bf{x}}}^{{\rm{ss}}})={\bf{0}}$ is unique for given values of N and M and under the condition that all concentrations be real and non-negative, and it is attracting. We give more details of the BC model in appendix A, where we list the 11 ODEs.

A separation of timescales exists between the dynamics of the non-tip subunit states and those of tip subunit states: the latter evolve much more rapidly than the former. Thus we partition the vector ${\bf{x}}$ into slow and fast variables: ${{\bf{x}}}_{{\rm{s}}}\equiv {({G}^{T},{G}^{D},{F}^{T},{F}^{{\rm{Pi}}},{F}^{D})}^{\top }$, ${{\bf{x}}}_{{\rm{f}}}\equiv {({T}_{+}^{T},{T}_{+}^{{\rm{Pi}}},{T}_{+}^{D},{T}_{-}^{T},{T}_{-}^{{\rm{Pi}}},{T}_{-}^{D})}^{\top }$, where the subscripts 's' and 'f' refer to 'slow' and 'fast'. Terms of comparable magnitude appear in the equations of motion for both ${{\bf{x}}}_{{\rm{f}}}$ and ${{\bf{x}}}_{{\rm{s}}}$, but the elements in ${{\bf{x}}}_{{\rm{f}}}$ are typically much smaller than those in ${{\bf{x}}}_{{\rm{s}}}$ because $N\ll M$, so ${{\bf{x}}}_{{\rm{f}}}$ reaches equilibrium more rapidly since the distance to equilibrium is not as large as it is for ${{\bf{x}}}_{{\rm{s}}}$. To see the separation timescales concretely, we refer the reader to the end of section 3.1, where we show the spectrum of the Jacobian matrix of the BC model evaluated at steady-state. With the new collective variables ${{\bf{x}}}_{{\rm{s}}}$ and ${{\bf{x}}}_{{\rm{f}}}$, equation (1) can be usefully rewritten as follows:

$$\begin{eqnarray}&&{\dot{{\bf{x}}}}_{{\rm{s}}}={{\bf{Ax}}}_{{\rm{s}}}+{{\bf{Bx}}}_{{\rm{f}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\dot{{\bf{x}}}}_{{\rm{f}}}={\bf{h}}({{\bf{x}}}_{{\rm{s}}},{{\bf{x}}}_{{\rm{f}}}),\end{eqnarray} \tag{ 3 }$$

where ${\bf{A}}$ and ${\bf{B}}$ are matrices whose off-diagonal elements are combinations of kinetic rate constants and whose columns sum to zero due to conservation of actin, and ${\bf{h}}$ is a nonlinear vector-valued function containing terms that are up to cubic products of variables.

The BC model provides a computationally accessible model of the dynamics of actin polymerization, however we might ask for several other features in such a model. We ask that it: (1) be low-dimensional, (2) capture the interesting and important timescales, and (3) be exactly solvable. To meet these goals, we make the decision to track only the vector ${{\bf{x}}}_{{\rm{s}}}$ explicitly. If we were to track the concentration of the tip subunit states, i.e. the elements in ${{\bf{x}}}_{{\rm{f}}}$, we would automatically increase the dimensionality of the model, and we will discuss reasons why we may assume that the tip subunits are evolving in such a way that keeping track of them explicitly is unnecessary. We make two approximations for how to treat ${{\bf{x}}}_{{\rm{f}}}$ in equation (2), first utilizing the fact that a separation of timescale exists between the dynamics of ${{\bf{x}}}_{{\rm{f}}}$ and ${{\bf{x}}}_{{\rm{s}}}$, and then utilizing the fact that $| {{\bf{Bx}}}_{{\rm{f}}}| \ll | {{\bf{Ax}}}_{{\rm{s}}}|$ except when the system is near its steady-state, at which point these terms have comparable magnitudes. This second approximation follows since ${{\bf{x}}}_{{\rm{f}}}$ contains terms up to ${ \mathcal O }(N)$ and ${{\bf{x}}}_{{\rm{s}}}$ contains terms up to ${ \mathcal O }(M)$.

It is also possible to demonstrate that a low-dimensional description of the slow dynamics is a valid approximation through the use of diffusion mapping on a stochastically generated data set based on the BC model. We describe this analysis in supporting information 1 which is available online at stacks.iop.org/NJP/19/125012/mmedia.

The QSSA relies on a separation of timescales between fast and slow variables. This separation allows one to assume that the fast variables ${{\bf{x}}}_{{\rm{f}}}$ are always in equilibrium with respect to the slow variables ${{\bf{x}}}_{{\rm{s}}}$, and therefore that the values of the slow variables determine the values of the fast variables at any moment. We can imagine that ${{\bf{x}}}_{{\rm{f}}}$ is effectively being 'dragged around' by the values of the elements in ${{\bf{x}}}_{{\rm{s}}}$. So, we can solve for functions ${{\bf{x}}}_{{\rm{f}}}^{{\rm{eq}}}({{\bf{x}}}_{{\rm{s}}})$ relating the quasi-equilibrated fast variables in terms of the slow variables by imagining holding ${{\bf{x}}}_{{\rm{s}}}$ fixed and finding the equilibrium values of ${{\bf{x}}}_{{\rm{f}}}$. This amounts to the condition ${\bf{h}}({{\bf{x}}}_{{\rm{f}}}^{{\rm{eq}}};{{\bf{x}}}_{{\rm{s}}})={\bf{0}}$. The functions ${{\bf{x}}}_{{\rm{f}}}^{{\rm{eq}}}({{\bf{x}}}_{{\rm{s}}})$ are then substituted in the equations of motion for the slow variables giving the closed system of equations

$$\begin{eqnarray}&&{\dot{{\bf{x}}}}_{{\rm{s}}}={{\bf{Ax}}}_{{\rm{s}}}+{{\bf{Bx}}}_{{\rm{f}}}^{{\rm{eq}}}({{\bf{x}}}_{{\rm{s}}}).\end{eqnarray} \tag{ 4 }$$

This subsystem is lower-dimensional, though it is nonlinear since ${{\bf{x}}}_{{\rm{f}}}^{{\rm{eq}}}({{\bf{x}}}_{{\rm{s}}})$ is nonlinear, and it describes the evolution of the system on the slow timescales.

In the BC model, the condition ${\bf{h}}({{\bf{x}}}_{{\rm{f}}}^{{\rm{eq}}};{{\bf{x}}}_{{\rm{s}}})={\bf{0}}$ implies the following algebraic systems of equations (see appendix A):

$$\begin{eqnarray} & & \displaystyle \frac{{{\rm{d}}{T}}_{\pm }^{T}}{{\rm{d}}{t}}=0,\\ & & \displaystyle \frac{{{\rm{d}}{T}}_{\pm }^{{\rm{Pi}}}}{{\rm{d}}{t}}=0,\\ & & \displaystyle \frac{{{\rm{d}}{T}}_{\pm }^{D}}{{\rm{d}}{t}}=0.\end{eqnarray} \tag{ 5 }$$

Only four of these six equations are linearly independent due to the conservation of number of plus and minus end filament tips, so we use the following supplementary equations to find a solution of the combined systems of equations:

$$\begin{eqnarray}&&{T}_{\pm }^{T}+{T}_{\pm }^{{\rm{Pi}}}+{T}_{\pm }^{D}=N.\end{eqnarray} \tag{ 6 }$$

System of equations (5), (6) can be solved numerically resulting in tabulated functions of the forms ${T}_{\pm }^{T}({G}^{T},{G}^{D})$, ${T}_{\pm }^{{\rm{Pi}}}({G}^{T},{G}^{D})$, and ${T}_{\pm }^{D}({G}^{T},{G}^{D})$. These functions do not depend on F^T, ${F}^{{\rm{Pi}}}$, and F^D because these variables do not enter into the the function ${\bf{h}}$. Subsequently, the nonlinear (slow) system described by equation (4) is numerically integrated. Figure 2 displays a comparison of the QSSA approximation to the original 11-dimensional BC model.

Zoom In Zoom Out Reset image size

This approximation succeeds in reducing the dimensionality of the model to 5. We note, however, that the dynamics of this model lie on a four-dimensional submanifold of the full five-dimensional manifold due to the fact that ${\bf{A}}$ and ${\bf{B}}$ are both singular, corresponding to conservation of actin. This model also captures the interesting timescales describing polymerization events and dynamics of the hydrolysis states of F and G actin, while it neglects the fast dynamics corresponding to events taking place at the tips. However the model is still nonlinear and analytically unsolvable, so we propose an alternative model reduction scheme.

In addition to having sufficiently fast dynamics as to be effectively described as in quasi-equilibrium with respect to ${{\bf{x}}}_{{\rm{s}}}$, the vector ${{\bf{x}}}_{{\rm{f}}}$ is also small in magnitude compared to ${{\bf{x}}}_{{\rm{s}}}$, and this fact can be utilized to drastically simplify the QSSA model. Although the tip dynamics are fast, one can profitably assume that the concentration of filament tip subunit states, that is, the elements in ${{\bf{x}}}_{{\rm{f}}}$, are constant in time. We refer to this assumption as the constant tip (CT) approximation. Certainly this assumption is not realistic since the tips have fast dynamics, but the effect of this error on the equations of motion of the other actin species is small in the regime where the number concentration of filaments is much smaller than the total amount of actin in the system, i.e. when $N\ll M$, or if the tip subunit states quickly attain their limiting values. This assumption allows the model to be reduced to a five-dimensional linear system of equations that can be solved analytically. The procedure is to replace the dynamical variables ${T}_{j}^{i}$ by constants $N{{\rm{\Gamma }}}_{j}^{i}$, chosen such that the steady-state of the CT model coincides with the steady-state of the BC model. We define ${{\rm{\Gamma }}}_{j}^{i}$ as the fraction of the filament tips at the j (plus or minus) end that are in the i (T, Pi, D) hydrolysis state. These constants determine the rate of depolymerization reactions, and replacing the variables ${T}_{j}^{i}$ with them causes the term ${{\bf{Bx}}}_{{\rm{f}}}$ to be a constant source and sink term in equation (2), which as a result becomes linear.

The steady-states of the BC and CT models will be the same if

$$\begin{eqnarray} & & N{{\rm{\Gamma }}}_{\pm }^{T}=\mathop{{\rm{lim}}}\limits_{t\to \infty }{T}_{\pm }^{T}(t),\\ & & N{{\rm{\Gamma }}}_{\pm }^{{\rm{Pi}}}=\mathop{{\rm{lim}}}\limits_{t\to \infty }{T}_{\pm }^{{\rm{Pi}}}(t),\\ & & N{{\rm{\Gamma }}}_{\pm }^{D}=\mathop{{\rm{lim}}}\limits_{t\to \infty }{T}_{\pm }^{D}(t).\end{eqnarray} \tag{ 7 }$$

In other words, the CT subunit state fractions in the CT model should be chosen as the steady-state values of the tip subunit state fractions in the BC model. As a result, only the approach to steady-state will be different between the two models. These limiting values can be found by numerically solving the algebraic system of equations

$$\begin{eqnarray} & & {\bf{f}}({{\bf{x}}}^{{\rm{ss}}})=0,\\ & & {T}_{\pm }^{T}+{T}_{\pm }^{{\rm{Pi}}}+{T}_{\pm }^{D}=N,\\ & & {G}^{T}+{G}^{D}+{F}^{T}+{F}^{{\rm{Pi}}}+{F}^{D}=M,\end{eqnarray} \tag{ 8 }$$

where ${{\bf{x}}}^{{\rm{ss}}}$ is the 11-dimensional steady-state vector of concentrations in the BC model, and taking the real non-negative solution. Equation (8) determines the values of ${{\rm{\Gamma }}}_{j}^{i}$ which will be unique for a given N and M.

We define

$$\begin{eqnarray}&&{\bf{b}}\equiv {{\bf{Bx}}}_{{\rm{f}}}^{{\rm{ss}}},\end{eqnarray} \tag{ 9 }$$

where ${{\bf{x}}}_{{\rm{f}}}^{{\rm{ss}}}$ contains the constants $N{{\rm{\Gamma }}}_{j}^{i}$ instead of the variables ${T}_{j}^{i}$. The equation of motion for ${{\bf{x}}}_{{\rm{s}}}$ is

$$\begin{eqnarray}&&{\dot{{\bf{x}}}}_{{\rm{s}}}={{\bf{Ax}}}_{{\rm{s}}}+{\bf{b}},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}{\bf{A}}=\left(\begin{array}{ccccc}-a-b & c & 0 & 0 & 0\\ b & -c-d & 0 & 0 & 0\\ a & 0 & -e & f & 0\\ 0 & 0 & e & -f-g & h\\ 0 & d & 0 & g & -h\end{array}\right),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\bf{b}}=\left(\begin{array}{c}i\\ j+k\\ -i\\ -k\\ -j\end{array}\right)\end{eqnarray} \tag{ 12 }$$

with

$$\begin{eqnarray}\begin{array}{ll}a=N({k}_{{\rm{on}},\,+}^{T}+{k}_{{\rm{off}},\,+}^{T}) & g={k}_{{\rm{phos}}}\\ b={k}_{{\rm{nex}}} & h={k}_{{\rm{rephos}}}\\ c={k}_{{\rm{renex}}} & i=N({{\rm{\Gamma }}}_{+}^{T}{k}_{{\rm{off}},\,+}^{T}+{{\rm{\Gamma }}}_{-}^{T}{k}_{{\rm{off}},\,-}^{T})\\ d=N({k}_{{\rm{on}},\,+}^{D}+{k}_{{\rm{off}},\,+}^{D}) & j=N({{\rm{\Gamma }}}_{+}^{D}{k}_{{\rm{off}},\,+}^{D}+{{\rm{\Gamma }}}_{-}^{D}{k}_{{\rm{off}},\,-}^{D})\\ e={k}_{{\rm{hyd}}} & k=N({{\rm{\Gamma }}}_{+}^{{\rm{Pi}}}{k}_{{\rm{off}},\,+}^{{\rm{Pi}}}+{{\rm{\Gamma }}}_{-}^{{\rm{Pi}}}{k}_{{\rm{off}},\,-}^{{\rm{Pi}}}).\\ f={k}_{{\rm{rehyd}}} & \end{array}\end{eqnarray} \tag{ 13 }$$

Equation (10) is a linear system of differential equations which can be solved exactly. One point of difficulty in solving this system is that the columns of ${\bf{A}}$ sum to zero due to conservation of actin, causing ${\bf{A}}$ to be singular. In chemical reaction network theory, one often has such systems with linear conservation laws. If the system is linear, with only first order or pseudo-first order reactions, then such a singular system can be solved cleanly by a method using the Drazin inverse ${{\bf{A}}}^{{ \mathcal D }}$ of the matrix ${\bf{A}}$. We believe that this method has certain advantages over other approaches to solving singular systems of differential equations. If ${\bf{A}}$ has Jordan decomposition

$$\begin{eqnarray}{\bf{A}}={\bf{V}}\left(\begin{array}{cc}{J}_{1} & 0\\ 0 & {J}_{0}\end{array}\right){{\bf{V}}}^{-1},\end{eqnarray} \tag{ 14 }$$

where J₁ and J₀ correspond to the non-zero and zero eigenvalues respectively, then

$$\begin{eqnarray}{{\bf{A}}}^{{ \mathcal D }}={\bf{V}}\left(\begin{array}{cc}{J}_{1}^{-1} & 0\\ 0 & 0\end{array}\right){{\bf{V}}}^{-1}.\end{eqnarray} \tag{ 15 }$$

In supporting information 2 in the supplementary material we go over the details of solving equation (10) as well as the advantages of this method [22]. The solution is

$$\begin{eqnarray}&&{{\bf{x}}}_{{\rm{s}}}=(-{{\bf{A}}}^{{ \mathcal D }}+{{\rm{e}}}^{{\bf{A}}t}{\bf{G}}){\bf{b}}.\end{eqnarray} \tag{ 16 }$$

The matrix

$$\begin{eqnarray}&&{\bf{G}}=\displaystyle \frac{1}{| {\bf{b}}{| }^{2}}{{\bf{x}}}_{{\rm{s}}}(0){{\bf{b}}}^{\top }+{{\bf{A}}}^{{ \mathcal D }}\end{eqnarray} \tag{ 17 }$$

encodes information about the initial conditions ${{\bf{x}}}_{{\rm{s}}}(0)$.

Perhaps the most important benefit of an exactly solvable model is the ability to formally determine the eigenvalues governing the linear dynamics. These eigenvalues characterize the relaxation times of the components of a perturbation from equilibrium, in the basis of the eigenvectors. The expressions for the eigenvalues thus shed light on the timescales that describe the different chemical processes. In our modeling we have included the reverses of kinetically dominant forward reactions, and we have set the rate constants of these reactions to be equal to something on the order of the corresponding forward reaction rate constants multiplied by a small parameter . Thus $b=\epsilon {b}^{* }$, where ${b}^{* }\sim c$, $f=\epsilon {f}^{* }$, where ${f}^{* }\sim e$, and $h=\epsilon {h}^{* }$, where ${h}^{* }\sim g$. By writing reverse rate constants this way, we can Taylor expand the expressions for the eigenvalues around the point $\epsilon =0$ to simplify the result. Doing so to first order in , we have

$$\begin{eqnarray}&&{\lambda }_{1}=0,\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}{\lambda }_{2} & = & \displaystyle \frac{1}{2}(-a-\epsilon {b}^{* }-c-d+\sqrt{{(a-d+\epsilon {b}^{* }-c)}^{2}+4\epsilon {b}^{* }c})\\ & \approx & -c-d+\displaystyle \frac{{b}^{* }c}{a-c-d}\epsilon .\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}{\lambda }_{3} & = & \displaystyle \frac{1}{2}(-a-\epsilon {b}^{* }-c-d-\sqrt{{(a-d+\epsilon {b}^{* }-c)}^{2}+4\epsilon {b}^{* }c})\\ & \approx & -a-\displaystyle \frac{{b}^{* }(a-d)}{a-c-d}\epsilon ,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}{\lambda }_{4} & = & \displaystyle \frac{1}{2}(-e-\epsilon {f}^{* }-g-\epsilon {h}^{* }+\sqrt{{(e-g+\epsilon {f}^{* }-\epsilon {h}^{* })}^{2}+4\epsilon {f}^{* }g})\\ & \approx & -g-\left({h}^{* }-\displaystyle \frac{{{gf}}^{* }}{e-g}\right)\epsilon ,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}{\lambda }_{5} & = & \displaystyle \frac{1}{2}(-e-\epsilon {f}^{* }-g-\epsilon {h}^{* }-\sqrt{{(e-g+\epsilon {f}^{* }-\epsilon {h}^{* })}^{2}+4\epsilon {f}^{* }g})\\ & \approx & -e-\displaystyle \frac{{{ef}}^{* }}{e-g}\epsilon .\end{eqnarray} \tag{ 22 }$$

The fact that the nonzero eigenvalues are negative implies the stability of the steady-state. These eigenvalues have no dependence on M, the total concentration of actin subunits, but they do depend on N, the number concentration of filaments, since this term appears in the expressions for a and d. For the parameterization used here and with $N=1\ {\rm{nM}}$, the eigenvalues can be ordered by magnitude as follows:

$$\begin{eqnarray}&&| {\lambda }_{5}| \gt | {\lambda }_{3}| \gt | {\lambda }_{2}| \gt | {\lambda }_{4}| \gt | {\lambda }_{1}| .\end{eqnarray} \tag{ 23 }$$

Equations (18)–(22) are in terms of reaction rate constants and can be interpreted as representing certain collective subprocesses in the chemical system corresponding to combinations of those reactions. The ordering indicates the comparative rates of those subprocesses. We interpret the zeroth order terms of the eigenvalues as follows:

• ${\lambda }_{1}=0$ because actin subunits are conserved in this system, causing ${\bf{A}}$ to be singular. Equivalently, one can say that the dynamics unfold on a four-dimensional submanifold of the five-dimensional manifold, and that this submanifold is determined by M.
• ${\lambda }_{2}$ represents the combination of the forward nucleotide exchange reaction (${G}^{D}\to {G}^{T}$) and the polymerization of G^D. Both of these reactions convert G^D into another species, so this eigenvalue represents the subprocess of G^D depletion.
• ${\lambda }_{3}$ represents the polymerization of G^T.
• ${\lambda }_{4}$ represents the release of phosphate by ${F}^{{\rm{Pi}}}$ to form F^D.
• ${\lambda }_{5}$ represents the hydrolysis of ATP converting F^T to ${F}^{{\rm{Pi}}}$.

Whether the magnitudes of these eigenvalues are increased or decreased due to the presence of reversible reactions (i.e. when $\epsilon \gt 0$) depends on the parameterization, since the sign of the first order terms depend on the comparative sizes of certain parameters.

In the full BC 11-dimensional model, one can numerically evaluate the Jacobian matrix of ${\bf{f}}({\bf{x}})$ at steady-state,

$$\begin{eqnarray}&&{{\bf{J}}}^{* }\equiv {\left.\displaystyle \frac{\partial {\bf{f}}}{\partial {\bf{x}}}\right|}_{{\bf{x}}={{\bf{x}}}^{{\rm{ss}}}}.\end{eqnarray} \tag{ 24 }$$

We find that, for the same parameterization, the smallest four non-zero eigenvalues of ${{\bf{J}}}^{* }$ are numerically equal to the non-zero eigenvalues of ${\bf{A}}$. This implies that the CT model has captured the slowest dynamics of the BC model by ignoring the processes involving the tip subunits. The main benefit of the CT approximation is that these slow linear dynamics can now be easily analyzed. These dynamics provide information about the polymerization process, the nucleotide composition of the filaments, and nucleotide exchange of globular actin. For most purposes, these aspects are of primary interest, and the processes at the tips are of lesser importance.

To quantitatively judge the separation of timescales of the BC model, we list the nonzero numerically evaluated eigenvalues of ${{\bf{J}}}^{* }$ when $N=1\ {\rm{nM}}$ and $\epsilon =0$, in decreasing magnitude $(-5.49,-1.41,-0.833,-0.393,-0.300,-0.0130,-0.0129,-0.002)$. For comparison, we do the same with the eigenvalues of ${\bf{A}}$ $(-0.300,-0.0130,-0.0129,-0.002)$. While there is no large spectral gap, there is certainly a wide range of timescales, and it would suffice to even keep the smallest 3 non-zero eigenvalues. We later show how, by combining certain species, the dimensionality can be reduced to 3.

One might expect that if we increase the amount of actin subunits in the system by a different choice of initial conditions, the concentrations of the different species at steady-state would change. An interesting feature of this system is that this is true only for some species, and which species it is true for depends on whether or not we have included reversible reactions (if $\epsilon \gt 0$). Additionally, this can be shown to be true in both the CT model and the BC model. We demonstrate it first in the CT model.

We find the steady-state vector of concentrations ${{\bf{x}}}_{{\rm{s}}}^{{\rm{ss}}}$ by taking the limit of equation (16) as $t\to \infty$,

$$\begin{eqnarray}&&{{\bf{x}}}_{{\rm{s}}}^{{\rm{ss}}}\equiv \mathop{\mathrm{lim}}\limits_{t\to \infty }{{\bf{x}}}_{{\rm{s}}}(t)=-{{\bf{A}}}^{{ \mathcal D }}{\bf{b}}+{\bf{CGb}},\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&{\bf{C}}\equiv \mathop{\mathrm{lim}}\limits_{t\to \infty }{{\rm{e}}}^{{\bf{A}}t}.\end{eqnarray} \tag{ 26 }$$

We eigendecompose ${\bf{A}}$ as ${{\bf{UDV}}}^{\dagger }$ and use it in the expression for ${\bf{C}}$,

$$\begin{eqnarray}&&{\bf{C}}={\bf{U}}\mathop{\mathrm{lim}}\limits_{t\to \infty }{{\rm{e}}}^{{\bf{D}}t}{{\bf{V}}}^{\dagger }.\end{eqnarray} \tag{ 27 }$$

With the exception of ${\lambda }_{1}$, which is zero, the eigenvalues of ${\bf{A}}$ are negative, so we have

$$\begin{eqnarray}{\bf{C}} & = & {\bf{U}}\ {\rm{diag}}(1,0,0,0,0)\ {{\bf{V}}}^{\dagger }\\ & = & \left(\begin{array}{ccccc}0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ \displaystyle \frac{{fh}}{{eg}} & 0 & 0 & 0 & 0\\ \displaystyle \frac{h}{g} & 0 & 0 & 0 & 0\\ 1 & 0 & 0 & 0 & 0\end{array}\right){{\bf{V}}}^{\dagger },\end{eqnarray} \tag{ 28 }$$

where ${\left(0,0,\tfrac{{fh}}{{eg}},\tfrac{h}{g},1\right)}^{\top }$ is the eigenvector corresponding to ${\lambda }_{1}$. Thus the top two rows of ${\bf{C}}$ are zero, and the top four rows would be zero if we exclude reversible reactions. This is also true of the term ${\bf{CGb}}$. Now, the matrix ${\bf{G}}$ is the only place where the initial conditions appear in equation (25). So if the top two rows ${\bf{CGb}}$ are zero, then the first two elements of ${{\bf{x}}}_{{\rm{s}}}^{{\rm{ss}}}$ cannot have any dependence on initial conditions. In other words, G^D and G^T always reach the same concentrations at steady-state no matter what the initial concentrations of any of the species are. If we have no reversible reactions, the same is also true for the species F^T and ${F}^{{\rm{Pi}}}$.

Consider the following thought experiment, assuming for simplicity $\epsilon =0$. We start with some initial amount M of actin in any form and let the system come to steady-state. We then add an amount ${\rm{\Delta }}M$ more actin to the system, in any form, and wait until the system has reached steady-state again. We would find that the only difference between the two steady-states is that the concentration of F^D had increased by ${\rm{\Delta }}M$. If $\epsilon \gt 0$, then we would find that the steady-state values of F^T, and ${F}^{{\rm{Pi}}}$, F^D had all increased, and the sum of these changes would be ${\rm{\Delta }}M$.

This lack of dependence of the steady-state concentrations of some species on M is not an artifact of the CT approximation; it is also the case in the BC model. It can not be shown to be true by taking the limit $t\to \infty$ as is the case here, but it can instead be shown by observing that a subset of the system of algebraic equations ${\bf{f}}({{\bf{x}}}^{{\rm{ss}}})=0$ is closed, and that a solution for the subset could be obtained without specifying M. This implies that the steady-state concentrations of the species represented by that solution have no dependence on M. We give the details of this argument in appendix B.