# Mass-conservation increases robustness in stochastic reaction-diffusion
models of cell crawling

Eduardo Moreno  Sergio Alonso<EMAIL_ADDRESS>Department of Physics,
Universitat Politècnica de Catalunya, Barcelona, Spain

(August 27, 2024)

###### Abstract

The process of polarization determines the head and tail of single cells. A
mechanism of this kind frequently precedes the subsequent cell locomotion and
it determines the direction of motion. The process of polarization has
frequently been described as a reaction-diffusion mechanism combined with a
source of stochastic perturbations. We selected a particular model of amoeboid
cell crawling for the motion of Dictyostelium discoideum and studied the
interplay between pattern formation and locomotion. Next, we integrated the
model in a two-dimensional domain considering the shape deformations of the
cells in order to characterize the dynamics. We saw that the condition of
pattern formation is finely tuned and we propose a modification based on the
use of a mass-conservation constraint to substantially increase the robustness
of the mathematical model.

bistability, stochastic partial differential equations, pattern formation,
dictyostelium discoideum, cell polarization, amoeboid motion

##  I Introduction

An intensive use of physical and mathematical theories on pattern formation in
extended systems [1, 2], gave rise to valuable arguments to explain the
formation of certain biological structures. Such mathematical mechanisms have
permitted the modeling of processes on very different spatial and temporal
scales: from the formation of the skin in fishes [3], to the definition of the
direction in embryonic developing [4].

Cell migration is a common phenomenon that is present both in prokaryotic and
eukaryotic cells. Living cells migrate to perform different tasks such as food
targeting, wound healing and immune response. Independently of the presence of
an external signal, before moving, cells need to define the direction to
follow. To do so, they first define the front and the back of the cell. The
process of formation of a polar direction inside a single cell is commonly
known as cell polarization [5] and it is a typical example of pattern
formation at the cellular level [6].

Several mathematical models have been developed to explain polarization of
single eukaryotic cells [7, 8, 9]. Some models rely on a local excitation,
which, combined with global inhibition makes the cell respond to external
gradients [10], while others rely on the accumulation of a certain biochemical
components to guide the motion of the single cells [11]. Normally the
accumulation responsible for this second mechanism is combined with a
conservation constraint because the process of polarization is fast in
comparison to the production of new biochemical components. With this
restriction, the mechanism of pattern formation inside living cells together
with the constraint of mass conservation is analogous to a process of
coarsening of the initial nucleus of components [12] and gives rise to phase
separation [13] and to models of pattern formation [14]. Several simple models
of intracellular pattern formation have appeared that include the conservation
restriction [11, 15, 16].

Once the axes and direction of movement are defined, in amoeboid cells small
projections (defined as protrusions) are formed in the membrane of the cells
[17]. These projections, which extend and retract periodically, are
responsible for pushing the cell to move in the typical amoeboid motion.
Dictyostelium discoideum is a species of soil-dwelling amoeba, frequently
employed to characterize amoeboid locomotion. Inside the projections of this
amoeba, several signaling events are triggered, for example the activation of
the Ras proteins and PI3K enzymes and accumulation of $PIP_{3}$ at the front
of the cell, while activation of PTEN and myosin occur at the rear of the cell
[18, 19]. Membrane areas where the protrusion activity is greater are
typically characterized by the presence of Ras-GTP protein zones, denoted as
patches [20]. The appearance of Ras regions around the membrane with high
protrusion activity [10, 21, 22] has been observed in the slime mold organism
Dictyostelium discoideum where it is related to cytoskeletal dynamics [23].

A reaction-diffusion model with bistable dynamics is one of the common models
of cell motility and particularly employed for Dictyostelium discoideum. Some
dimensional bistable models are based on the formation of a finite lifetime
and localized patches of high protein concentration [24, 25] while others are
based on membrane dynamics of moving connected points [26]. Bistable
conditions for cellular processes can be obtained by the coupling of a mass
control regulating condition in the biochemical components present around the
cytosol and the membrane of the cell, for example proteins, phospholipids and
enzymes present or other substances in the cytoskeleton [14, 11, 27].

A common technique to model pattern formation inside a cell and the shape
evolution of the membrane is the addition of a phase field with a sharp
interface to distinguish the interior and exterior of the cell. This field
maintains the correct boundary conditions while the borders are moving [28].
Some studies have applied phase field modeling to study keratocyte motility
[29, 30, 31] and amoeboid motility [32] which can be divided into diffuse and
persistent migration depending on the starvation level [33, 34]. There are
also some intermediate cases observed in Dictyostelium discoideum cells for
certain types of genetic variants [35] which have also been modelled with a
phase field and variation of some parameters [36, 37, 38]. Other properties of
Dictyostelium discoideum cells modeled employing a phase field are the
viscoelasticity of the cells [39] and cell division [40]. Finally, we would
like to mention that interactions among cells can also be considered for
keratocyte [41, 42] and amoeboid cells [43].

Here, first, we transform a one dimensional model of the polarization at the
membrane of a single Dictyostelium discoideum cell [21] into a two-dimensional
domain for the waves in the basal membrane, in contact with the surface, using
an additional phase field for the shape of the cell. We observe a strong
dependence of the numerical dynamics on the explicit parameter values of the
original model and on the intensity of the stochastic fluctuations. Such
dependence indicates that the model is not robust to small variations of the
parameters. Therefore, we propose a constraint on the conservation of a
component of the signal pathway, controlling the autocatalytic mechanism, in
order to increase the robustness of the model in representing changes in
parameter values [44]. Similar types of conservation have previously been
employed in other models of motion of Dictyostelium discoideum cells [32, 36,
45, 46], and we show here that a constraint of this kind is an useful and
reasonable condition to systematically increase the robustness of the crawling
mechanism, increasing the window of parameter values allowing cell migration.

##  II Model and Methods

###  II.1 Biochemical model for Ras activation and pseudopod extension

