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 at the front triggering push of the cytoskeleton and the motion of the amoeba.

We investigate a simple reaction-diffusion model [23] for Ras-GTP (R) patches that consists on the next partial differential equation: ∂ R ∂ t = 1 − R k 0 + k 1 S + k 2 R n 1 R n 1 + K R n 1 − k 3 G R − k 4 L R − k 5 R 1 + k 6 P + D R ∇ 2 R + ξ R x , t , (1) where we have a basal production, the stimulation production R by local occupied receptor S, an autocatalytic stimulation of R, an inhibition R by global and local inhibitors (G _R and L _R ), the degradation of R, the diffusion of R and, a Gaussian spatio-temporal distributed white noise ξ _R (x, t) with zero mean ⟨ξ _R (x, t)⟩ = 0 and correlation ⟨ξ _R (x, t)ξ _R (x′ , t′)⟩ = 2σ _R δ(x − x′ )δ(t − t′). The external S can significantly increase the response of Ras [23], see diagram in Figure 1, therefore, we keep here S = 0.

The variables L _R and G _R correspond, respectively, to local and global inhibitors of R: ∂ G R ∂ t = k 7 R ̄ − k 9 G R , (2) ∂ L R ∂ t = 1 − L R k 10 R − k 11 L R + D L R ∇ 2 L R ; (3) where both are produced by R and degrade, and only the local inhibitor diffuses. Note that the quantity R ̄ corresponds to the spatial average in the whole space of the field R.

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 Figure 1. This variable P was coupled with its respective inhibitors: ∂ P ∂ t = 1 − P k 12 + k 13 R + k 14 P n 2 P n 2 + K P n 2 + k 15 M n 3 M n 3 + K M n 3 − k 16 G P − k 17 L P − k 18 P + D P ∇ 2 P + ξ P x , t , (4) where we have a basal production, the stimulation production by R, an autocatalytic stimulation of P, a production by an extra term M, an inhibition of P by global and local inhibitors (G _P and L _P ), the degradation of P, the diffusion of P and, a Gaussian spatio-temporal distributed white noise ξ _P (x, t) with zero mean ⟨ξ _P (x, t)⟩ = 0 and correlation ⟨ξ _P (x, t)ξ _P (x′ , t′)⟩ = 2σ _P δ(x − x′ )δ(t − t′).

Similar to the previous set of equations, the variables L _P and G _P correspond, respectively, to local and global inhibitors of P: ∂ G P ∂ t = k 19 P ̄ − k 20 G P , (5) ∂ L P ∂ t = 1 − L P k 21 P − k 22 L P + D L P ∇ 2 L P ; (6) where both are produced by P and degrade, and only the local inhibitor diffuses. Note that the quantity P ̄ corresponds to the spatial average in the whole space 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: ∂ M ∂ t = k 23 P − k 24 M + D M ∇ 2 M . (7) which is produced by P, degrades and diffuses. 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 [49, 50].

A more exhaustive description of the model is shown in the original study [23]. However, 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 [51]; for more details see Supplementary Material.

One dimensional simulations were made using periodic boundary conditions in a grid of 120 points. The cell was considered as circular and having a radius of 6.25 μm. The pixel size for this case was set at Δx = 0.32 μm and the time step Δt = 0.03 s. The definition and the value of the parameters of the model can be found in Supplementary Table S1.

The original model did not include deformable cells. Therefore, we expanded the model to 2D geometry and introduced an auxiliary phase field ϕ 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 ϕ = 1 inside and ϕ = 0 outside of the cell.

The phase field equation is the result of a force balance involving several types of forces of different origins acting in the cell body. The equation for the phase field is as follows τ ∂ ϕ ∂ t = γ ∇ 2 ϕ − G ′ ϕ ϵ 2 − β ∫ ϕ d A − A 0 ∇ ϕ + α ϕ P ∇ ϕ , (8)

The first term on the right side in Eq. 8 is related to the surface energy of the cell membrane, where γ is the surface tension and G(ϕ) = 18ϕ ²(1 − ϕ)² is a double well potential. The second term keeps the area close to the value of A ₀. And the third term represents the active force generated by the Rac-GTP (P) molecules when pushing on the cell membrane [34]. The inclusion of the phase field permits to mimic the shape of the ventral membrane of crawling cells.

We consider a deformable cell of 122 μm ² in area corresponding to a circle with radius equal to 6.25μm. The pixel size used was the half as in the 1D case Δx = Δy = 0.16 μm. Also, to increase accuracy, the time discretization was reduced to Δt = 0.003s. We integrated Eqs 1–8 using periodic boundary conditions and standard finite differences. The rest of parameters are as displayed in Supplementary Table S1.

The mathematical model represented in Eqs 1–7 has several conservation terms which affect the dynamics of the system; see for example Eqs.2–5. However, we further modify the original model to include a new feedback through parameters k ₁₄ to create a more robust model by the control of the bistable properties, which are discussed in following sections. Therefore, the parameter k ₁₄ is dynamically controlled depending on the total amount of the protein of the inducer P at the cell. A new control term related to k ₁₄ is incorporated, following previous similar approaches in simple models [32, 38, 47]: k 14 = k 14 ∗ + η P ̄ − C P , (9) where k 14 * is a new constant that take the values of the original value of k ₁₄. It is important to note that the new constrain precludes the cell to be fully covered with the component P or emptying of the component P. Both extreme conditions are rare during cell crawling.

