## 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.

This is a preview of subscription content, access via your institution.

## References

Bui J, Conway DE, Heise RL, Weinberg SH (2019) Mechanochemical coupling and junctional forces during collective cell migration. Biophys J 117:170–183

Camley BA, Zhao Y, Li B, Levine H, Rappel WJ (2013) Periodic migration in a physical model of cells on micropatterns. Phys Rev Lett 111(15):158102

Crampin EJ, Gaffney EA, Maini PK (1999) Reaction and diffusion on growing domains: scenarios for robust pattern formation. Bull Math Biol 61(6):1093–1120

Cusseddu D, Edelstein-Keshet L, Mackenzie JA, Portet S, Madzvamuse A (2019) A coupled bulk-surface model for cell polarisation. J Theor Biol 481:119–135

Dawes AT, Edelstein-Keshet L (2007) Phosphoinositides and Rho proteins spatially regulate actin polymerization to initiate and maintain directed movement in a one-dimensional model of a motile cell. Biophys J 92(3):744–768

Dierkes K, Sumi A, Solon J, Salbreux G (2014) Spontaneous oscillations of elastic contractile materials with turnover. Phys Rev Lett 113(14):148102

Drasdo D, Höhme S (2005) A single-cell-based model of tumor growth in vitro: monolayers and spheroids. Phys Biol 2(3):133–47

Gerisch A (2001) Numerical methods for the simulation of taxis diffusion reaction systems. Dissertation, Martin-Luther-Universität Halle-Wittenberg

Holmes WR, Edelstein-Keshet L (2016) Analysis of a minimal Rho-GTPase circuit regulating cell shape. Phys Biol 13(4):046001

Holmes WR, Lin B, Levchenko A, Edelstein-Keshet L (2012) Modelling cell polarization driven by synthetic spatially graded Rac activation. PLoS Comput Biol 8(6):e1002366

Jilkine A, Edelstein-Keshet L (2011) A comparison of mathematical models for polarization of single eukaryotic cells in response to guided cues. PLoS Comput Biol 7(4):e1001121

Liu Y (2019) Analysis of pattern formation in reaction–diffusion models for cell polarization. Master thesis, University of British Columbia

Marée AF, Jilkine A, Dawes A, Grieneisen VA, Edelstein-Keshet L (2006) Polarization and movement of keratocytes: a multiscale modelling approach. Bull Math Biol 68(5):1169–1211

Marée AF, Grieneisen VA, Edelstein-Keshet L (2012) How cells integrate complex stimuli: the effect of feedback from phosphoinositides and cell shape on cell polarization and motility. PLoS Comput Biol 8(3):e1002402

Meinhardt H (1999) Orientation of chemotactic cells and growth cones: models and mechanisms. J Cell Sci 112(17):2867–2874

Michaelson D, Silletti J, Murphy G, D’Eustachio P, Rush M, Philips MR (2001) Differential localization of Rho GTPases in live cells: regulation by hypervariable regions and RhoGDI binding. J Cell Biol 152(1):111–126

Mori Y, Jilkine A, Edelstein-Keshet L (2008) Wave-pinning and cell polarity from a bistable reaction–diffusion system. Biophys J 94(9):3684–3697

Mori Y, Jilkine A, Edelstein-Keshet L (2011) Asymptotic and bifurcation analysis of wave-pinning in a reaction–diffusion model for cell polarization. SIAM J Appl Math 71(4):1401–1427

Murray JD (1989) Mathematical biology. Springer, Berlin

Neilson MP, Veltman DM, van Haastert PJ, Webb SD, Mackenzie JA, Insall RH (2011) Chemotaxis: a feedback-based computational model robustly predicts multiple aspects of real cell behaviour. PLoS Biol 9(5):e1000618

Otsuji M, Ishihara S, Kaibuchi K, Mochizuki A, Kuroda S et al (2007) A mass conserved reaction–diffusion system captures properties of cell polarity. PLoS Comput Biol 3(6):e108

Palsson E (2008) A 3-D model used to explore how cell adhesion and stiffness affect cell sorting and movement in multicellular systems. J Theor Biol 254(1):1–13

Peyret G, Mueller R, dAlessandro J, Begnaud S, Marcq P, Mege RM, Yeomans JM, Doostmohammadi A, Ladoux B (2019) Sustained oscillations of epithelial cell sheets. Biophys J 117:464–478

Ryan GL, Watanabe N, Vavylonis D (2012) A review of models of fluctuating protrusion and retraction patterns at the leading edge of motile cells. Cytoskeleton 69(4):195–206

Weiner R, Schmitt BA, Podhaisky H (1997) ROWMAP-a ROW-code with Krylov techniques for large stiff ODEs. Appl Numer Math 25:303–319

Wolgemuth CW, Zajac M (2010) The moving boundary node method: a level set-based, finite volume algorithm with applications to cell motility. J Comput Phys 229(19):7287–7308

Zehnder SM, Suaris M, Bellaire MM, Angelini TE (2015) Cell volume fluctuations in MDCK monolayers. Biophys J 108(2):247–250

Zmurchok C, Holmes WR (2019) Modeling cell shape diversity arising from complex Rho GTPase dynamics. bioRxiv p 561373

Zmurchok C, Bhaskar D, Edelstein-Keshet L (2018) Coupling mechanical tension and GTPase signaling to generate cell and tissue dynamics. Phys Biol 15(4):046004

## Acknowledgements

LEK and coauthors were (partially) supported by NSERC (Natural Sciences and Engineering Research Council) Discovery Grant 41870. YL was partially supported by an NSERC postgraduate fellowship. AB was partially supported by an NSERC PDF Fellowship. We are grateful to the Pacific Institute for Mathematical Sciences for providing space and resources for AB’s postdoctoral research. The authors thank C. Zmurchok and W. R. Holmes and E. Rens for valuable discussions.

## Author information

### Affiliations

### Corresponding author

## Additional information

### Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Appendices

### A Mass Balance on Growing or Shrinking Domain

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

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)]\)

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

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

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

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

### 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

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

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

Canceling two terms on the LHS and rearranging leads to

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

### 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:

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

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.

### 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

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:

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).^{Footnote 1} The ROWMAP integrator (written in Fortran) is wrapped for easy use in a Python environment using f2py^{Footnote 2} (f2py provides a connection between the Python and Fortran languages) and integrated into scipy’s^{Footnote 3} integrator class (a package providing several different techniques to integrate ODEs). The spatial discretization (right hand side of ODE) is implemented using NumPy.^{Footnote 4} 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

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:

where

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\).

### 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

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):

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

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).

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

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:

where

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

## Rights and permissions

## About this article

### Cite this article

Buttenschön, A., Liu, Y. & Edelstein-Keshet, L. Cell Size, Mechanical Tension, and GTPase Signaling in the Single Cell.
*Bull Math Biol* **82, **28 (2020). https://doi.org/10.1007/s11538-020-00702-5

Received:

Accepted:

Published:

### Keywords

- Cell size
- Motility
- GTPases
- Rac and Rho
- Cell polarization
- Reaction–diffusion equations
- Growing 1D domain