Biochemical components inside the Dictyostelium discoideum amoeba self-
organized to follow chemical signals in the exterior by the accumulation of
Ras-GTP in the front of the cell. Next, F-actin molecules accumulates also in
the front triggering push of the cytoskeleton and the motion of the amoeba.

We investigate a simple reaction-diffusion model [21] for Ras-GTP (R) patches
that consists of a local activator and a global inhibitor enhanced by the
activation of the Ras-GTP; see diagram in Fig. 1. The partial differential
equations employed are:

$\displaystyle\frac{\partial R}{\partial t}$
$\displaystyle=(1-R)(k_{1}+k_{2}\frac{R^{n1}}{R^{n1}+K_{R}^{n1}}-k_{3}G_{R}-k_{4}L_{R})$
(1) $\displaystyle-
k_{5}R\frac{1}{1+k_{6}P}+D_{R}\nabla^{2}R+\xi_{R}(\vec{x},t),$

$\frac{\partial G_{R}}{\partial t}=k_{7}\langle R\rangle-k_{9}G_{R},$ (2)

$\frac{\partial L_{R}}{\partial
t}=(1-L_{R})k_{10}R-k_{11}L_{R}+D_{LR}\nabla^{2}L_{R};$ (3)

where $\xi_{R}(\vec{x},t)$ is a Gaussian spatio-temporal distributed white
noise with zero mean $\langle\xi_{R}(\vec{x},t)\rangle=0$ and correlation
$\langle\xi_{R}(\vec{x},t)\xi_{R}(\vec{x^{\prime}},t^{\prime})\rangle=2\sigma_{R}\delta(\vec{x}-\vec{x^{\prime}})\delta(t-t^{\prime})$.
The variables $L_{R}$ and $G_{R}$ correspond, respectively, to local and
global inhibitors of R. Finally, note that the quantity $\langle R\rangle$
corresponds to the spatial integration in the whole system of the field R.

Figure 1: Sketch of the biochemical model of polarization and locomotion
modules. Blue solid lines corresponds to the interactions and red dashed lines
to auto-catalytic processes. The mass-conservation constrains control the
auto-catalytic interactions.

At the same time we included a quantity related to the formation of
protrusions (P) such as F-actin and Rac-GTP molecules; see diagram in Fig.1.
These variables were coupled with their respective inhibitors:

$\displaystyle\frac{\partial P}{\partial t}$
$\displaystyle=(1-P)(k_{12}+k_{13}R+k_{14}\frac{P^{n2}}{P^{n2}+K_{P}^{n2}}$
$\displaystyle+k_{15}\frac{M^{n3}}{M^{n3}+K_{M}^{n3}}-k_{16}G_{P}-k_{17}L_{P})-k_{18}P$
$\displaystyle+D_{P}\nabla^{2}P+\xi_{P}(\vec{x},t),$

$\frac{\partial G_{P}}{\partial t}=k_{19}\langle P\rangle-k_{20}G_{P},$ (5)

$\frac{\partial L_{P}}{\partial
t}=(1-L_{P})k_{21}P-k_{22}L_{P}+D_{LP}\nabla^{2}L_{P};$ (6)

where $\xi_{P}(\vec{x},t)$ is a Gaussian spatio-temporal distributed white
noise with zero mean $\langle\xi_{P}(\vec{x},t)\rangle=0$ and correlation
$\langle\xi_{P}(\vec{x},t)\xi_{P}(\vec{x^{\prime}},t^{\prime})\rangle=2\sigma_{P}\delta(\vec{x}-\vec{x^{\prime}})\delta(t-t^{\prime})$.
As in the previous set of equations, the variables $L_{P}$ and $G_{P}$
correspond, respectively, to local and global inhibitors of P; while the
quantity $\langle P\rangle$ corresponds to the spatial integration in the
whole system of the field P.

Finally, the model takes into account the inclusion of a variable of memory
(M) which is coupled with P; this variable stimulates the formation of new
protrusion zones and represents the results observed in some experiments where
a certain relation between the formation of new pseudopods in previous
pseudopod locations is present [47, 48]:

$\frac{\partial M}{\partial t}=k_{23}P-k_{24}M+D_{M}\nabla^{2}M.$ (7)

The memory term is related to the probability of identifying when and where
the signaling cascade is activated to generate new pseudopods. At a molecular
level, such term is identified with the mechanisms of how information is
collected, stored and used to bias future pseudopods [48].

For a more exhaustive description of the model, check the original study [21].
Note that the noise description in the original study was different and we
have adapted to an equivalent description based on physical derivations of
stochastic fluctuations [49]; for more details see appendix A.

One dimensional simulations were made using periodic boundary conditions and
using 120 points for the grid, see an example in Fig.(2,A) where R and P are
plotted as function of time. The cell was considered as circular and having a
radius of $6.25$ $\mu m$. The pixel size for this case was set at $\Delta
x=0.32$ $\mu m$ and the time step $\Delta t=0.03s$ The definition and the
value of the parameters of the model can be found in Table 1 in appendix B.

###  II.2 Physical phase field model for cell shape deformations

The original model did not include deformable cells. Therefore, we expanded
the model to 2D geometry, see Fig.(2,B) for a simple extension, and coupled it
to an auxiliary phase field $\phi$ with the purpose of describing the
evolution of the cell shape. The use of a phase field permits a smooth
variation between the values of $\phi=1$ inside and $\phi=0$ outside of the
cell, respectively. See an example in Fig.(2,C) where the field P deforms the
shape. Note that the fluctuations in P will also affect the shape, see
Fig.(2,D).

The phase field equation is the result of a force balance involving several
types of forces of different natures acting in the cell body. The equation for
the phase field is as follows

$\displaystyle\tau\frac{\partial\phi}{\partial t}$ $\displaystyle=$
$\displaystyle\gamma\left(\nabla^{2}\phi-\frac{G^{\prime}(\phi)}{\epsilon^{2}}\right)$
$\displaystyle-\beta\left(\int\phi\,dA-
A_{0}\right)\left|\nabla\phi\right|+\alpha\,\phi\,P\left|\nabla\phi\right|\,,$