Parameter C _P is the target value of the fraction of the cell area occupied by patches of pseudopod inducer P. When, the total concentration P ̄ is larger (lower) than C _p the parameter value of k ₁₄ is larger (lower) than k 14 * and the bistable conditions implies the reduction (increase) of the P ̄ . It corresponds to a feedback control of the parameter to keep a particular stable solution in Eq. 9. Note that for η = 0 we recover the original model.

In a one-dimensional system, mimicking the membrane of a crawling cell, the generation of a single 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 [23], can be numerically integrated into a one-dimensional domain as in the original study, where kymographs are shown for a certain window of parameter values.

If we vary the parameters k ₁₄ and k ₁₈, we obtain the phase diagram; see Figure 2. For large values of k ₁₄ and small values of k ₁₈ the membrane is completely covered by P, while for small values of k ₁₄ and large values of k ₁₈ the concentration of P strongly decreases. Locomotion is, therefore, expected for intermediate values of these two parameters as shown in Figure 2. For the evaluation of the mechanism of the simulations shown in Figure 2, we fixed the stochastic source of perturbations of the additive noises in Eqs 1–7 to σ _R = 0.04 and σ _P = 0.025.

For the evaluation of the mechanism of pattern formation we found the stationary solutions for the deterministic version of Eqs 1–7 by removing the additive noise. We have slowly increased the parameter k ₁₄ by small intervals (0.05 s ⁻¹) and integrated the model for a considerable amount of time to permit the model reach the equilibrium (180 s). Subsequent increase of the parameter permits the observation of the evolution of the stationary value of P by obtaining a value for P in every interval in k ₁₄; see Figure 3A. Once the system of equations saturated to a certain value of P ∼ 0.75 − 0.80 we stopped and reduced the value of k ₁₄, giving rise to a hysteresis cycle; see Figure 3A. The hysteresis cycle shows that the system actually is bistable for certain window of parameter values. Equivalent dynamics of P can be observed under a similar change in k ₁₈; see Figure 3B. Therefore for low and large values of k ₁₄ and k ₁₈ there is a single stable solution and for intermediate values the model shows bistability.

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 [52, 53]. This formation is the mechanism responsible for the formation of localized domains of P. The localized pattern observed in Figure 2 is not due to excitable dynamics but rather to a delicate equilibrium between bistable dynamics and noise.

For the coupling of the polarization mechanism to the cell motion 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 membrane of a crawling cell.

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 Figure 4 we showed snapshots of the in silico cells with the parameter values corresponding to the values shown in Figure 2. We reproduced the dynamics expected from the one dimensional simulations, and the crawling dynamics can be clearly studied. Again, for high values of parameter k ₁₄ and small values of k ₁₈ the concentration of P is mostly high and homogeneously distributed along the membrane; see Figure 4. On the other hand, in the opposite limit, for small values of k ₁₄ and high values of k ₁₈ the concentration of P 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 Figure 4 we obtain fan-shaped cells, a typical mode of Dictyostelium discoideum cell crawling [37].

The dependence of locomotion on the parameters k ₁₄ and k ₁₈ is strong as shown in Figure 4. Persistent motion is observed only for an intermediate window of values of the two parameters. For very small or very large amounts of P ̄ cells do not move. Such quantity is determined by the parameters k ₁₄ and k ₁₈, and small changes on their values prevent the cell motion, making the whole model non-robust.

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 P needed to give rise to cell movement; see Figure 5. 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 Figure 5. This type of motion is reminiscent of the fan-shaped amoeboid cells previously reported [37] and also reproduced with simpler models [38]. 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 Figure 5 for higher noise intensities. 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 Figure 4, small variations in the parameter values k ₁₄ and k ₁₈ may completely change the movement dynamics.

In order to reduce the high sensitivity of the appearance of motion on the parameter values, we included a conservation constraint on the total inducer P, see Eq. 9, in the biochemical rates responsible for the bistable transition in Eq. 4. Following similar previous approaches [34, 47], we included this conservation as a global feedback condition in the reaction-diffusion equations of the model, see Section 2.3.

In Figure 6, we show several realizations incorporating the conservation constraint of P proteins, for different k ₁₄ and k ₁₈ parameter. A direct comparison can be made with the results displayed in Figure 4 showing the effects of the inclusion of the mass conservation feedback. Results in Figure 6 show that for a wide range of parameter values, the variation of the parameters does not affect either the bistability of the model or the quantity of inducer P 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.

The shape and type of motion were affected by the noise intensity. Using the same parameters of Figure 5 we studied the effects of varying the noise variance in the new mass-conserved model. When the noise intensity related to the pseudopod inducer P increases, changes the shape from persistent fan-shape to random amoeboid phenotype, see Figure 7. This change is present because while small values of σ _P generates less blurred patches, an increment of σ _P drives the opposite effect, the appearance of more blurred and distributed patches for P.

We measured 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 Figure 8.

For the original model the speed and shape drastically changes when varying the parameters, and only a particular combination of parameters gives rise to cell crawling (see Figure 4). On the other hand, 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 Figure 8. Such results were obtained by changing k ₁₈ to correspond to simulations shown in the third row of Figure 4 for the no mass conservation case and for Figure 6 for the mass-conservation system.

In Figure 9 we show the dependency of the speed on the noise σ _P for different values of σ _R . The speed of the cells employed for the calculation comes from the dynamics shown in Figures 5, 7. We observe that the speed obtained from both conditions decreases for greater noise intensities. However, for the original model almost null velocities appear for low noise magnitude. As we increase the noise intensity, values of the speed are more similar.