Cell Size, Mechanical Tension, and GTPase Signaling in the Single Cell

Abstract

Cell polarization requires redistribution of specific proteins to the nascent front and back of a eukaryotic cell. Among these proteins are Rac and Rho, members of the small GTPase family that regulate the actin cytoskeleton. Rac promotes actin assembly and protrusion of the front edge, whereas Rho activates myosin-driven contraction at the back. Mathematical models of cell polarization at many levels of detail have appeared. One of the simplest based on “wave-pinning” consists of a pair of reaction–diffusion equations for a single GTPase. Mathematical analysis of wave-pinning so far is largely restricted to static domains in one spatial dimension. Here, we extend the analysis to cells that change in size, showing that both shrinking and growing cells can lose polarity. We further consider the feedback between mechanical tension, GTPase activation, and cell deformation in both static, growing, shrinking, and moving cells. Special cases (spatially uniform cell chemistry, the absence or presence of mechanical feedback) are analyzed, and the full model is explored by simulations in 1D. We find a variety of novel behaviors, including “dilution-induced” oscillations of Rac activity and cell size, as well as gain or loss of polarization and motility in the model cell.

Appendices

A Mass Balance on Growing or Shrinking Domain

Consider a prototypical reaction diffusion equation for some concentration u(x, t) in conservative form,

$$\begin{aligned} u_t = - \nabla \cdot {\mathbf {j}}(u) + f(u), \quad \text{ in }\ \Gamma (t), \end{aligned}$$

where ${\mathbf {j}}(u)$ is the particle flux. Further, assume no-flux boundary conditions on the boundary of the domain, $\partial \Gamma (t)$. In integral form, we write a conservation law for a small volume element $V(t) = [x(t), x(t) + \Delta x(t)]$

$$\begin{aligned} \frac{\, \mathrm {d}}{\, \mathrm {d}t} \int _{V(t)} u(x, t) \, \mathrm {d}x = \int _{V(t)} \left[ \, -\nabla \cdot {\mathbf {j}} + f(u) \,\right] \, \mathrm {d}x. \end{aligned}$$

Using the Reynolds transport theorem, the left hand side is equivalent to

$$\begin{aligned} \frac{\, \mathrm {d}}{\, \mathrm {d}t} \int _{V(t)} u(x, t) \, \mathrm {d}x = \int _{V(t)} \left[ \, u_t + \nabla \cdot \left( {\mathbf {v}} u\right) \,\right] \, \mathrm {d}x, \end{aligned}$$

where ${\mathbf {v}}$ is the velocity of the material point x(t). We then obtain the following reaction–transport–diffusion equation

$$\begin{aligned} u_t + \nabla \cdot \left( {\mathbf {v}} u \right) = D \Delta u + f(u), \qquad \text{ in }\ [x_1(t), x_2(t)]. \end{aligned}$$

(14)

Similarly, the boundary conditions on $\Gamma (t)$ become

$$\begin{aligned} -D \nabla u + {\mathbf {v}} u = 0, \quad \text{ on }\ \partial \Gamma (t). \end{aligned}$$

We have two new terms

• 1. An advection term, ${\mathbf {v}} \cdot \nabla u$, which corresponds to elemental volumes moving with the flow due to local growth.

• 2. A dilution term, $u \nabla \cdot {\mathbf {v}}$, due to local volume changes.

B Numerical Implementation

1.1 B.1 Mapping to a Constant Domain

It is most convenient to discretize and numerically solve the PDEs on a constant domain. Hence, a mapping is used to transform the growing/shrinking cell to the interval $[-1,1]$. This can be achieved by the following change of variables

$$\begin{aligned} (x,t) \mapsto ({{\bar{x}}}, {{\bar{t}}}) = \left( \frac{x - x_c(t)}{R(t)}, t \right) . \end{aligned}$$

Note that ${\bar{x}}$ is the location in the new domain $[-1,1]$, whereas x is located in the cell domain $[x_1(t), x_2(t)]$. Carefully applying the chain rule leads to