The first term on the right side in eq. (II.2) is related to the surface
energy of the cell membrane, where $\gamma$ is the surface tension and
$G(\phi)=18\phi^{2}(1-\phi)^{2}$ is a double well potential.The second term
helps keep the area close to the value of $A_{0}$. And the third term
represents the active force generated by the Rac-GTP (P) molecules when
pushing on the cell membrane [32].

The inclusion of the phase field as previously mentioned implies a change of
geometry from 1D to 2D. Here we consider a deformable cell of $122\mu m^{2}$
in area corresponding to a circle with radius equal to $6.25\mu m$. The pixel
size used was the half as in the 1D case $\Delta x=\Delta y=0.16\mu m$. Also,
to increase accuracy, the time discretization was reduced to $\Delta
t=0.003s$. We integrated Eqs. (1-II.2) using periodic boundary conditions and
standard finite differences. The rest of parameters are as displayed in Table
1 in appendix B.

Figure 2: Explanatory box of the modelling techniques. A) Kymograph of a
typical one dimensional simulation showing the formation of patches of R and P
(horizontal direction corresponds to position in the membrane, and vertical
direction to the time evolution ). B) Extension to two dimensions with the
addition of a circular phase field (Value 1 inside and 0 outside of the
circular border). Patches of R and P form in the basal surface of the cell. C)
Addition of a dynamic phase field producing the motion and the deformation of
the cell. D) Addition of the fluctuations, random two dimensional patches are
distributed inside the phase field, driving the motion of the cell.

###  II.3 Mathematical description of the mass-conservation constraint

The mathematical model represented in Eqs. (1-7) has several conservation
terms which affect the dynamics of the system; see for example Eq.(2) and
Eq.(5).

However, we modify the original model to include two new feedback through
parameters $k_{2}$ and $k_{14}$ to create a more robust model by the control
of the bistability properties of the model, which are discussed in following
sections. Therefore, the parameters $k_{2}$ and $k_{14}$ are dynamically
controlled depending on the total amount of the protein R on inducer P at the
cell, respectively.

The new control term related to $k_{2}$ is

$\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle
k_{2}^{*}+\eta\left(\langle R\rangle-C_{R}\right)\,,$ (9)

Also, in a similar way, the control term for $k_{14}$ is as follows:

$\displaystyle k_{14}$ $\displaystyle=$ $\displaystyle
k_{14}^{*}+\eta\left(\langle P\rangle-C_{P}\right)\,,$ (10)

where $k_{2}^{*}$ and $k_{14}^{*}$ are new constants that take the values of
the original value of $k_{2}$ and $k_{14}$, respectively. Note that for
$\eta$=0 we recover the original model.

Parameter $C_{R}$ is the target value of the fraction of the cell area
occupied by patches of Ras $R$. When, the total concentration $\langle
R\rangle$ is larger (lower) than $C_{R}$ the parameter value of $k_{2}$is
larger (lower) than $k*_{2}$ and the bistable conditions implies the reduction
(increase) of the R and, therefore, of $\langle R\rangle$. It corresponds to a
feedback control of the parameter to keep a particular stable solution in
Eq.(9). The parameter $C_{P}$ is the target fraction occupied by the pseudopod
inducer $P$ and produces an equivalent dynamics on Eq.(10).

##  III Results

###  III.1 Mechanism of cell motion in the mathematical model is based in
stochastic generation of patches

In a one-dimensional system, mimicking the membrane of a crawling cell, the
generation of a local patch of high biochemical concentration is equivalent to
cell polarization. A stable domain in a specific location of the membrane is
related with actin accumulation and produces persistent motion of the cell in
that direction. The random appearance and disappearance of small domains is
related with amoeboid motion, where two or three pseudopods compete for a
certain time. The alternation of direction gives rise to random motion of the
virtual cell.

Stochastic reaction-diffusion equations, as in Eqs.( 1-7), and previously
developed [21], can be numerically integrated into a one-dimensional domain as
in the original study. Kymographs are shown next for a certain window of
parameter values, the spatio-temporal dynamics of concentration R displays
relatively stable patches which appear at different locations; see spatio-
temporal plots in Fig. 3.

In the simulations shown in Fig. 3, the parameters $k_{14}$ and $k_{18}$ are
kept constant and the response of P along the membrane is quite similar for
the parameters considered. There are some persistent patches during a long
period of the simulations. Under such conditions the cell may develop its own
locomotion. Note that the activation of P by R is weak; a strong coupling,
controlled by parameter $k_{13}$, would correlate the patches of both fields.
We keep the original values of this parameter [21]. For small values of the
parameter $k_{5}$ and large values of $k_{2}$, the concentration of R is high
and homogeneously distributed along the membrane; see corresponding spatio-
temporal plots in Fig.3. The concentration of R increases the probability of
the generation of patches of P, and therefore, depending on the particular
realization, some direct motion is expected under such conditions. On the
other hand, in the opposite limit, large values for parameter $k_{5}$ and
small values for $k_{2}$, the concentration of R is low, and the probability
of directed motion is smaller.

If we fix $k_{2}$ and $k_{5}$ values and vary the parameters controlling for
the evolution of P, we obtain the inverse behaviour; see Fig.4. For large
values of $k_{14}$ and small values of $k_{18}$ the membrane is completely
covered by P, while for small values of $k_{14}$ and large values of $k_{18}$
the concentration of P strongly decreases. Locomotion is, therefore, expected
for intermediate values of these two parameters as shown in Fig. 4. However,
the effect of changing such parameters on the dynamics of R is small.

Figure 3: Spatio-temporal plots of R (in green) and P (in red) for different
values of $k_{5}$ and $k_{2}$. For each panel times goes from bottom to top
and horizontal direction corresponds to membrane position in the cell
circumference. Variances of noise intensity were fixed at $\sigma_{R}=0.04$
and $\sigma_{P}=0.025$. Figure 4: Spatio-temporal plots of R (in green) and P
(in red) for different values of $k_{14}$ and $k_{18}$. For each panel times
goes from bottom to top and horizontal direction corresponds to membrane
position in the cell circumference. Variances of noise intensity were fixed at
$\sigma_{R}=0.04$ and $\sigma_{P}=0.025$.

