Waves in a Stochastic Cell Motility Model

In Bhattacharya et al. (Science Advances, 2020), a set of chemical reactions involved in the dynamics of actin waves in cells was studied. Both at the microscopic level, where the individual chemical reactions are directly modelled using Gillespie-type algorithms, and on a macroscopic level where a deterministic reaction-diffusion equation arises as the large-scale limit of the underlying chemical reactions. In this work, we derive, and subsequently study, the related mesoscopic stochastic reaction-diffusion system, or Chemical Langevin Equation, that arises from the same set of chemical reactions. We explain how the stochastic patterns that arise from this equation can be used to understand the experimentally observed dynamics from Bhattacharya et al. In particular, we argue that the mesoscopic stochastic model better captures the microscopic behaviour than the deterministic reaction-diffusion equation, while being more amenable for mathematical analysis and numerical simulations than the microscopic model.


Introduction
In order to move around, an amoeboid cell can change its shape by polymerising actin to curve the cell membrane. The actin polymerisation is controlled by signalling molecules and experiments in Dictyostelium discoideum have shown that activation of these signalling molecules happens at localised patches that can move along the membrane like a wave [1,21]. In wild-type (WT) cells, these waves move fast and die out, creating familiar-shaped pseudopods, while in cancerous cells these waves stick to a point, creating elongated protrusions [1], see Figure 1.1. In absence of a signal, the formation of pseudopods happens at random places on the cell membrane, resulting in random motion. In contrast, when a cell senses a chemical signal, it can concentrate the random protrusions at the side of the cell where the signal comes from, leading to movement in the direction of the signal [6]. As cells are small, the difference in signal strength between the front and the back of the cell (the gradient) is small as well. Furthermore, the cell can only use discrete points at the membrane where the receptors are to estimate the direction of the signal [6]. Therefore, one of the main questions is "How can a cell use a small gradient in the signal to concentrate the actin activity in the front?". This question has been studied intensively, but no complete description of all the microscopic chemical processes involved has been given yet, see [8] for a review.
In [1], the choice is made to describe the highly complex actin dynamics with a conceptual activator u and inhibitor v that diffuse and react with each other as summarised in Table 1. The species u and v are an abstraction of the dozens of components that regulate the actual cell movement, but the activator u can be thought of as Ras activity [1], which plays an important role in cell growth and differentiation [28]. In particular, u is being activated by Reaction #3 and Reaction #4, while being inhibited by Reaction #1 and Reaction #2, with propensities as indicated in the table. In addition, v is inhibited by Reaction #5, while Reaction #6 activates the inhibitor.
The information on the chemical reactions, in combination with the diffusion of both species, is generally used in one of two ways. First, there is a Gillespie-type algorithm [15,16] which can be used to simulate the involved chemical reactions on a microscopic level. For these simulations, (u k (t n ), v k (t n )) (the solution at time t n at grid cell k) is treated as the number of molecules of type u and v at time t n in a grid cell with finite size. For all these individual molecules the probabilities of diffusing to other grid cells or taking part in a chemical reaction are prescribed as by Table 1. To be precise, Reaction #1 implies that the time to the next reaction that degrades a u molecule in grid cell k is exponentially distributed with rate parameter (a 1 u k (t n )) −1 . See the panels on the left of Figure 1.1 for examples of these simulations. This Gillespie-type algorithm approach takes the stochastic nature of a single cell into account. However, it is computationally very expensive and difficult to analyse mathematically. Hence, it is hard to use this type of modelling approach to make valuable predictions.
A second way to use the reactions in Table 1 is to derive an average large-scale limit macroscopic equation. Hence, we assume that u and v are densities on a continuous domain, described by a reaction-rate equation with diffusion, also known as a Reaction-Diffusion Equation (RDE). In particular, the RDE 1 related to the chemical reactions in Table 1 is given by which is a specific version of the general RDE we will encounter in §2. This model is a variation on the classic FitzHugh-Nagumo model for neuron spiking [12,30]. Protrusions are formed at places with high activator u and u is inhibited by the terms −a 1 u and −a 2 uv, see Reaction #1 and Reaction #2 in Table 1. This implies that an increase in u or v leads to a decrease in u, unless the increase is high enough such that activation from Reaction #3, modelled by a nonlinear Hill function a 3 u 2 /(a 4 +u 2 ), takes over and negates the inhibiting effects. Effectively, this means that a small increase in u can lead to a much larger increase in u, that is, the system is locally activated. Once u is large and the Hill function levels off at a fixed value a 3 , the amount of inhibitor v increases via the term εc 2 u (related to Reaction #6), leading to a fast decay in u by the −a 2 uv term (related to Reaction #2). The inhibitor v then decays via Reaction #5 to the rest state and activation can happen again. In addition, both species diffuse with diffusion coefficient D u , respectively D v , where it is assumed that D u < D v . It is important to realise that, in both approaches, the modelled actin waves happen on the surface of the cell, and, as in [1], we only study a slice of this surface. Therefore, the spatial domain must be thought of as an (approximate) circle. For deterministic RDEs like (1.1), a plethora of analytical tools are available (see, for instance, Appendix B) and numerical simulations are relatively straightforward. However, being a deterministic equation, this RDE does not show the same stochastic dynamics as the Gillespie simulations and experiments. A crucial difference between the macroscopic RDE model (1.1) and the Gillespie simulations revolves around the duration of the patterns. In the RDE, an established pattern, e.g. a standing or travelling wave, will, if uninterrupted, remain there for a very long time, while these patterns are destroyed quickly both in stochastic simulations and experiments. Furthermore, when the rest state of the RDE (1.1) is stable, activation cannot come from the RDE itself, but it needs an external signal large enough to activate the nonlinear term a 3 u 2 /(a 4 + u 2 ) related to Reaction #3. We generally refer to the activation of these patterns as activation events.
It is important to realise that the dynamics of the different chemical processes in the cell are inherently stochastic and at the size of a single cell chemical reactions are not well approximated by large-scale approximations, as Figures 1.1 and 1.2 show. In other words, treating the relevant enzymes and receptors like a continuous medium of infinitely many, infinitely small, particles is invalid, and the stochastic nature of reactions between individual molecules becomes important. This so-called internal noise can serve as a signal to activate the dynamics if it is large enough at a certain point in space and time. As we noted before, the cell hence executes a random walk in the absence of a signal 2 . This implies that an external signal does not necessarily activate the dynamics at a certain point on the membrane, but rather changes the random walk of the cell into a biased random walk in the direction of the signal. Using a more extended model than presented here, it is shown in [2] that coupling an external signal to the stochastic dynamics of the cell indeed can lead to movement in the direction of that signal.
Instead of studying the complex internal dynamics of the cell, it can be advantageous to perturb the deterministic RDE (1.1). For instance, in [1], an external source of noise is applied to the RDE (1.1), turning it into a Stochastic RDE (or Stochastic Partial Differential Equation (SPDE)). While this approach can indeed activate the dynamics and make long-term deterministic waves collapse, it is inherently ad hoc and not a priori based on any of the involved biologically relevant processes.
In between the macroscopic level of the RDE and the microscopic level of the chemical reactions, one can derive a mesoscopic SPDE, known as a Chemical Langevin Equation (CLE) [18], that also incorporates the internal noise of the cell. In §2, we will show that the SPDE associated with the chemical reactions as described in Table 1 plus diffusion is given by Here, (dW 1 t , dW 2 t ) and (dW 1 t , dW 2 t ) are two independent noise vectors with space-time white noise (each component is also independent of the other) and σ is a measure for the strength of the noise. Indeed, in the no-noise limit σ → 0 the mesoscopic SPDE (1.2) reduces to the macroscopic RDE (1.1). In that sense, σ serves as a scale parameter.
The main advantage of the SPDE description is, on one hand, that the solutions still show the rich dynamics of the Gillespie models, i.e. the activation and destruction of waves, but are computationally significantly less expensive. On the other hand, since the SPDE in the no-noise limit reduces to the deterministic RDE model (1.1), we can use well-developed Partial Differential Equation (PDE) theory to gain insight into the dynamics of the RDE (1.1) and use this to study the closely related SPDE, see for instance [19,25]. To give an idea of the differences between the deterministic and stochastic models we plot two simulations in Figure 1.2 that will be discussed later in §3. It is clear that the simulation of the SPDE paints a much more dynamic picture than the deterministic one, which is more in line with the inherently noisy nature of the cell's chemical processes. Hence, SPDEs are an invaluable tool in unravelling the dynamics of a cell.
This article is now organised as follows. In §2 we explain how to derive the SPDE (1.2) from Table 1. Subsequently, in §3 we study both the SPDE (1.2) and the RDE (1.1) numerically in different parameter regimes and qualitatively compare the observed dynamics to the Gillespie simulations from [1]. In §4, we discuss the results and how they relate to the questions posed in this introduction.