$$\begin{aligned} \begin{aligned} u_{{\bar{x}}}&= \frac{\partial u}{\partial x} \frac{\partial x}{\partial {{{\bar{x}}}}}= u_x R(t) \\ u_{{\bar{t}}}&= \frac{\partial u}{\partial t} \frac{\partial t}{\partial {{{\bar{t}}}}}+ \frac{\partial u}{\partial x} \frac{\partial x}{\partial t}\frac{\partial t}{\partial {{{\bar{t}}}}} = u_t + u_x \left( {\bar{x}} R'(t) + x_c'(t) \right) , \end{aligned} \end{aligned}$$

where we have used $x={{{\bar{x}}}}R(t)+x_c(t)$ and $t={{{\bar{t}}}}$ in the transformation. Substituting these results into Eq. (2), we obtain

$$\begin{aligned}&u_{{{\bar{t}}}} - u_{{{\bar{x}}}} \frac{{\bar{x}} R'(t) + x_c'(t)}{R(t)} \\&\quad + \frac{x_c'(t) + {\bar{x}} R'(t)}{R(t)} u_{{{\bar{x}}}} + u \frac{R'(t)}{R(t)} = \frac{D}{R^2(t)} { u_{{{\bar{x}}}{{\bar{x}}}} }+ f(u), \ \text{ in }\ [-1, 1]. \end{aligned}$$

Canceling two terms on the LHS and rearranging leads to

$$\begin{aligned} u_{{{\bar{t}}}} = \frac{D}{R^2(t)}u_{{\bar{x}}{\bar{x}}} + f(u) - u \left( \frac{\dot{R}(t)}{R(t)}\right) ,\ \text{ in }\ [-1,1]. \end{aligned}$$

The total amount of GTPase on the interval $[-1,1]$ becomes time dependent through

$$\begin{aligned} G_\mathrm{{T}} = \int _{-R(t)}^{R(t)} \left( G + G_i \right) \, \mathrm {d}x = R(t) \int _{-1}^{1} \left( G + G_i \right) \, \mathrm {d}y. \end{aligned}$$

1.2 B.2 Changing Rest Length is Equivalent to Imposing Additional Force

It is easy to see from Eq. (10) that adjusting the cell’s rest length, $L_0$ is equivalent to assuming that additional forces act on the cell:

$$\begin{aligned} (\mu +2\eta ) \frac{\, \mathrm {d}L}{\, \mathrm {d}t} = (F_2-F_1 ) -2 E(L - L_0) = -2 E(L - (L_0 + \Delta L_0)), \end{aligned}$$

(15)

where $\Delta L_0 = (F_2 - F_1) / 2E$. In other words, a change of $\Delta L_0$ is equivalent to an additional force of $-2E\Delta L_0$.

In Zmurchok et al. (2018) and Liu (2019) only the spatially uniform G model variants were considered. There, cell length was modeled by (15), with

$$\begin{aligned} \Delta L_0 (G)= \pm b \frac{G^m}{G_c^m + G^m}. \end{aligned}$$

Rac (+) is assumed to increase $L_0$ and Rho (-) to decrease it, so the sign of the Hill function depends on the given GTPase.

As shown here, the assumption that GTPases work through modification of the cell’s effective rest length is mathematically equivalent to assumptions we make about GTPase-dependent load forces on the cell ends.

1.3 B.3 Numerical Methods

Each cell introduces a set of mechanical equations for its endpoints, i.e., Eq. (9). Next, each cell contains K GTPases, that is, we include K sets of reaction–diffusion Eq. (1). The cell interior is divided into a cell-centered grid with uniform length $h = 1 / N$, where N is the number of grid cells per unit length. Here, we use $N = 2^{10}$. The second-order derivative is discretized using a second-order cell-centered finite difference. For more detailed information on the numerical method, see Gerisch (2001). This results in M total grid points per GTPase species and K GTPase species. The resulting state vector is given by

$$\begin{aligned} y = \begin{pmatrix} x_1 \\ x_2 \\ G_{1, 1} \\ G_{1, 2} \\ \vdots \\ G_{K, M} \end{pmatrix} \end{aligned}$$

where $G_{j, k}$ denotes the concentration of the k-th GTPase at the j-th spatial grid point. The temporal evolution of the state of the cell can be written as:

$$\begin{aligned} \frac{\, \mathrm {d}y}{\, \mathrm {d}t} = f(t, y), \end{aligned}$$

where f(t, y) is a function whose first two entries contain Eq.(9) and the final KM entries are the result of the spatial discretization of the GTPase PDEs.

The resulting system of ordinary differential equations are integrated using the ROWMAP integrator (Weiner et al. 1997), an specialized integrator for large stiff ODEs (such systems are commonly the result of spatially discretizing PDEs). Here, we use the ROWMAP implementation provided by the authors of Weiner et al. (1997).¹ The ROWMAP integrator (written in Fortran) is wrapped for easy use in a Python environment using f2py² (f2py provides a connection between the Python and Fortran languages) and integrated into scipy’s³ integrator class (a package providing several different techniques to integrate ODEs). The spatial discretization (right hand side of ODE) is implemented using NumPy.⁴ The error tolerance of the integrator is set to $v_{tol} = 10^{-6}$.

C Linear Stability of Spatially Uniform GTPase Model

In this section, we consider the steady states and their linear stability of Eq. (13).

The Jacobian of the ODE system (13) at any steady state $(L_\mathrm{{eq}}, G_{eq})$ is

$$\begin{aligned}&J(L_\mathrm{{eq}}, G_{eq}) \\&\quad = \begin{pmatrix} -\frac{2E}{\mu + 2 \eta } &{} \pm \frac{2F_{\text {max}}}{\mu +2\eta } \frac{mG_c^m G_\mathrm{{eq}}^{m-1}}{\left( G_c^m + G_\mathrm{{eq}}^m\right) ^2} \\ -\frac{G_\mathrm{{T}}}{L_\mathrm{{eq}}^2}\left( b + \gamma \frac{G_\mathrm{{eq}}^n}{1 + G_\mathrm{{eq}}^n}\right) &{} \frac{\gamma n G_\mathrm{{eq}}^{n-1}}{\left( G_\mathrm{{eq}}^n + 1 \right) ^2}\left( \frac{G_\mathrm{{T}}}{L_\mathrm{{eq}}} - G_\mathrm{{eq}}\right) - \left( b + I + \gamma \frac{G_\mathrm{{eq}}^n}{1 + G_\mathrm{{eq}}^n}\right) \end{pmatrix}. \end{aligned}$$

Model parameters are all positive, and $L_\mathrm{{eq}}, G_{eq} > 0$ so in the Rac and Rho cases, respectively, the Jacobian has the sign patterns:

$$\begin{aligned} J_{\text {Rac}} = \begin{pmatrix} - &{} + \\ - &{} {\text {sgn}}(d) \end{pmatrix}, \qquad J_{\text {Rho}} = \begin{pmatrix} - &{} - \\ - &{} {\text {sgn}}(d) \end{pmatrix}. \end{aligned}$$

where

$$\begin{aligned} d :=\frac{\gamma n G_\mathrm{{eq}}^{n-1}}{\left( G_\mathrm{{eq}}^n + 1 \right) ^2}\left( \frac{G_\mathrm{{T}}}{L_\mathrm{{eq}}} - G_\mathrm{{eq}}\right) - \left( b + I + \gamma \frac{G_\mathrm{{eq}}^n}{1 + G_\mathrm{{eq}}^n}\right) . \end{aligned}$$

Note that since $(L, G) \in \Delta$, where $\Delta$ is the invariant region in Result 1, we have that d is the difference of two positive numbers. It is straightforward to conclude the following (see, for instance, Murray (1989)). In the Rho case, the steady states are either stable nodes or saddles and no oscillations are possible. In the Rac case, the steady states are stable when ${\text {sgn}}(d) < 0$, and oscillations are possible when ${\text {sgn}}(d) > 0$.

1.1 C.1 Sharp-Switch Limit

Using a sharp-switch approximation, we can easily compute two steady states:

• 1. When both nonlinear switches are off, i.e., $G \ll 1$, the steady states are:

  $$\begin{aligned} L_\mathrm{{eq}} = L_0, \qquad G_{1} = \frac{b G_\mathrm{{T}}}{L_0 (b + I)}. \end{aligned}$$

  Since the linearization of (13) is lower triangular, the eigenvalues are readily found

  $$\begin{aligned} \lambda _1 = \frac{-2E}{\mu + 2\eta }, \quad \lambda _2 = -(b + I). \end{aligned}$$

  Since $b, I, E, \mu , \eta >0$, we conclude that this steady state is a stable node. Finally, since $G_\mathrm{{eq}} \ll 1$, we obtain an existence condition

  $$\begin{aligned} G_\mathrm{{T}} < L_0 \left( 1 + \frac{I}{b} \right) . \end{aligned}$$

  Therefore, the steady state of a large cell, with low total active GTPase, can only exist if the total amount of GTPase is small.

• 2. In the case $G_c< G < 1$, i.e., one of the nonlinear switches is turned on, the steady states are:

  $$\begin{aligned} L_\mathrm{{eq}} = L_0 \pm \frac{F_{\text {max}}}{E}, \qquad G_{2} = \frac{b G_\mathrm{{T}}}{L_\mathrm{{eq}}(b + I)}. \end{aligned}$$

  Since the linearization of (13) is lower triangular as before, the eigenvalues are readily found

  $$\begin{aligned} \lambda _1 = \frac{-2E}{\mu + 2\eta }, \quad \lambda _2 = -(b + \gamma + I). \end{aligned}$$

  Since $b, I, \gamma , E, \mu , \eta >0$, we conclude that this steady state is a stable node. A condition for the existence of this steady state is

  $$\begin{aligned} \frac{L_0 \left( b + \gamma + I\right) }{b + \gamma }< G_\mathrm{{T}} < \frac{G_c L_0 \left( b + \gamma + I \right) }{b + \gamma }. \end{aligned}$$

  Thus, we must have some intermediate amount of GTPase.

• 3. In the case $1< G < G_c$, i.e., one of the nonlinear switches is turned on, the steady states are:

  $$\begin{aligned} L_\mathrm{{eq}} = L_0, \qquad G_{2} = \frac{(b + \gamma ) G_\mathrm{{T}}}{L_0(b + \gamma + I)} \end{aligned}$$

  Eigenvalues are readily found, as before,

  $$\begin{aligned} \lambda _1 = \frac{-2E}{\mu + 2\eta }, \quad \lambda _2 = -(b + \gamma + I). \end{aligned}$$

  From $b, I, \gamma , E, \mu , \eta >0$ we conclude that this steady state is a stable node. The existence condition for this steady state can be written in either of two forms

  $$\begin{aligned} \frac{L_0 \left( b + \gamma + I\right) }{b + \gamma }< & {} G_\mathrm{{T}}< \frac{G_c L_0 \left( b + \gamma + I \right) }{b + \gamma }, \quad \frac{L_0}{G_\mathrm{{T}}-L_0}I-b< \gamma \\< & {} \frac{L_0 G_c}{G_\mathrm{{T}}-L_0 G_c}I-b. \end{aligned}$$

  For these ranges to be non-empty, we require that $G_\mathrm{{T}} > L_0 G_c$. Thus, we must have some intermediate amount of GTPase.

• 4. In the case, when both nonlinear switches are turned on, i.e., $G \gg G_c$, the steady states are

  $$\begin{aligned} L_\mathrm{{eq}} = L_0 \pm \frac{F_M}{E}, \qquad G_{3} = \frac{(b + \gamma ) G_\mathrm{{T}}}{L_\mathrm{{eq}}(b + \gamma + I)}. \end{aligned}$$

  Then the eigenvalues are

  $$\begin{aligned} \lambda _1 = \frac{-2E}{\mu + 2\eta }, \quad \lambda _2 = -(b + \gamma + I). \end{aligned}$$

  Since $b, I, \gamma , E, \mu , \eta >0$, this steady state is a stable node. The condition for existence of this steady state is given by

  $$\begin{aligned} G_\mathrm{{T}}> G_c L_\mathrm{{eq}} \left( 1 + \frac{I}{b + \gamma }\right) , \ \text{ or }\ \ \gamma > \frac{L_b G_c}{G_\mathrm{{T}}-L_b G_c}I-b. \end{aligned}$$

  Thus, only cells with a sufficiently amount of GTPase can be large or small.

In the bifurcation diagram of Fig. 7, the regions of existence of $G_1, G_2$ or $G_3$ are identified.

D The Case of a Deforming Cell with Conserved Volume

Here, we consider the case that the volume of a cell is roughly constant, while the cell deforms and its membrane size changes. This idea suggests a slight modification to our model. The modified mass conservation reads

$$\begin{aligned} G_\mathrm{{T}} = G_i V + \beta G L, \end{aligned}$$

where V is the constant volume of the cell and $\beta$ is a constant of proportionality with units of [length]$^2$. We can solve for $G_i$ to get a relation analogous to (6):

$$\begin{aligned} G_i = \frac{G_\mathrm{{T}}}{V}-\frac{\beta GL}{V}. \end{aligned}$$

Notice that after substituting this into the well-mixed system, we can set $V=1$ by rescaling b and $\gamma$. Therefore, we get

$$\begin{aligned} \frac{\, \mathrm {d}G}{\, \mathrm {d}t} = \left( b + \gamma \frac{G^n}{1+G^n}\right) \left( G_\mathrm{{T}} - \beta G L\right) - IG - G\left( \frac{{\dot{L}}}{L}\right) . \end{aligned}$$

(16)

This, together with (13a) and (13b), forms our new system.

Much of the analysis of this system is similar to the earlier case. The changes to the nullclines in the $G-L$ phase plane correspond to a vertical shift in the sharp-switch limit, which can be compensated by adjusting the parameters. This suggests that the behaviors of the new system are similar to the original (13). The equilibrium branches and the conditions for their existence in the sharp-switch limit are:

• 1. For $G \ll 1$:

  $$\begin{aligned}L_\mathrm{{eq}}=L_0, \quad G_\mathrm{{eq}}=\frac{b G_\mathrm{{T}}}{b \beta L_0+ I} \quad \text{ exists } \text{ provided } \quad G_\mathrm{{T}} < \beta L_0 + \frac{I}{b}. \end{aligned}$$

• 2. For $G_c<G<1$:

  $$\begin{aligned}L_\mathrm{{eq}}= & {} L_0\pm \frac{F\mathrm{{max}}}{E}, \\ \quad G_\mathrm{{eq}}= & {} \frac{b G_\mathrm{{T}}}{b \beta L_\mathrm{{eq}}+ I} \quad \text{ exists } \text{ provided } \quad G_c \left( \beta L_\mathrm{{eq}}+\frac{I}{b} \right)<G_\mathrm{{T}} <\beta L_\mathrm{{eq}}+\frac{I}{b}. \end{aligned}$$

• 3. For $1<G<G_c$:

  $$\begin{aligned} L_\mathrm{{eq}}= & {} L_0, G_\mathrm{{eq}}=\frac{(b+\gamma ) _\mathrm{{T}}}{(b+\gamma ) \beta L_0+ I} \quad \text{ exists } \text{ provided } \quad \beta L_0 + \frac{I}{b+\gamma }\\< & {} G_\mathrm{{T}} < G_c \left( \beta L_0 + \frac{I}{b+\gamma } \right) . \end{aligned}$$

• 4. For $G \gg G_c$:

  $$\begin{aligned} L_\mathrm{{eq}}= & {} L_0 \pm \frac{F\mathrm{{max}}}{E}, G_\mathrm{{eq}}=\frac{(b+\gamma ) G_\mathrm{{T}}}{(b+\gamma ) \beta L_\mathrm{{eq}}+ I} \quad \text{ exists } \text{ provided } \quad \\ G_\mathrm{{T}}> & {} G_c \left( \beta L_\mathrm{{eq}} + \frac{I}{b+\gamma } \right) . \end{aligned}$$

In the Rho (contraction) case, tristability still occurs. In this model, we can more readily find parameter sets where the middle ($1<G<G_c$) and high branch persist as $\gamma \rightarrow \infty$, whereas the low branch terminates at a fold point (Fig. 11). (This is also possible in the original model.) This suggests that even with a very large positive feedback in GTPase activity, and when GTPase activity is high enough to trigger that feedback, the cell can still stay at a large size. In contrast, in the original model, the regime where only the middle and high branch coexist is usually bounded above by some maximum $\gamma$, as seen in the regime labeled “2,3” in Fig. 7 (top right).

Fig. 11

Bifurcation diagrams as in Fig. 7, but for the constant volume spatially uniform GTPase model, (16). Top: the contractile (Rho) case, with parameters $k=1,I=5,L_0=3.5,\frac{2E}{\mu +2\eta }=3,\frac{2F\mathrm{{max}}}{\mu +2\eta }=10,G_c=4,\beta =1,G_\mathrm{{T}}=5,m=n=10$. Bottom: the protrusive (Rac) case, with parameters $k=1,I=10,L_0=2.3,\frac{2E}{\mu +2\eta }=8,\frac{2F\mathrm{{max}}}{\mu +2\eta }=55,G_c=0.9,\beta =1,G_\mathrm{{T}}=10,n=20,m=10$. While the Rac case is still tristable, a fold point that previously connected the middle and high branch (in the original model) is now missing. These branches persist as $\gamma \rightarrow \infty$, whereas the lower branch exists only for low $\gamma$. In the Rac case, the Hopf curve in the two-parameter bifurcation diagram (lower right) differs significantly from the corresponding panel in Fig. 7, which yield parameter regimes where the limit cycle exists for $\gamma \rightarrow 0$ (Color figure online)

In the Rac (expansion) case, periodic orbits are still possible. The regime in which limit cycles exist is much wider that in the original model. We can identify parameter sets where a limit cycle exists even at $\gamma =0$, which hints that in this model, a positive feedback in GTPase activity is not necessary for an oscillating cell.

Overall, the model variant with constant cell volume has the same regimes as the original model. However, the location and size of these regimes have shifted.

The Jacobian of the ODE system at a steady state $(L_\mathrm{{eq}}, G_{eq})$, is

$$\begin{aligned} J(L_\mathrm{{eq}}, G_{eq}) = \begin{pmatrix} -\frac{2E}{\mu + 2 \eta } &{} \pm \frac{2F_{\text {max}}}{\mu +2\eta } \frac{mG_c^m G_\mathrm{{eq}}^{m-1}}{\left( G_c^m + G_\mathrm{{eq}}^m\right) ^2} \\ -\frac{2E}{\mu + 2 \eta } \frac{1}{L_\mathrm{{eq}}} &{} \frac{\gamma n G_\mathrm{{T}} G^{n-1}}{\left( 1+G^n \right) ^2} - \beta L\left( b + \gamma \frac{G^n}{1+G^n} + G \right) - I \end{pmatrix}. \end{aligned}$$

Model parameters are all positive, and $L_\mathrm{{eq}}, G_{eq} > 0$ so in the Rac and Rho cases, respectively, the Jacobian has the sign patterns:

$$\begin{aligned} J_{\text {Rac}} = \begin{pmatrix} - &{} + \\ - &{} {\text {sgn}}(d) \end{pmatrix}, \qquad J_{\text {Rho}} = \begin{pmatrix} - &{} - \\ - &{} {\text {sgn}}(d) \end{pmatrix}. \end{aligned}$$

where

$$\begin{aligned} d :=\frac{\gamma n G_\mathrm{{T}} G^{n-1}}{\left( 1+G^n \right) ^2} - \beta L\left( b + \gamma \frac{G^n}{1+G^n} + G \right) - I. \end{aligned}$$

The sign patterns are the same as in the original model, and thus, the same conclusions apply.