###  III.2 Mechanism of amoeboid motion is based in bistability

For the evaluation of the mechanism of the spatio-temporal plots in Figs. 3
and 4, we fixed the stochastic source of perturbations of the additive noises
in eqs (1-7) to $\sigma_{R}=0.04$ and $\sigma_{P}=0.025$. No changes were
expected in the average behaviour or in the quantification the value of the
total concentration R in the face of the parameters $k_{2}$ and $k_{5}$
because the fluctuations are additive [49]. Stationary solutions were found
for the deterministic version of Eqs. (1-7). Slowly increasing the parameter
$k_{2}$ we observed the evolution of the stationary value of R; see Fig. 5.
Once the system of equations saturated to a certain value of R$\sim 0.85-0.90$
we reduced the value of $k_{2}$, giving rise to a hysteresis cycle; see Fig.
5A. The hysteresis cycle shows that the system actually follows bistable
dynamics. Equivalent dynamics of R can be observed under a similar change in
$k_{5}$; see Fig. 5B.

In contrast, the same analysis was made for the constants $k_{14}$ and
$k_{18}$, which are related to the variable of the pseudopod inducer P. The
resulting dynamics revealed bistable behavior for two stable states as shown
in the hysteresis curve for parameters $k_{14}$ (Fig. 5C) and $k_{18}$ (Fig.
5D).

Figure 5: Hysteresis curves of the deterministic model for R dynamics varying
(A) $k_{2}$ and (B) $k_{5}$. Similar hysteresis curves of the deterministic
model for P dynamics were obtained when varying (C) $k_{14}$ and (D) $k_{18}$.

It is known that the combination of a bistable dynamics with an appropriate
noise intensity, produces the formation of localized patches in reaction-
diffusion equations [50, 51]. This formation is the mechanism responsible for
the formation of localized domains of R and P. The localized pattern observed
in Fig. 3 is not due to excitable dynamics but rather to bistable dynamics
with noise.

###  III.3 Stochastic bistability provides a mechanism of cell crawling motion
in two dimensional systems

For the coupling of the polarization mechanism to the proper cell motion shown
in the previous section we used a phase field. This additional field is
employed for the definition of the interior of the cell. In this case, the two
dimensional phase field permits the use of different biochemical
concentrations in the ventral part of the cell in a surface.

The extension of the stochastic reaction-diffusion described in the previous
section to two dimensions permits numerical simulations of the shape of the
crawling cell and the corresponding motion responding to the dynamics of the
patches. In Fig. 6 we see snapshots of the in silico cells with different
parameter values corresponding to equivalent parameter values shown in Fig. 3.
We reproduced the dynamics expected from the one dimensional simulations and
the crawling dynamics can be clearly studied. Again, for small values of
parameter $k_{5}$ and large values of $k_{2}$ the concentration of R is mostly
high and homogeneously distributed along the membrane; see Fig. 6. On the
other hand, in the opposite limit, for large values of $k_{5}$ and small
values of $k_{2}$ the concentration of R disappears from the ventral membrane.
However, for a small window of values of the parameters the cell moves and
inspects the surrounding region. The shape of the resulting cell depends on
the particular realization and the parameter values; in the snapshots shown in
Fig. 6 we obtain fan-shaped cells, a typical mode of Dictyostelium discoideum
cell crawling [35].

Figure 6: Snapshots of the virtual cells showing R (in green) and P (in red)
obtained from computer simulations and applying the phase field technique to
the model for different values of $k_{2}$ and $k_{5}$. Variances of noise
intensity were fixed at $\sigma_{R}=0.04$ and $\sigma_{P}=0.025$.

From the snapshots resulting from Fig. 6 different scenarios show different
types of movement together with the concentrations of R and P. The model
exhibits sensitivity to variations in parameters $k_{2}$ and $k_{5}$.
Different realizations of simulations with the same parameters could bring the
cell from quite persistent movement into an alternation between an almost
static state for low concentrations of $\langle R\rangle$, and transitory
motion in random directions for large values of $\langle R\rangle$. As in
Fig.3 the small influence of R in P is evident in these two-dimensional
simulations. Note that the presence of patches of P appearing for high values
of $k_{2}$ and low values of $k_{5}$ is more common because $\langle R\rangle$
slightly increases the probability of forming new $\langle P\rangle$ patches.
In contrast, for low values of $k_{2}$ and high values of $k_{5}$ the presence
of $\langle P\rangle$ patches is less frequently observed.

In contrast to this, the dependence of locomotion on the parameters $k_{14}$
and $k_{18}$ is much stronger as shown in Fig.7, where persistent motion is
observed only for intermediate values of the two parameters. For very small or
very large amounts of $\langle P\rangle$ cells do not move, independently of
the behaviour of R.

Figure 7: Snapshots of the virtual cells showing R (in green) and P (in red)
obtained from computer simulations and applying the phase field technique to
the model for different values of $k_{14}$ and $k_{18}$. Variances of noise
intensity were fixed at $\sigma_{R}=0.04$ and $\sigma_{P}=0.025$.

###  III.4 Noise intensities determine the type of motion

The variation of the type of cell motion is determined by the intensity of the
noise. Low levels of both noises do not permit the formation of the domains of
R and P needed to give rise to cell movement; see Fig.8. As we increase the
amplitude of noise for P, we reach a single domain that fills the front part
of the cell and persistent motion is shown; see middle column and right hand
snapshots in Fig.8. This type of motion is reminiscent of the fan-shaped
amoeboid cells previously reported [35] and also reproduced with simpler
models [36]. An increase in the noise intensity produces an increase in the
appearance and disappearance of pseudopods and therefore a transition to
amoeboid movement; see fifth column in Fig.8 for higher noise intensities.
Another evident change the effects of which will be discussed in the following
sections is the variance associated with the formation of Ras patches
($\sigma_{R}$) becoming varied. There are some noise intensities which produce
patterns and dynamics similar to the characteristic patterns observed in the
experiments, however, as previously mentioned, and shown in Fig.10, small
variations in the parameter values may completely change the movement
dynamics.