Derivation of the SPDE
Our starting point to derive (1.2) is the set of chemical reactions as laid out in Table 1. First, we introduce the column vector X(t) = (u(t), v(t))) T , where T indicates that we transpose the row vector, and the column vector R(X(t)) with the propensities of the six reactions: The associated stoichiometric matrix S, which describes the change in X(t) for each reaction, is then given by see the last two columns of Table 1. On top of these reactions, we assume that both variables also diffuse, so for a well-mixed solution in a large container we find the classic PDE where D is a diagonal diffusion matrix with coefficients D u and D v on the diagonal [3]. This PDE is identical to the RDE (1.1) and describes the dynamics of X(t), averaged over many individual reactions. When the number of reacting molecules is large enough, and when we zoom out far enough such that all individual molecules become effectively a density, the macroscopic PDE gives a good approximation of the microscopic behaviour. Statistically speaking, this means that the probability distribution of all possible states must be very sharply peaked around the average value described by the PDE, so the deviations from the mean can be ignored.

Motivating Example
The assumption that we can ignore deviations from the mean is not always valid. For example, in population dynamics, we can write down birth-death models for several hundred individuals and with this number of individuals, random deviations from the mean are actually significant. To further exemplify, and to set the stage for the upcoming derivation, let us study such a simple discrete birth-death process: suppose a population is at time t in state X(t). In the next timestep dt, there are three possible outcomes: (i) the population grows by one individual with probability b(X(t))dt, (ii) the population decreases by one individual with probability d(X(t))dt, or (iii) nothing happens to the population with probability 1 − b(X(t))dt − d(X(t))dt. Now, assume we have a continuous Stochastic Differential Equation (SDE) where β t is Brownian motion, i.e. we can think of dβ t as a random step with average zero and variance dt. We now ask the question: "When is this continuous SDE a good approximation of the described discrete birth-death process?". Or, more precisely, "What should f (x) and g(x) be such that (2.3) is a good approximation of the described discrete process?". Given a solution x of the SDE, we see that the average expected value at x(t + dt) is approximated, at lowest order in dt, by For the described birth-death process, we have that the expectation is Hence, the average expected jump size in population is identical for the SDE (2.3) and the birthdeath process if we take f ( Next, we compute the deviation from the mean of the SDE (2.3) while this deviation for the birth-death process is Therefore, to make these deviations coincide at first order in dt, we must take g( Hence, the process x(t) described in (2.3), which is continuous in population size and time, is a good approximation of the discrete process X(t) when The stochastic process x(t) shares the average and variance with X(t) but differs in other points. Higher order moments of x(t) and X(t) will not be identical and x(t) can become negative, even when b and d are chosen such that this is not possible in the discrete model. In order to link the SDE above to chemical reactions, we make the following observation. The birth of an individual can be thought of as the chemical reaction ∅ → X with propensity b(X) and stoichiometric value 1, while the death of an individual can be seen as the chemical reaction X → ∅ with propensity d(X) and stoichiometric value −1. Next, we make an assumption which is called the leap condition [3]. That is, we assume that, given a state X(t), enough reactions happen in the interval [t, t+dt] to describe the average jump size in [t, t+dt] by a Poisson process whose parameters depend on X(t). With this leap condition assumption, we implicitly also assume that X(t) is a good approximation of the solution in the whole time interval [t, t + dt]. We now turn the discrete process X(t) into a continuous process x(t) by approximating the discrete Poisson process by a continuous Gaussian, see [24] for details. This approach results in an SDE similar to the SDE (2.4): for two independent Brownian motions β 1 t and β 2 t . Although visually different from (2.4), both SDEs have a noise term that is Gaussian with identical average and variance. Therefore, both SDEs describe the same stochastic process and hence we can say that (2.4) and (2.5) are equivalent.

Derivation of the CLE
We have now gained some intuition for linking more general discrete chemical reactions to continuous S(P)DEs: if we have M different molecules in a vector X(t) with diffusion matrix D, N reactions given by a vector R(X(t)) and a stoichiometric matrix S, then the continuous SPDE for X(t) is given by see [3,24]. The equation is made of two parts, a local equation that describes the kinetics as in SDE (2.4) and a stochastic diffusion equation as derived in [10]. Here, dW t and dW t are two independent vectors with space-time white noise. The vector dW t has N components coming from the N reactions, while dW t has the dimension M of X(t). SPDE (2.6) is known as the Chemical Langevin Equation (CLE) [18]. The vector X(t) now describes the densities of the molecules involved, not the actual number of molecules. How well the discrete number of molecules is approximated by a density is determined by the scale parameter Ω and is in that sense a measure for the noisiness of the system. In the no-noise limit Ω → ∞, we recover the classic RDE (2.2). In contrast, for small Ω the dynamics of the discrete process is dominated by random events and the discrete process should be described in full detail by a chemical master equation [17]. The CLE can be understood as the lowest order approximation of the chemical master equation for large Ω, see for more details [3]. For an overview of all different paths leading from molecular kinetics to (S)PDEs, see [26,Fig. 3.4]. It is important to realise that SPDE (2.6) does not necessarily inherit all the statistical properties of the chemical master equation, only averages and variances. Another potential issue is that it does not necessarily ensures positivity of the solutions. Just as (2.4) and (2.5) are identical, we can rewrite (2.6) in the following way: This time, the noise vector dW t has just M components, reducing the number of random vectors that must be generated (when M < N ). The downside is that the computation of Sdiag(R(X))S T is in general numerically more expensive then the computation of S diag(R(X)). However, in the present setting, there are no connections between the two variables in the stoichiometric matrix S (2.1) and the matrix Sdiag(R(X))S T is thus diagonal, making the computation of the square root trivial.
Note that once we have the CLE (2.9), it can be applied to any set of chemical reactions and can therefore have widespread use. For example, we can now return to Table 1 and apply the CLE to these reactions, which results in (2.10) For notational convenience, we replaced 1/ √ Ω by a small parameter σ, resulting in the SPDE (1.2) from the introduction. In the remainder of this work, we will study the SPDE above, mainly using numerical techniques.

Remark 1.
It is important to realise that the SPDE above does not have a function-valued solution in general. The term ∂ x 2DX(t)dW t makes the equation ill-posed and solutions can only be understood in terms of distributions. Therefore, it is not a priori clear if the numerical solutions shown in the next section converge to a solution of the SPDE when the spatio-temporal discretisations dx and dt are sent to 0. In §3.4, we will discuss the implications of omitting this term on the wave dynamics.

Simulations
In this section, we will numerically investigate the PDE (1.1) and SPDE (2.10). We investigate three of the main building blocks of the PDE dynamics: localised standing waves, localised travelling waves and time-periodic solutions, together with their counterparts in the SPDE. However, before we can investigate the dynamics, we must first establish some basic properties of the (S)PDE, like the existence, uniqueness and stability of the background state(s).
Since we are interested in localised waves and expect the activator to be in rest otherwise, we need for the existence of these localised waves that the spatially homogeneous background state is stable. In contrast, for the time-periodic solutions, we expect the background state to be unstable such that continuous excitations of the background state can happen. The possible background states (u * , v * ) of (1.1) are given by the positive real solutions of the u-nullcline and v-nullcline See Figure 3.1a for a typical representation of the shape of the nullclines. Since the system parameters are all assumed to be positive, this is equivalent to finding the positive solutions u * of with v * = c 2 u * /c 1 . Due to the complexity of the general solution formula for quartic polynomials, it is not feasible to write down its solutions explicitly. However, by Descartes' rule of signs [7] we know that there is only one positive real root if c 1 (a 3 + a 5 ) < a 2 a 4 c 2 and one or three positive real roots otherwise 3 . The stability of a background state (u * , v * ) is then determined by the eigenvalues of the associated Jacobian matrix Since we do not have an explicit formula for (u * , v * ), we must compute these eigenvalues numerically.
For example, when we allow one free parameter, e.g. c 1 , and fix the other values, then we can compute the background states and the associated eigenvalues of the Jacobian matrix. Taking the parameter values a 1 = 0.167, a 2 = 16.67, a 3 = 167, a 4 = 1.44, a 5 = 1.47, ε = 0.52 and c 2 = 3.9 from [1] and letting c 1 range from 0.18 to 0.35, such that c 1 (a 3 + a 5 ) < a 2 a 4 c 2 , results in one admissible positive background state ranging from (u * , v * ) ≈ (0.077, 1.669) to (u * , v * ) ≈ (0.142, 1.586). Initially, for the lower values of c 1 , the eigenvalues are real and negative, resulting in a stable background state. Increasing the value of c 1 to approximately 0.25 results in complex eigenvalues, still with negative real parts. When we further increase the value of c 1 to approximately 0.29, both eigenvalues cross the imaginary axis, i.e. the background state undergoes a Hopf bifurcation and we expect to see timeperiodic solutions. See Figure 3.1b for a visual representation of the evolution of the eigenvalues. In Figure 3.1a we show the nullclines for c 1 = 0.18 and c 1 = 0.35. The unique background state moved along the fold in the u-nullcline and as long as the background state is in between the two folds, the fixed point is unstable.
In the next sections, we will study localised standing and travelling waves for the same parameter set with c 1 < 0.25 and for time-periodic solutions with c 1 > 0.29. The complex dynamics of pulse adding for c 1 -values in the intermediate regime between these two boundary values, where the eigenvalues of the Jacobian are stable but complex-valued, is outside the scope of this work, see for example [4] for more information.
So far, we only looked at background states, which are spatially homogeneous. However, we are interested in spatially nonhomogeneous patterns. By definition, a localised wave is a fixed profile (Φ u , Φ v ) that moves with a fixed speed c (possibly zero). Therefore, when we change the spatial coordinate x to ξ = x − ct using the chain rule, the profile (Φ u , Φ v ) is a stationary solution of the following shifted Ordinary Differential Equation (ODE): This ODE problem can be solved using numerical fixed-point algorithms. For these algorithms, a crude starting point is needed for the profile and the value of c, which can come from a PDE simulation. Note that this problem is translation invariant, meaning that we find a one-dimensional family of travelling waves, all shifted versions of each other. Hence, for the solver to converge, an extra condition to fix the location of the wave is necessary.

Standing Waves
In this section, we will study standing waves, which means we look for solutions of (3.3) with c = 0. A solution to this ODE is shown in Figure 3.2a. We observe that both components u and v indeed start at and return to their background state (u * , v * ) ≈ (0.0523, 2.0394). We observe that the activator u changes rapidly in a small region in the spatial domain and we, therefore, call the activator u the fast variable. On the other hand, the inhibitor v is the slow variable as it changes more gradually over a larger spatial distance. Figure 3.2b shows the corresponding phase plane. The majority of the spatial dynamics happens near the lower branch of the u-nullcline before it has a fast excursion from the lower branch to the upper branch of this nullcline and, by the symmetry x → −x of the ODE (3.3), it then returns back to the lower branch in a similar fashion. The fact that both components of the standing pulse evolve on a different spatial scale allows us to mathematically analyse this standing pulse, see Appendix B. For instance, the valuev at which the activator u makes a sharp transition (approximately 3.8 in Figure 3.2b), can be approximated by the algebraic relation (B.10). The analysis also explains why the solution trajectory in the phase plane closely follows the lower branch of the u-nullcline for the most part of the trajectory. By assumption, the standing wave in Figure 3.2a is a stationary solution of the PDE (1.1). This can be confirmed by using the wave from the ODE as the initial condition for a PDE simulation (not shown). However, we are not likely to find this single standing wave in a PDE simulation  without a fine-tuned initial condition. As an example, we use for the simulation the initial condition u 0 = u * + e −x 2 and v 0 = v * + 2/ cosh 2 (5x) as a crude approximation of the wave. The resulting simulation is shown in Figure 3.3. This initial condition splits in, what appears to be, two wellseparated localised standing waves 4 . However, the plot of the slow v-component makes clear that this is not the case, and that the two standing waves are connected through the slow component, i.e. the slow component is not in its rest state in between the two standing waves. For more details on the numerics of the (S)PDE simulations, see Appendix A.
The interaction between the two standing waves in Figure 3.3 through the slow v-component makes that the two standing waves repel each other on a very slow timescale as is made clear by taking long integration times, see Figure 3.4b. On an infinite domain, the two standing waves slowly drift apart forever, but on a periodic domain, we can expect them to stabilise once they are at an equal distance on both sides. On the timescales of biological processes, this slow continuous splitting is probably not relevant and on short timescales, the term 'standing waves' for the solution at later times in Figure 3.3 is biologically justifiable. Furthermore, note that for our understanding of the presented dynamics, it is essential to look at both components simultaneously. In other words, for our understanding of Figure 3.3a it is essential to also look at Figure 3.3b.
We now take a closer look at the short-time dynamics presented in Figure 3.4a. In [1], this splitting of the initial condition is described as two counter-propagating travelling waves, sometimes called trigger waves [14]. By the formal mathematical definition, a travelling wave is a fixed profile moving with a fixed speed, i.e. a solution of (3.3). Therefore, mathematically speaking, these do not classify as travelling waves. Instead, what we observe here would be classified as transient dynamics and pulse splitting. However, it is clear that at t = 0, the activity of u is around x = 0, and after some time it moved to two different places, justifying the term 'travelling'. If we adopt the terms 'standing' and 'travelling', it is clear from Figure 3.3a that around t = 3 a transition occurs from travelling to standing.
Standing waves with noise For the same parameter values as in the previous paragraph, we now study the full SPDE (2.10). In Figure 3.5, we plot realisations of the SPDE for different noise intensities. For low noise levels, we see two quasi-stationary waves appear, like in Figure 3.3, before they are destroyed at different points in time by the noise. Since the noise is low, no new activation events happen. When we increase the noise intensity, the noise is able to activate the stable background state, but the waves are also destroyed more quickly, resulting in a constant appearance and disappearance of waves. Note the comparison between Figures 3.5c and the figures in [2], where a similar model is studied using Gillespie algorithms. This activation of the background state is not possible in the deterministic PDE (1.1) without an external force. In Figures 3.5b and 3.5c, we see that in the first instances, many patterns are generated, causing the inhibitor to increase everywhere which blocks new activation events. After this initial phase, new activation events appear, and significantly more for higher values of the noise as expected. When we increase the noise even further, it becomes impossible to form patterns as every activation event is destroyed instantly. Therefore, pattern formation happens at intermediate values of the noise. The idea that there is some 'optimal' value of the noise resulting in complex dynamics has been observed before in, for instance, the context of nerve impulses [13]. In order to quantify this notion of optimality in the noise intensity we must first quantify the size and shape of the patterns in Figures 3.5b and 3.5c. Using Matlab's regionprops algorithm we can automatically detect the patches with a high value for the activator u (see Appendix A for details), giving us the possibility to compute the number of activation events and determine the width and duration of each event, see Figure 3.6a. In Figure 3.6b we show the statistics for a range of σ values. This figure shows that there is a clear cutoff for when activation events are likely to happen. For values of σ < 0.035, the average number of events is lower than 1, and the number of activation events increases sharply after this value. We observe that the width, the length and the maximum height of the events are all higher when the number of excitation events is low, but the variability in these values is also larger. in Figure 3.7, we look at the statistics of the events for the specific value σ = 0.046. The value of the maximum is sharply peaked. This is something we expect, as the maximum is mainly determined by the deterministic dynamics after the excitation. The width and length of the events are much more spread out. Especially for the width, we see a heavy tail towards zero. This is also expected because activation events come in two forms. Most events result in two waves, but a small part of the events has the shape of just a single wave, which has a width of 0.87 in the deterministic case. We checked whether or not these histograms are well approximated by a  Figure 3.5, but with σ = 0.45 and a homogeneous initial condition with (u * , 4v * ). The red boxes are the result of the pattern finding algorithm regionprops in Matlab; it identifies all the regions of excitations which we would also find by eye, see Appendix A for details. In Figure (b), we used this algorithm to find the length, width and maximum of these pulses (left axis), as well as the total number of activation events (right axis). For each value of σ, the number of events is averaged over 100 simulations, and the length, width and maximum are averaged over all events in the 100 simulations. We plot the average together with the standard deviation.
Using the statistics on the width, length, and maximum, we can compare the solutions of SPDE (2.10) to SPDEs with the same deterministic part but different noise terms. First, we can set the ∂ x √ 2DXdW t term coming from the diffusion to zero. As noted in Remark 1, this term makes the mathematical analysis of the SPDE (2.10) significantly harder. Figures 3.7d-3.7f show that the statistics of the solutions do not change significantly when we delete this term. This indicates that the noise coming from the reaction terms plays a more influential role in determining the shape of the patterns.
We are now also in the position to compare the CLE approach with the more ad hoc approach of adding additive white noise the to u-component to mimic the inherent noisiness of the system, see

Travelling Waves
In order to find a travelling wave solution of (1.1), understood as a solution of (3.3) with c = 0, we must ensure that the dynamics starting from the initial condition does not reach the standing phase or returns to the background state. This can be achieved by increasing the value of c 1 . Increasing c 1 results in a faster exponential decay of v back to the background state after an excitation, see Table 1, preventing the inhibitor from glueing the two waves together like in Figure 3.3. Simulations for an increased value of c 1 , from 0.1 to 0.2 5 , are shown in Figure 3.8. Note that the PDE (1.1) still only has one stable background state (u * , v * ) ≈ (0.0833, 1.625). The initial condition splits into two counter-propagating travelling waves, but opposite to what happened with the standing wave before, they keep separating and move away from each other at a fixed speed until they collide and cancel each other out due to the periodicity of the domain, see Figure 3.8.
To find a single travelling wave, we again need to properly tune the initial condition. This can be done by selecting one of the two waves in Figure 3.8 and using it as the initial condition of the PDE simulation (not shown). In Figure 3.9 we show the travelling wave profile and its associated phase plane. As with the standing pulse, the dynamics around the u-nullcline is essential. The solution trajectory starts from near the background state and follows the lower branch of the u-nullcline, jumps towards the upper branch of the nullcline and keeps following it until it falls off and returns to the lower branch to slowly evolve back towards the stable background state. In contrast to the standing pulse, see Figure 3.2, the travelling wave is no longer symmetric and it jumps back to the lower branch by falling off the edge of the upper branch. These travelling wave solutions could be analysed further using techniques similar to Appendix B.
It is important to realise that we do not expect to see travelling waves in practice as the travelling wave gets destroyed when it collides with another wave. Therefore, in the stochastic simulations, it might not always be clear if we are looking at a travelling wave that collapses or at the transient dynamics towards a double pulse that subsequently gets destroyed by the noise.
Travelling waves with noise When we now return to SPDE (2.10), there are now four regimes for the same parameters as in the previous section. For high values of the noise, we, as before, do not observe any patterns (not shown). For low values of the noise, we just find the travelling wave (if the simulation is initiated by an appropriate initial condition) since the noise is not strong enough to destroy the wave, nor to activate another pattern, on the timescales of the simulation (not shown). The interesting dynamics happens again at the intermediate levels of the noise. As Figure 3.10a shows, the noise activates the dynamics, resulting in many counter-propagating travelling waves. A travelling wave is subsequently annihilated when it collides with a travelling wave coming from the other direction. Hence, the collision dynamics of Figure 3.8 is repeated many times on smaller spatial-temporal scales. We see in Figure 3.10 that after the annihilation of the travelling waves, the slow inhibitor v initially remains high preventing the activation of new counterpropagating travelling waves. Only when after a certain time the inhibitor has sufficiently decayed, do we see the activation of new counterpropagating travelling waves by the noise. The creation and annihilation of travelling waves happen at a shorter time scale than the decay of the inhibitor, which makes the dynamics look synchronised, or even periodic. In Figure A.2a we plot the approximate period versus the intensity of the noise. As expected, the period decreases with the intensity of the noise. It differs however significantly from the true time periodic motion we will discuss in §3.3. When we increase the noise, the quasi-periodic pattern is broken up, as the counter-propagating travelling waves are destroyed before they collide and annihilate each other, so no synchronised patterns emerge, see The red dashed line in (a) has a slope of 2.05, close to the deterministic wave speed, but given the short time interval the wave exists, precise estimates are difficult to obtain. We observe that there is a quasi-periodic behaviour with a period of roughly 20. In Figures (c) and (d), the quasi-periodic structure is destroyed. The same parameters and initial condition are used as in Figure 3.8.

Time Periodic Solutions
In the previous sections, it was essential that the background state of the system was stable, because this allowed the dynamics to return to the rest state after an activation event. When we increase the value of c 1 , the background state becomes unstable through a Hopf bifurcation, see Figure 3.1b. In the phase plane, this transition is characterised by the fact that the background state is no longer located on the lower branch of the u-nullcline, as in Figures 3.2b and 3.9b, instead, it lies on the middle branch of the u-nullcline, see Figure 3.12b. Hence, after an excursion, the solution cannot return to the unstable background state and is exited again, resulting in time-periodic motion. When we start with a spatial homogeneous initial condition, the PDE simulation shows periodic oscillations in time, see Figures 3.11 and 3.12. Both components still display slow-fast behaviour, however, this time not in the spatial variable x but in the temporal variable t. In the case of nonhomogeneous initial conditions, it takes several oscillations before they are all synced up spatially (not shown). The observed behaviour has the characteristics of a relaxation oscillation as studied intensively for the Van der Pol equation [32]. This is not a surprise as the Van der Pol equation formed the foundation for the classic FitzHugh-Nagumo model and PDE (1.1) can be seen as a variation on this classic model.

Time periodic solutions with noise
For small values of the noise σ, the observed period is close to the deterministic version, but when the value of σ increases, the period also decreases monotonically, as is expected. Note that after excitation, the inhibitor remains high preventing activation events. When the noise is too high no patterns are observed. We can investigate the relation between the reduction of the period and the intensity of the noise. In FigureA.2b, we plot the estimated period versus the noise intensity. We indeed see that the period decreases monotonically When we average over the x-direction and measure the distance between the maxima, we find T ≈ 7.87. Same parameters as in Figure 3.11 with σ = 0.01. with the noise.

Wild-Type versus PTEN-null Cells.
Now that we have studied several different fundamental patterns, we can focus on understanding the different cell shapes. In [1], two sets of parameters are compared, representing WT cells (i.e. healthy cells) and PTEN-null cells where the tumour-suppressing gene PTEN has been switched off [5]. First, we simulate the deterministic PDE (1.1) for both sets of parameters, see Figure 3.14.
We observe that in both parameter regimes, there are two counter-propagating travelling waves but the specific profiles and speeds are different. Especially, note that the wave in Figure 3.14b is significantly broader and higher than the wave in Figure 3.14a. When noise is applied, the statistics of the dynamics shows a clear difference. In Figure 3.15, we compare the SPDE simulations of (2.10) to the Gillespie simulations from [1]. Focusing on the typical shape of the excitations, there is a clear qualitative correspondence between the two types of simulation. Furthermore, in both types of simulation, the average pulse duration is longer in the case of the PTEN-null cell simulations. Note that we show the SPDE simulations on a larger spatiotemporal scale to get a better idea of the distribution in shapes and the zoom-boxes highlight the detailed structure of a typical single activation event. In the case of PTEN-null cells, the background state can be excited for much lower noise values (σ ≈ 0.007), while for WT cells, the noise needs to be twice as large (σ ≈ 0.014) as a result of the increased values of c 2 and a 3 . Hence, in PTEN-null cells, an already existing pattern can more easily sustain itself, leading to the elongated shapes of Figure 3.15d.

Discussion & Outlook
We set out to show how Stochastic Partial Differential Equations (CLE), or more specifically, Chemical Langevin Equations, can be used to gain more insight into the dynamics of models for cell motility. We have shown for an exemplary set of chemical reactions (see Tabel 1) that the CLE approach, combined with a basic analysis of the corresponding deterministic PDE, allows us to study the different possible patterns with relative ease, both qualitative and quantitative, while remaining close to the underlying chemical processes. To understand differences in cell behaviour, like the difference between wild-type and cancerous cells as in [1], the study of the statistical properties of the observed dynamics is essential. For instance, an essential characteristic differentiating wild-type cells from cancerous cells is how long a pattern can survive after activation. The simulations in the previous section show that the answer not only depends on the parameters of the system but crucially on the interplay between the parameters and the noise. The CLE can be used to study this interplay. A natural question to ask is if all the stochastic terms introduced in the CLE (2.9) are really necessary. Could we, for example, ignore the noise term coming from the diffusion or forget the derivation of the CLE altogether and just naively add an additive white noise term to the equation for u? The histograms in Figure 3.7 indicate that the effects of the terms that come from the diffusion are minimal (for the parameter values studied here) and therefore that these terms do not contribute meaningfully to our understanding of the cell dynamics. Note that this would solve the problem of the equation being ill-posed, see Remark 1, and would open up the possibilities for more rigorous mathematical analysis based on the results in [19]. We also noted that adding just additive white noise changes the statistics significantly, which indicates that completely abandoning the CLE approach throws away too much detail.
In this paper, we studied a basic activator-inhibitor system with only a limited number of chemical reactions. However, the derivation of the CLE (2.9) in §2 holds for any number of molecules and for any number of chemical reactions. As such, one can see this paper as a proof of concept and the methodology of this paper can be directly applied to more complex regulating systems, such as the eight-component system designed in [2]. In subsequent work, we aim to work on these type of more complex model to better understand the stochastic dynamics that causes the cell to move robustly in one specific direction. Furthermore, as shown in detail in Appendix B, the underlying deterministic RDE (1.1) is amenable for rigorous mathematical analysis by using Geometric Singular Perturbation Theory [11, 20, 22, 23, e.g.]. We derived a first-order approximation for the jump location where, under certain conditions, the standing wave has a sharp transition in its activator. This methodology could also be used to, for instance, further analyse the travelling waves to derive approximations for the speed of the waves. In other words, questions about the existence of localised solutions of (1.1) and bifurcations can thus be reduced to understanding relatively simple ODEs and the connections between them. The details of these computations are left as future work.

A.1 (S)PDE Simulations
All the (S)PDE simulations in this paper were done using a semi-implicit Euler-Maruyama method from [29, §10.5]. For the spatial directions, a standard 2nd-order central difference is used and for the time stepping Euler-Maruyama. The deterministic linear part is evaluated at the next timestep, making it semi-implicit. To be concrete, we study an SPDE of the form where u is a vector, L is a linear differential operator, f, g are functions and dW t a white noise vector. In comparison with the main text, the vector u equals (u, v) and L = D∂ xx .
When we denote the numerical approximation of the linear part L with A, and the spatial discretisation of u at time t with u(t), we find The white noise step dW t is a vector where each random element is distributed as N (0, dt/h). Hence, the approximation for the new value u(t + dt) becomes The equation for u(t + dt) is now a matrix equation and can be solved using standard solvers. It is important to realise that the algorithms from [29] only work for Lipschitz noise terms. Hence, when the term under the square root in (2.10) becomes close to zero, the algorithms become unstable. To correct this, we take after every timestep the maximum of u(t + dt) and 0. The specific models studied in the main text, even the PDE (1.1) can be very sensitive to the size of the spatial discretisation h and temporal discretisation dt in certain parameter regimes. For example, when c 1 = 0.15 and the remaining parameters are equal to those in For the values used in the main text, c 1 = 0.1 and c 1 = 0.2, such a discrepancy was not observed for reasonable discretisations. This is possibly related to the co-existence of travelling and standing waves in this regime.

A.2 Pattern Recognition
For Figures 3.6 and 3.7, Matlab's regionprops algorithm is used to identify the activation events automatically. This proceeds in the following steps. First, we smooth the data using Matlab's Gaussian filter. Without smoothing, the algorithm detects multiple objects in a single event. Next, we transform the data to a binary value by comparing it with a certain threshold: we say that u is activated when it is five times its stationary value u * . Then, the regionprops algorithm is applied with the option BoundingBox.
One needs to take care of which initial condition to use. When we start with a spatial homogeneous initial condition (u * , v * ), there is a lot of activation in the first instances of the simulation, see Figure 3.5c, and it is not possible to define and detect individual activation events. Therefore, we start not on the fixed point (u * , v * ), but on (u * , 4v * ) plus a small perturbation. The result is that activation events only appear when v has decayed enough for excitations to happen. As the decay is stochastic, and therefore not spatially homogeneous, the activation events start to appear more spread out, making it possible to determine individual events.

B Analysis
The deterministic PDE (1.1) has two components and ten parameters, making it difficult to directly analyse mathematically, even for the simplest of localised structures simulated in the main text. However, these simulations do reveal that the profiles of the two components of the PDE evolve on a different spatial scale: the spatial changes in the slow v-component are more gradual than these of the fast u-component, see, for instance, Figures 3.2a, 3.9a, and 3.12a. Furthermore, these simulations also revealed that a large part of the spatial dynamics centres around the lower and upper branch of the u-nullcline in the phase plane, with the u-profile making fast jumps in between. To amplify (and exploit) this scale separation, we set D v = 1 (as in [1]) and D u =ε 2 , whereε is a small parameter that can be taken arbitrary small. Furthermore, we assume that our spatial domain is no longer periodic but instead unbounded 6 . This transforms the PDE model (1.1) into The small parameterε allows us to use Geometric Singular Perturbation Theory (GSPT) [11, 20, 22, 23, e.g.] to construct solutions that, to leading order in the small parameter, approximate the localised structures of the main text.
In GSPT, the observation that the dynamics centres around the branches of u-nullcline is taken to the extreme and we construct solutions whose slow dynamics in the singular limit, i.e. in the limit of the small parameterε to zero, is confined to this nullcline, which we will refer to as the slow or critical manifold. In contrast, during the fast jump in u, the slow component will not change in this singular limit. These assumptions simplify the computations and allow us to compute parts of the solution in the singular limit. The main theorems of GSPT [11, 20, 22, 23, e.g.], sometimes called Fenichel Theorem 1-3, allow us to conclude that if the small parameter is small enough 7 , then there indeed is a true solution of the PDE close to the one constructed in the singular limit.
Here, we only show the construction of the standing waves we found in §3.1. That is, we are interested in the fixed points of the PDE dynamics 0 =ε 2 ∂ xx u − a 1 u − a 2 uv + a 3 u 2 a 4 + u 2 + a 5 0 = ∂ xx v + ε(−c 1 v + c 2 u). We refer to this set of equations as the slow system. This system should be understood in the following sense: on a large spatial scale, the dynamics of (v, q) is approximated by lines 3 and 4 of the ODE above, and this approximation is valid in the region of the phase plane given by the algebraic equations in lines 1 and 2. We refer to the solution of these algebraic equations as the slow or critical manifold. When we try to explicitly compute the critical manifold as a function u(v), we encounter a third-order polynomial, which can be solved exactly. However, this is not practical as the graph u(v) cannot be represented by a single function, but it has three branches, the upper, middle and lower branch. We will denote the upper branch with u + (v) and the lower branch with u − (v). Hence, system (B.4) now becomes v =q q =ε(c 1 v − c 2 u ± (v)). (B.5) We refer to this equation as the reduced slow system. In both figures, the red dots are a solution of (B.2), found by using Matlab's bvp4c solver, with an initial condition coming from a PDE simulation and the small parameterε 2 was set to 10 −4 . Note that this value corresponds to Du = 0.01, which is a factor ten smaller than in the main text. In Figure ( When we are interested in the dynamics of u instead of v, we must zoom in to a smaller length scale. Therefore, we defineεξ = x and use˙to denote the derivative with respect to ξ. System (B.3) now becomesu =ṗ p =a 1 u + a 2 uv − a 3 u 2 a 4 + u 2 − a 5 v =εq q =εε(c 1 v − c 2 u).

(B.6)
This system is called the fast system and is still equivalent to (B.3), but when we setε = 0, it is no longer equivalent and reduces tou =ṗ p =a 1 u + a 2 uv − a 3 u 2 a 4 + u 2 − a 5 v =0 q =0.

(B.7)
This shows that in the fast limit, the value of v is constant. When we denote the unknown value bȳ v, the system reduces tou =ṗ p =a 1 u + a 2 uv − a 3 u 2 a 4 + u 2 − a 5 .

(B.8)
This system is known as the reduced fast system.
How can we use both reduced systems to understand the dynamics of Figure B.1b? We observe slow dynamics on the upper and lower branch of the critical manifold and a fast jump in between. The reduced fast system describes the fast jump between the upper and lower branch of the critical manifold. Therefore, a standing wave exists when this system has a heteroclinic orbit between the upper and lower branch. The reduced fast system (B.8) is a Hamiltonian system with Hamiltonian H(u, p) = 1 2 p 2 − 1 2 (a 1 + a 2v ) u 2 + a 3 (u − √ a 4 arctan(u/ √ a 4 )) + a 5 u. We cannot solve this algebraic equation exactly, but it is a straightforward numerical problem. Note thatv only depends on the parameters a 1 , ..., a 5 and not on the parameters of the equation for v.
For this value ofv, the Hamiltonian overlaps with the fast dynamics, as is shown in Figure B.1a. Furthermore, from Figure B.1b, it is clear that the value forv is a good approximation for the location of the jump forε = 10 −2 . Now we have all the ingredients to construct the standing wave. We start at x = −∞ in the background state of the reduced slow system on the lower branch. We follow the dynamics of the reduced slow system (B.5) until we reach the valuev where we jump to the upper branch following the reduced fast system (B.8). We will follow the slow (v, q)-dynamics on the upper branch until we return to the valuev, but with the opposite sign for the derivative, i.e. we trace a curve from (v, q(v)) to (v, −q(v)) in the reduced slow system. Then, we jump down again to the lower branch, which we now trace back to the background state. We implicitly assume here that the maximum value of v remains below the fold of the critical manifold (which is not the case for travelling waves, see Figure 3.9).
This example shows how GSPT can be used to construct localised solutions of (1.1) and also how to understand these solutions. Questions about the existence of localised solutions of (1.1) and bifurcations can thus be reduced to understanding relatively simple ODEs and the connections between them.