Figure 8: Map of snapshots of the virtual cell taken from different variance
values of noise intensity for the original model [21] . For every column,
left/green snapshots refer to R distribution while right/red snapshots
indicate P distribution. The value of the parameter of the simulations are
$k_{2}=9.50$, $k_{5}=7.0$, $k_{14}=10.20$ and $k_{18}=7.0$. The rest of them
are consistent with Table 1.

###  III.5 Inclusion of mass-conservation stabilizes cell motion

In order to reduce the high sensitivity of the movement on the parameter
values, we included a conservation constraint on the total quantities of the
protein Ras R and the inducer P in the biochemical rates responsible for the
bistable transition. Following similar previous approaches [32, 45], we
included this conservation as a global feedback condition in the reaction-
diffusion equations.

In Fig. 9 several realizations incorporating the conservation constraint of R
and P proteins, for different parameter values are shown; see Eqs.9 and 10. A
direct comparison can be made with the results displayed in Fig. 6. Results in
Fig. 9 show the effects of the inclusion of the mass conservation feedback,
where a wide change in the parameter values does not affect either the
bistability of the model or the quantity of protein R inside the cell phase
field, which was kept constant. The inclusion of this global condition permits
the use of a larger range of parameter values, increasing the robustness of
the mechanism in comparison with the original model.

As with the cases without mass conservation, the shape and type of motion were
affected by the noise intensity. Using the same parameters of Fig. 8 we
studied the effects of varying the noise variance in the new mass-conserved
model. In Fig. 10 we see two different types of transitions. The first is
observed when the variance related to the pseudopod inducer P increases,
changing the shape from persistent fan-shape to random amoeboid phenotype.
This change is present because while small values of P variance generates less
blurred patches, an increment in the P variance drives the opposite effect,
which is the appearance of more blurred and distributed patches for P.

For a second transition, we follow the same comparison line between the
simulations without and with mass conservation from Fig. 8 and 10,
respectively; now we analyze the effects for different values of the variance
corresponding to the dynamics of the Ras concentration R. Similar to the
patterns of P, small values of $\sigma_{R}$ induce the formation of less
diffuse and distributed Ras patches. In contrast, greater values of
$\sigma_{R}$ stimulate more diffuse and distributed Ras patches. Despite the
close relationship between the effects of $\sigma_{R}$ and $\sigma_{P}$, the
real difference is in the dynamics of the movement. For example, if we focus
on the top of the first column of Fig. 10 for small values of the noise
variances we see for both snapshots a domain in the front part of the cell
which translates into persistent motion. Now, if we look in detail at the
bottom of the first column corresponding to higher a value for the variance
$\sigma_{R}$, we see a notable difference in the left hand snapshots, moving
from a more compact domain in the front part of the cell to a more diffuse
domain inside the cell. This transition only takes place at the Ras
concentration (R) level, while maintaining consistent dynamics of the
pseudopod inducer P and thus persistent motion of the cell. The same analogy
could be made for the last column in Fig. 10 where high values of $\sigma_{P}$
lead to a more random motion for all the cases regardless of the value of
$\sigma_{R}$, the effect of which is to promote the appearance of less and
more distributed Ras patches but without affecting the dynamics.

Figure 9: Map of snapshots of the virtual cell taken from different values of
$k_{2}$ and $k_{5}$ for the new mass-conserved model. For every column,
left/green snapshots refer to R distribution while right/red snapshots
indicate P distribution. Variances of noise intensity were fixed at
$\sigma_{R}=0.04$ and $\sigma_{P}=0.025$. Figure 10: Map of snapshots of the
virtual cell taken from different variance values of noise intensity for the
new mass-conserved model. For every column, left/green snapshots refer to R
distribution while right/red snapshots indicate P distribution. The value of
the parameter of the simulations are $k_{2}=9.50$, $k_{5}=7.0$, $k_{14}=10.20$
and $k_{18}=7.0$. The rest of them are consistent with Table 1.

###  III.6 Model with mass-conservation constraint is more robust

The use of the mass conservation condition has been previously employed for
the generation of stable pattern formation on the polarization process in
single cells. We employed a version of the mechanism to modify previous
reaction-diffusion equations to incorporate restrictions on the conservation
of R and P. We show that this condition increases the window of parameters of
the model to develop typical dynamics of cell crawling.

For the characterization of the mechanism we calculated the speed of the
resulting cells for different parameter conditions for the original and the
new mass-conserved models. We observed that the dynamics of the simulated
cells are more robust for the mass-conserved version of the equations; see
Fig. 11.

For the original model the speed and shape drastically changes when varying
the parameters, and every simulation reproduces different dynamics of crawling
cells (see Fig. 6). Furthermore, for the new mass-conserved version both
velocities and shapes are more similar to each other for a much larger window
of parameter values, as we can see in the speeds calculated in Fig. 11.
Results displayed in Fig. 11A were obtained by changing $k_{5}$ to correspond
to simulations shown in the third row of Figs. 6 and 9. The dynamics are fixed
because of the conservation condition. In contrast, Fig. 11B was obtained by
changing $k_{18}$ to correspond to simulations shown in the third row of Figs.
7 for the no mass conservation case (results not shown for the mass-
conservation system). Here one can see the preserving of the speed in the new
mass-conserved system in comparison with the original model due to the small
window of parameters, without the mass conservation condition that shows
crawling dynamics.

Figure 11: Box plot representation for the cell speed of the original model
(blue color) and the new mass-conserved model (red color). Results are for
when A) $k_{5}$ and B) $k_{18}$ are varied. For both cases it was used a fixed
value of $k_{2}=9.5$ and $k_{14}=10.2$, $\sigma_{R}=0.04$ and
$\sigma_{P}=0.025$. Ten realizations for each case were performed.

In Fig. 12 we see the dependency of the speed on the noise $\sigma_{P}$ for
different values of $\sigma_{R}$. The speed of the cells employed for the
calculation comes from the results in Figs. 8 and 10. For the outcome of the
original model, which is displayed in blue bars, we do not observe a clear
tendency in the bars. However, in contrast, in the new mass-conserved model
results, represented by red bars, a clearer tendency in the bars is noted.

We observe that the speed obtained from the new mass-conserved model decreases
for greater noise intensities. However, for the set of blue bars almost null
velocities appear for low noise magnitude, and as we increase the noise
variance, different values close to the magnitude of the red bars are reached,
thereby showing a similarity in the speed measurements as the variance is
increased.

Figure 12: Speed bar plots measured by the original model (blue bars) and the
new mass-conserved model (red bars) by varying the noise variance
$\sigma_{P}$. In each panel the variance associated to $\sigma_{R}$ was fixed
for different values.

##  IV Discussion

We studied the dynamics of a one-dimensional Dictyostelium Discoideum cell
model which is mainly made up of coupled Ras activation (R), pseudopod inducer
(P) and stochasticity. The deterministic model by itself reveals the emergence
of bistability. This is from the visible hysteresis transition under the
change of a set of parameters. In contrast, the model with the stochastic part
permits the transition between two stable regimes which are induced by the
parameters and the noise intensity. The extension of the model into two
dimensions with the use of the addition of a phase field, allowing cell
deformation, shows similar behaviour when parameters are varied to the one-
dimensional model, in agreement with [21].

Since for most of the realistic scenarios the mass barely changes despite
environmental perturbations, we chose to incorporate global feedback for the
Ras and pseudopod inducer inside the cell. This increases the robustness of
the model and expands the parameter range without affecting the shape and
dynamics.

Note that other global feedback were already considered in the original
version of the model [21]; see the dynamics of the global inhibitors, defined
by GR and GP in Eq.(2) and Eq.(5) respectively, similarly to other models of
the motion of Dictyostelium discoideum cells based on global cytosolic
quantities [46, 10]. The addition of global control quantities which can
determine the available concentration inside the cell has previously been
employed in other models [30, 36, 45].

Instant speed was measured, highlighting that for the original model a
considerable set of parameters driven to velocities of close to zero, while
for the set of parameters for the new mass-conserved model the speed remained
almost constant.

Note that we have shown that the original model is based on a bistable
conditions with the appropriate noise intensity. We have previously used
simple stochastic bistable models for the description of the motion of
Dictyostelium discoideum cells [32, 36] based on the same concept together
with a mass-conservation constraint. There are also other models based on
similar constraints [45]. Therefore, this constraint is a convenient mechanism
to increase robustness in the delicate equilibrium between stochastic
fluctuations and a bistable condition.

The coupling of R on the pseudopod inducer P is weak. The coupling is produced
by the term proportional to $k_{13}$ and it increases the probability of the
generation of patches of P. A chemotactic concentration can enhance R and,
therefore, polarize the field P and direct the locomotion [21]. In the absence
of chemotaxis the coupling is small, while more physiological conditions may
require larger couplings. In Fig.13 we increased parameter $k_{13}$ to
demonstrate the increase in the coupling of both field R and P through the
parameter. Intermediate values can be compared with experimental measures to
fit an appropriate value; however such evaluation is outside the scope of this
work.

Figure 13: Spatio-temporal plots of R (in green) and P (in orange) for
different values of the coupling parameter between R and P: $k_{13}=0.1$ (A),
$k_{13}=0.50$ (B), $k_{13}=1$ (C), and $k_{13}=2$ (D) in the new mass-
conserved model. Results are for a fixed values of $\sigma_{R}=0.04$ and
$\sigma_{P}=0.025$.

An extension of the new mass-conserved model is possible with the inclusion of
an external chemo-attractant gradient to analyze the response, known as
chemotaxis and already considered in the original model [21]. Future
simulations should study the relation of the chemotactic motion of
Dictyostelium discoideum cells with the stochastic fluctuations and the rest
of the parameters. They should also measure the effects on the shape and speed
caused by the gradient.

In summary, we have shown that the inclusion of a mass-conservation constraint
in the correct position substantially increases the robustness of a
computational model of a crawling cell.

## Acknowledgments

We thank grant PGC2018-095456-B-I00 funded by MCIN/AEI/10.13039/501100011033
and by “ERDF A way of making Europe”, by the “European Union”. E.M.
acknowledges also financial support from CONACYT.

##  Appendix A Integration of stochastic partial differential equations

The Gaussian white noises employed in the model described here are always
additive, and therefore it is enough to use a Euler-Maruyama method for the
integration of the stochastic partial differential equations (SPDE) considered
here[49]. Assuming a generic version of the equations considered here:

$\displaystyle\frac{\partial C(\vec{x},t)}{\partial
t}=F(C(\vec{x},t))+\xi(\vec{x},t),$ (11)

the standard method to integrate this SPDE consists in the discretization in
space ($\Delta\vec{x}$): from $C(\vec{x},t)$ to $C_{i}(t)$, and time ($\Delta
t$) and to account for the particularities of the noise terms:

$\displaystyle C_{i}(t+\Delta t)=C_{i}(t)+\Delta
tF(C_{i}(t))+\sqrt{\frac{2\sigma^{2}\Delta t}{\Delta x^{d}}}N(0,1),$ (12)

where d is the spatial dimensionality of the domain considered (here d=1 or
d=2), and N(0,1) corresponds to a random number generated by a normal
distribution of mean zero and variance 1.

Such standard method for the integration of SPDE is equivalent under some
conditions with the version employed in the original publication of the model
[21] employed here:

$\displaystyle C_{i}(t+\Delta t)=C_{i}(t)+\Delta
tF(C_{i}(t))+N(0,\sigma_{2}),$ (13)

which under correct tuning of the parameter $\sigma_{2}$ can give rise to an
equivalent description.

##  Appendix B Parameters and constants used in the model.

The meaning of the biochemical parameters and of the parameters associated to
the phase field can be obtained, respectively from the original studies [21]
and [36].

Note also that the parameters related to a chemoattractant gradient [21] were
excluded in this study.

Parameter | Value | Parameter | Value
---|---|---|---
$k_{0}$ | 0.01 $s^{-1}$ | $D_{R}$ | 0.75 $\mu m^{2}/s$
$k_{1}$ | 1.5 $s^{-1}$ | $D_{LR}$ | 0.25 $\mu m^{2}/s$
$k_{2}-k_{2}^{*}$ | 7.5-11.5 $s^{-1}$ | $D_{P}$ | 0.75 $\mu m^{2}/s$
$k_{3}$ | 1.5 $s^{-1}$ | $D_{LP}$ | 0.25 $\mu m^{2}/s$
$k_{4}$ | 1.0 $s^{-1}$ | $D_{M}$ | 0.25 $\mu m^{2}/s$
$k_{5}$ | 6.0-8.0 $s^{-1}$ | $n_{1}$ | 2.0
$k_{6}$ | 0.1 $s^{-1}$ | $n_{2}$ | 2.0
$k_{7}$ | 0.05 $s^{-1}$ | $n_{3}$ | 2.0
$k_{9}$ | 0.15 $s^{-1}$ | $K_{R}$ | 0.4
$k_{10}$ | 0.1 $s^{-1}$ | $K_{P}$ | 0.4
$k_{11}$ | 0.1 $s^{-1}$ | $K_{M}$ | 0.4
$k_{12}$ | 0.1 $s^{-1}$ | $\sigma_{R}$ | $0.02-0.06$
$k_{13}$ | 0.1 $s^{-1}$ | $\sigma_{P}$ | $0.015-0.055$
$k_{14}-k_{14}^{*}$ | 9.2-11.2 $s^{-1}$ | $\tau$ | 2.0 $pNs\mu m^{-2}$
$k_{15}$ | 0.5 $s^{-1}$ | $\gamma$ | 2.0 $pN$
$k_{16}$ | 1.5 $s^{-1}$ | $\epsilon$ | 0.75 $\mu m$
$k_{17}$ | 1.8 $s^{-1}$ | $\beta$ | 18.64 $pN\mu m^{-3}$
$k_{18}$ | 6.0-8.0 $s^{-1}$ | $A_{0}$ | 122.7 $\mu m^{2}$
$k_{19}$ | 0.1 $s^{-1}$ | $\alpha$ | 10.0 $pN\mu m^{-1}$
$k_{20}$ | 0.15 $s^{-1}$ | $C_{R}$ | 30.6 $\mu m^{2}$
$k_{21}$ | 0.2 $s^{-1}$ | $C_{P}$ | 30.6 $\mu m^{2}$
$k_{22}$ | 0.1 $s^{-1}$ | $\eta$ | $2.0s^{-1}$
$k_{23}$ | 0.1 $s^{-1}$ | $\Delta t$ |  $0.03-0.003$ $s$
$k_{24}$ | 0.03 $s^{-1}$ | $\Delta x$ |  $0.32-0.16$ $\mu m$

## References

  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Reviews of modern physics 65, 851 (1993).
  * Koch and Meinhardt [1994] A. Koch and H. Meinhardt, Reviews of modern physics 66, 1481 (1994).
  * Kondo and Miura [2010] S. Kondo and T. Miura, science 329, 1616 (2010).
  * Gross  _et al._ [2019] P. Gross, K. V. Kumar, N. W. Goehring, J. S. Bois, C. Hoege, F. Jülicher, and S. W. Grill, Nature Physics 15, 293 (2019).
  * Mogilner  _et al._ [2012] A. Mogilner, J. Allard, and R. Wollman, Science 336, 175 (2012).
  * Camley and Rappel [2017] B. A. Camley and W.-J. Rappel, Journal of physics D: Applied physics 50, 113002 (2017).
  * Rappel and Edelstein-Keshet [2017] W.-J. Rappel and L. Edelstein-Keshet, Current opinion in systems biology 3, 43 (2017).
  * Jilkine and Edelstein-Keshet [2011] A. Jilkine and L. Edelstein-Keshet, PLoS Comput Biol 7, e1001121 (2011).
  * Beta and Kruse [2017] C. Beta and K. Kruse, Annual Review of Condensed Matter Physics 8, 239 (2017).
  * Xiong  _et al._ [2010] Y. Xiong, C.-H. Huang, P. A. Iglesias, and P. N. Devreotes, Proceedings of the National Academy of Sciences 107, 17079 (2010).
  * Mori  _et al._ [2008] Y. Mori, A. Jilkine, and L. Edelstein-Keshet, Biophysical journal 94, 3684 (2008).
  * Bergmann and Zimmermann [2019] F. Bergmann and W. Zimmermann, PloS one 14, e0218328 (2019).
  * Brauns  _et al._ [2020] F. Brauns, J. Halatek, and E. Frey, Physical Review X 10, 041036 (2020).
  * Halatek and Frey [2018] J. Halatek and E. Frey, Nature Physics 14, 507 (2018).
  * Alonso and Baer [2010] S. Alonso and M. Baer, Physical biology 7, 046012 (2010).
  * Trong  _et al._ [2014] P. K. Trong, E. M. Nicola, N. W. Goehring, K. V. Kumar, and S. W. Grill, New Journal of Physics 16, 065009 (2014).
  * Bosgraaf and Van Haastert [2009] L. Bosgraaf and P. J. Van Haastert, PloS one 4, e5253 (2009).
  * Kae  _et al._ [2004] H. Kae, C. J. Lim, G. B. Spiegelman, and G. Weeks, EMBO reports 5, 602 (2004).
  * Etzrodt  _et al._ [2006] M. Etzrodt, H. C. Ishikawa, J. Dalous, A. Müller-Taubenberger, T. Bretschneider, and G. Gerisch, FEBS letters 580, 6707 (2006).
  * Sasaki  _et al._ [2007] A. T. Sasaki, C. Janetopoulos, S. Lee, P. G. Charest, K. Takeda, L. W. Sundheimer, R. Meili, P. N. Devreotes, and R. A. Firtel, Journal of Cell Biology 178, 185 (2007).
  * van Haastert  _et al._ [2017] P. J. van Haastert, I. Keizer-Gunnink, and A. Kortholt, Molecular biology of the cell 28, 922 (2017).
  * Shi  _et al._ [2013] C. Shi, C.-H. Huang, P. N. Devreotes, and P. A. Iglesias, PLoS computational biology 9, e1003122 (2013).
  * Huang  _et al._ [2013] C.-H. Huang, M. Tang, C. Shi, P. A. Iglesias, and P. N. Devreotes, Nature cell biology 15, 1307 (2013).
  * Hecht  _et al._ [2010] I. Hecht, D. A. Kessler, and H. Levine, Physical review letters 104, 158301 (2010).
  * Hecht  _et al._ [2011] I. Hecht, M. L. Skoge, P. G. Charest, E. Ben-Jacob, R. A. Firtel, W. F. Loomis, H. Levine, and W.-J. Rappel, PLoS Comput Biol 7, e1002044 (2011).
  * Neilson  _et al._ [2011] M. P. Neilson, D. M. Veltman, P. J. van Haastert, S. D. Webb, J. A. Mackenzie, and R. H. Insall, PLoS Biol 9, e1000618 (2011).
  * Otsuji  _et al._ [2007] M. Otsuji, S. Ishihara, K. Kaibuchi, A. Mochizuki, S. Kuroda, _et al._ , PLoS Comput Biol 3, e108 (2007).
  * Camley  _et al._ [2013] B. A. Camley, Y. Zhao, B. Li, H. Levine, and W.-J. Rappel, Physical review letters 111, 158102 (2013).
  * Shao  _et al._ [2012] D. Shao, H. Levine, and W.-J. Rappel, Proceedings of the National Academy of Sciences 109, 6851 (2012).
  * Camley  _et al._ [2017] B. A. Camley, Y. Zhao, B. Li, H. Levine, and W.-J. Rappel, Physical Review E 95, 012401 (2017).
  * Nickaeen  _et al._ [2017] M. Nickaeen, I. L. Novak, S. Pulford, A. Rumack, J. Brandon, B. M. Slepchenko, and A. Mogilner, PLoS computational biology 13, e1005862 (2017).
  * Alonso  _et al._ [2018] S. Alonso, M. Stange, and C. Beta, PloS one 13, e0201977 (2018).
  * Stankevicins  _et al._ [2020] L. Stankevicins, N. Ecker, E. Terriac, P. Maiuri, R. Schoppmeyer, P. Vargas, A.-M. Lennon-Duménil, M. Piel, B. Qu, M. Hoth, _et al._ , Proceedings of the National Academy of Sciences 117, 826 (2020).
  * Ecker and Kruse [2021] N. Ecker and K. Kruse, Plos one 16, e0246311 (2021).
  * Miao  _et al._ [2017] Y. Miao, S. Bhattacharya, M. Edwards, H. Cai, T. Inoue, P. A. Iglesias, and P. N. Devreotes, Nature cell biology 19, 329 (2017).
  * Moreno  _et al._ [2020] E. Moreno, S. Flemming, F. Font, M. Holschneider, C. Beta, and S. Alonso, Physica D: Nonlinear Phenomena , 132591 (2020).
  * Cao  _et al._ [2019a] Y. Cao, E. Ghabache, Y. Miao, C. Niman, H. Hakozaki, S. L. Reck-Peterson, P. N. Devreotes, and W.-J. Rappel, Journal of the Royal Society Interface 16, 20190619 (2019a).
  * Cao  _et al._ [2019b] Y. Cao, E. Ghabache, and W.-J. Rappel, Elife 8, e48478 (2019b).
  * Moure and Gomez [2017] A. Moure and H. Gomez, Computer Methods in Applied Mechanics and Engineering 320, 162 (2017).
  * Flemming  _et al._ [2020] S. Flemming, F. Font, S. Alonso, and C. Beta, Proceedings of the National Academy of Sciences 117, 6330 (2020).
  * Löber  _et al._ [2015] J. Löber, F. Ziebert, and I. S. Aranson, Scientific reports 5, 9172 (2015).
  * Kulawiak  _et al._ [2016] D. A. Kulawiak, B. A. Camley, and W.-J. Rappel, PLoS computational biology 12, e1005239 (2016).
  * Moreno  _et al._ [2021] E. Moreno, R. Grossmann, C. Beta, and S. Alonso, Preprint (2021).
  * Kitano [2007] H. Kitano, Towards a theory of biological robustness (2007).
  * Imoto  _et al._ [2020] D. Imoto, N. Saito, A. Nakajima, G. Honda, M. Ishida, T. Sugita, S. Ishihara, K. Katagiri, C. Okimura, Y. Iwadate, _et al._ , bioRxiv (2020).
  * Taniguchi  _et al._ [2013] D. Taniguchi, S. Ishihara, T. Oonuki, M. Honda-Kitahara, K. Kaneko, and S. Sawai, Proceedings of the National Academy of Sciences 110, 5016 (2013).
  * Cooper  _et al._ [2012] R. M. Cooper, N. S. Wingreen, and E. C. Cox, PloS one 7, e33528 (2012).
  * van Haastert [2021] P. J. van Haastert, PloS one 16, e0246345 (2021).
  * García-Ojalvo and Sancho [2012] J. García-Ojalvo and J. Sancho,  _Noise in spatially extended systems_ (Springer Science & Business Media, 2012).
  * Hanggi  _et al._ [1985] P. Hanggi, T. J. Mroczkowski, F. Moss, and P. V. McClintock, Physical Review A 32, 695 (1985).
  * Kramers [1940] H. A. Kramers, Physica 7, 284 (1940).