A simple reaction–diffusion system as a possible model for the origin of chemotaxis

Chemotaxis is a directed cell movement in response to external chemical stimuli. In this paper, we propose a simple model for the origin of chemotaxis – namely how a directed movement in response to an external chemical signal may occur based on purely reaction–diffusion equations reflecting inner working of the cells. The model is inspired by the well-studied role of the rho-GTPase Cdc42 regulator of cell polarity, in particular in yeast cells. We analyse several versions of the model to better understand its analytic properties and prove global regularity in one and two dimensions. Using computer simulations, we demonstrate that in the framework of this model, at least in certain parameter regimes, the speed of the directed movement appears to be proportional to the size of the gradient of signalling chemical. This coincides with the form of the chemical drift in the most studied mean field model of chemotaxis, the Keller–Segel equation.


Introduction
Chemotaxis is a directed cell movement in response to external chemical stimuli.Chemotaxis is ubiquitous in biology; for example, it plays a role in organism morphology [1,2], reproduction processes [4,5,6,7] and workings of immune system [8,9].There are many mathematical models of chemotaxis; the most studied is the Keller-Segel equation and its variants.Virtually all of these models incorporate a transport term driven by the concentration of the external chemical, which may be produced by the cells themselves (see e.g.[10,11], where further references can be found).Yet we are not aware of any mathematical models that would aim to explain the origin of the transport based on reaction-diffusion processes taking place inside cells.The way chemotaxis happens, at least for eukaryotic cells, is that cells translate chemical environmental cues into amplified intracellular signaling that results in elongated cell shape, actin polymerization toward the leading edge, and movement along the gradient.In this paper, instead of presenting chemotaxis as an explicit transport term, we explore model that aims to explain the origin the chemotactic ability of cells.Inspired by [12], we look at sexual yeast reproduction and simplify the polarization process into understanding active rho-GTPase Cdc42 concentration in one yeast under chemical gradient produced by another yeast.This is certainly just an element of a more complex picture involved into producing chemotactic response in cells, but we limit consideration to this one stage.Our first goal is to explore the well-posedness properties of the model and its variants and to understand analytic features involved.Our second goal is to get more information on the nature of transport generated by the model reacting to external chemical stimuli.In particular a natural question is how the speed of transport, which we measure via the coordinate moment of a density, depends on the gradient of the attractive chemical.This question we approach through numerical simulations, and find that for certain reasonable ranges of parameters, this dependence is linear.
The two-species mass-conserved activator-substrate (MCAS) system that is the basis of our model consists of two partial differential equations (PDEs) governing the kinetics of the slowly diffusing activator u (GTP-bound GTPase on the membrane) and the rapidly diffusing substrate v (GDP-bound GTPase in the cytoplasm).In general, this system has the following form in 1D (see [12]): (1.1) Here, k refers to the diffusion of u, k v refers to the diffusion of v.These two diffusion constants usually differ by two orders of magnitude.F (u, v) describes the biochemical interconversions between u and v, given in the form: Goryachev's (see [13]) and (see [14]) In [12], Turing-type instability for these types of reaction have been analyzed; it was also shown that steady states with more than one peak are unstable for many kinds of F (u, v).This analysis is in agreement with experiments as one usually only observe one bud in yeast asexual production [12].
Several complicated computational models have been developed to mimic gradient-induced polarization toward the pheromone source [15,16] and have shown the rate of movement is dependent on pheromone concentration [17].Here, we propose a simpler system that is capable of capturing such gradient tracking ability, specifically, in the context of chemotactic reaction of a single yeast cell to an external pheromone signal.We apply a modification to the Turing-type model described above and add a pheromone density profile αf (x) that depends on location in the form similar to [15] -and we obtain the following system: In (1.3), ã, b, α, k, k v are constants.ã is the reaction activation constant, b is the reaction depletion constant, α is the overall pheromone strength, k is the diffusion coefficient for u, and k v is the diffusion coefficient for v. f (x) is a bounded smooth function that describes the pheromone level at different locations.
We assume that rho-GDPase diffuses infinitely fast, i.e, k v approaches ∞.Since the total mass of rho-GTPase and rho-GDPase is conserved, M = (u(•, t) + v(•, t))dx is a constant.Then we can obtain the following equation (1.4) that describes the activatorsubstrate reaction between these two substances.The setting we have is x ∈ T d when dimension d = 1, 2, with periodic boundary condition: In (1.4), |T d | is the measure of the domain, M is the total mass.We are interested in the non-negative solution u with udx ≤ M for all time.By rescaling space and time, we can simplify the equation (1.4) as follows: (1.5) Here depletion rate and diffusion coefficient are normalized to 1, and  Our main results are as follows.On the rigorous level, we are able to establish global regularity results for equation (1.5) in one and two dimensions for all non negative initial data.To better understand the structure of the equation, we consider (1.5) in the absence of regularizing diffusion, and prove that for all finite 0 ≤ t ≤ T ≤ ∞, the profile for active rho-GTPase u stays smooth even when there is no diffusion term.When we reintroduce diffusion in one dimension, uniform in time global bound on derivatives of u is shown.With diffusion in two dimension, global regularity with possible growth is proved.In our numerical experiments, we observe the profile of active rho-GTPase u move towards higher concentration of pheromone, and stops moving once it reaches the location with maximum pheromone concentration.In addition, we explore the speed of such movement through tracking the center of mass of rho-GTPase profile.If pheromone concentration is linear, the center of mass moves with a constant speed towards the pheromone peak.More importantly, we verify that the movement speed depends linearly on the pheromone gradient in a natural parameter range similar to that used in [12].Note that such linear dependence of chemotactic drift on the gradient of the attractive chemical density f (x) is a common feature of chemotaxis models, including the most studied Keller-Segel equation which in its simplest form reads (see e.g.[11]) The emergence of transport mean field equations such as (1.6) from kinetic equations has been extensively studied (see e.g.[18,19,20,21]).However the existence of chemotactic transport is already built in on the kinetic level.As far as we know, the equation (1.5) is the first simple reaction-diffusion model that aims to rigorously analyze the emergence of chemotaxis from the inner cell workings, even if it is focused on just one stage of the process that can be quite complex.
The paper is organized as follows: in the next section, we introduce the general set up and key parameters of the model in more detail and present our numerical scheme.We then proceed to describe results of the numerical experiments.After this we state the rigorous results that we are able to prove, and proceed with the proofs.

General Set Up and Numerical Scheme
We want to explore the origin of chemotactic ability of cells with simulation in 1D using the following equation (dropping the scripts in (1.4) and denote the total mass of u(x, t) as The parameters we use are the same as in [12] and are shown in Table 1 with some basic conversions.We used the method of lines to turn spatially discretized PDE into a system of ODEs, then we use a robust ODE solver ODE15s in Matlab to solve.Note that since we assume rho-GDPase v is rapidly diffusing, we use Simpson's method to numerically integrate rho-GTPase to obtain u(x, t)dx, and calculate rho-GDPase as follows: For the computational part, we restrict ourselves to one dimensional surface and assume the pheromone profile is generated by another yeast cell.If this cell is some distance away, one reasonable model of the two-dimensional pheromone distribution is a solution to the heat equation ∂ t ω − ∆ω = 0 with a δ function initial data, that is, as a fundamental solution of 2D heat equation.Then we can derive the pheromone profile on the cell boundary as shown in Figure 3, and in general, it has the form: where φ = x−x peak r and φ ∈ [−π, π).We use γ to denote the diffusion coefficient for the source, β as the response strength to the source, and without loss of generality, we assume t = 1.With γ = 10, t = 1, L = 10, β = 1500, we can plot f h (x) in Figure .4. We also plot a linear pheromone profile f (x) that is similar to f h (x) with a peak at x peak defined below in Figure 4 as well.The reason for this choice will be explained later. (2.4) Parameters from [12] To obtain the initial profile of u, we start with a uniformly distributed v and a bump function u, and we run the simulation without pheromone until it stabilizes according to (2.5) Switching on the pheromone, we track the movement using center of mass since the profile of u is relatively stable over time.The center of mass as a function of time is defined using: We also measure the movement of the profile of u with the time derivative of the center of mass: dCMu(t) dt .

Numerical Results: Pheromone Induced Movement
While there is no explicit transport term in (2.1), we observe movement of the rho-GTPase u profile over time in response to "reallocation of resources" to more favorable reaction regions induced by pheromone αf h (x) and αf (x) as shown in Figure 5.As our numerical simulations show, at least in certain parameter ranges, this transport appears quite similar to the Keller-Segel-type transport with speed proportional to the concentration gradient.
As we can see in Figure 5, one does not expect that the solution will be an exact traveling pulse since the background level of the pheromone affects the shape of the bump of u density.With time t and recorded CM u (t), we can compute the movement speed of the center of mass, dCMu(t) dt and the corresponding pheromone profile f h and f at the center of mass.From Figures 6 and 7, we propose the hypothesis that the movement speed dCMu(t) dt depends on the derivative of the pheromone.
We have tried some other pheromone profiles with similar results.While the graphs on 7 appear close to lines, they are not quite lines -but perhaps because the density bump of u has spatial scale, and so is exposed to a range of concentration slopes (we take the slope at the center of mass as the basis of the functional relationship pictured in this figure).As one of our goals is to study the dependence of the movement speed of the center of mass CM u (t) and on the gradient of the pheromone concentration αf (x), we will also consider the piecewise linear pheromone profile f (x) : it has extended regions with the constant slope (A) Pheromone profile given by f h (x).
that makes it easier to capture the effect more precisely.We illustrate this transport picture in Figure 8 by calculating the center of mass using (2.6) of the initial profile of u and the profile at t = 10000s pheromone profile f (x) with pheromone strength α = 2.
As one can expect, the profile of u slows down once its center of mass starts to approach x = x peak = 5µm.We can plot the center of mass a function as a time of time for different pheromone strength α.As presented in Figure 9, the center of mass of u stays at x = x peak = 5µm after t = 7000s for the pheromone strength α = 3.In fact, if we run the simulation long enough, the center of mass of u appears to get arbitrarily close to x = x peak = 5µm for all α > 0.  We then continue to explore the movement speed of the center of mass as a function of pheromone strength α, which controls the slope of the pheromone concentration.From Figure 9, we can see a constant movement speed of the center of mass when the profile of u is far away from x peak .Moreover, if we plot the movement speed (before the profile of u is too close to x peak ) as a function of α as shown in Figure 10, we observe linear dependence of transport speed on pheromone strength.Such linear dependence corresponds to the mean field chemotaxis models in [10,11].

Our goal in this section is to initiate rigorous mathematical analysis of the equation (1.5).
There is much literature on regularity of solutions for reaction-diffusion type equations, but we could not locate references dealing directly with (1.5) due to nonlocality produced by our assumption of infinite diffusion for v.In [22], the author presented global existence result to a similar equation for sufficiently small and non-negative initial data.In this section, we will present regularity result for (1.5) in T d , with d = 1, 2 for any non-negative smooth initial data with diffusion and without diffusion.When there is no diffusion with d = 1, 2, we can still show that u is smooth for all times 0 ≤ t < ∞.With diffusion in 1D, uniform in time global bound is proven while with diffusion in 2D, global regularity with possible growth is shown.At this time, we are unable to prove rigorous results that would provide  more qualitative features of the evolution, and in particular establish the conjecture on linear dependence of the transport speed on pheromone gradient at least in some regimes.Part of the difficulty is that, as we mentioned before, the solution cannot be expected to take form of a traveling pulse of the fixed shape on an unbounded domain, and lack of such framework complicates analysis.This paper can be viewed as creating a foundation for such further investigation that remains an intriguing open problem.4.1.Without diffusion.We start with the case without diffusion.Consider the following equation: Theorem 1. Suppose u(t, x) is a non-negative solution to (4.1) with dimension d = 1, 2 and periodic boundary condition; a, α, M are constant parameters, and f (x) is a smooth nonnegative function.If u 0 (x) is a smooth initial profile, then u(t, x) stays smooth for all times 0 ≤ t < ∞.
Proof.We first show some a-priori bounds on u, namely that all Sobolev norms of u are controlled by L ∞ norm.It is clear that on a bounded domain, L ∞ norm controls all other L p norms.In the estimates that follow, D stands for any partial derivative (just ∂ x in one dimension), and W k,p is the usual Sobolev space.Multiplying (4.1) by (−∆) s u and integrating, we get The last term can by controlled by αM ||f || W 2s,1 ||u|| L ∞ (moving all derivatives to f ).Term I can be represented by a sum of integrals of the type D l uD s−l uD s udx, where l = 0, . . ., s.
Then with Hölder's inequality and Gagliardo-Nirenberg interpolation inequality in dimension d = 1, 2, we can bound them by: p(−1+2s) and β = 2(−1+(s−l))q q(−1+2s) , and α + β = 2( 1 p + 1 q ) = 1.In 2D, α = −2+lp p(−1+s) and β = −2+(s−l)q q(−1+s) , and Substituting back into (4.2) gives: [3], we see that to show global regularity, it suffices to prove that T 0 ||u(•, t)|| ∞ dt remains bounded.We will show a stronger constraint that u(•, t) ∞ remains finite for all times.Via contradiction, denote T the first time of blow up of ||u|| ∞ .Consider ρ = e −t 1+u , then ρ satisfies the following equation: Then at time T , ρ will reach 0 at some point.For simplicity, let us focus on d = 1; the argument for d = 2 is similar.From (4.3), we see that ρ is decreasing for all times, therefore ρ is bounded from above by ||ρ 0 || ∞ .Now we take a derivative with respect to x to obtain: Then it is clear that for all x, By Gronwall inequality, we can see that |∂ x ρ| is bounded for any finite time t ≤ T .We can take another derivative in x and obtain: Boundedness of ρ along with control of |∂ x ρ| give: xx ρ| is bounded for any finite time t ≤ T .One can effectively continue this calculation and get that all derivatives in space are bounded for t ≤ T .Since blow up happens for the first time at time T , then ρ(x B , T ) = 0 at some point x B .There can be many such points, but let us focus on one of them.Due to control of ∂ x |ρ| and ∂ 2 xx |ρ|, we have ρ(x, T ) ≤ C(T )|x − x B | 2 by Taylor expansion.Observe that ρ(x, t) → ρ(x, T ) monotonically for every x.Therefore, as u = 1 e t ρ − 1, we have u(x, t) → u(x, T ) (including when u(x, t) = ∞).Then we have u(x, T ) = 1 e T ρ − 1 ≥ Then by Fatou's lemma, we have (4.4) which is a contradiction.Therefore, we cannot have such finite time blow up.Note that the argument above works both in 1D and 2D, only the computation yielding control of the derivatives of ρ needs a minor adjustment.Note that the size of the Sobolev norms of the solution may depend on time.

4.2.
With diffusion in 1D.Now we turn out attention to the system with diffusion in 1D: udx.
We will prove global regularity for (4.5) as well.
Theorem 2. Suppose u is a non-negative solution to (4.5) in dimension d = 1 and periodic boundary condition.Let a, α, M be constant parameters, and f (x) a smooth function.If u 0 (x) is a smooth initial profile, then u(x, t) stays smooth for all time; in particular, all Sobolev norms ||u|| H s with s > 0 are bounded uniformly for all time.
Proof.First we show that ||u|| 2 is bounded: multiplying both sides by u and integrating, we have Using Gagliardo-Nirenberg-Sobolev inequality (see eg. [23]) gives: Substituting (4.7) into (4.6)yields (note that the constant C changes from line to line and may depend on M ): Note that in the second inequality above, we used Young's inequality ab ≤ a p p + b q q with p = 3, q = 3 2 .The calculation (4.8) implies that ||u|| L 2 is globally bounded by Gronwall inequality.Moreover, from (4.8), we can see that ||u|| 2 is in fact uniformly bounded by M 2 ||f || ∞ + C(M ) since if ||u|| 2 ever crosses this value for the first time at t 0 , ∂ t ||u|| 2 becomes negative, which implies that before t 0 , ||u|| 2 lies above M 2 ||f || ∞ + CM 5 already.We arrive at a contradiction.Therefore, ||u|| 2 is uniformly bounded for all time.Next, we estimate the higher order Sobolev norms.In fact, in 1D, this can be done using just control of the L 1 norm, but we are going to use L 2 norm for convenience as we have shown it remains bounded.Multiplying (1.5) by (−∆) s u and integrating, we get The last term can by controlled by αM ||f || H 2s ||u|| L 2 (moving all derivatives to f ).Term I can be represented by a sum of integrals of the type D l uD s−l uD s udx, where l = 0, . . ., s.In 1D, D = ∂ x .Then with Hölder's inequality and Gagliardo-Nirenberg interpolation inequality with dimension d = 1, we can bound them by: Here we deploy the Gagliardo Nirenberg inequalities . Substituting gives (note that the constant C changes from line to line): where we use Young's inequality ab ≤ a p p + b q q with p = 4s+4 4s+1 , and q = 4s+4 3 .Given that we proved u 2 is bounded uniformly in time, (4.11) implies that ||u|| H s is also bounded uniformly for all time.

With Diffusion in 2D.
The equation that we are interested in is given by: (4.12) Theorem 3. Suppose u is a non-negative solution to (4.12) with dimension d = 2 and periodic boundary condition.Let a, α, M be constant parameters, and f (x) a smooth function.If u 0 (x) is a smooth initial profile, then u(x, t) stays smooth for any finite time, that is, Sobolev norms ||u|| H s with s > 0 are bounded for any time 0 ≤ t < ∞.The bound on the Sobolev norms may now depend on time.
Proof.First we derive an a-priori estimate.In the 2D case, the analog of the estimate (4.8) is not available, as the exponents do not allow to control u L 2 uniformly in time in this way.Therefore, we need a more nuanced argument.Note that The integral I is a sum of terms of the form: D l uD s−l uD s udx, and l = 0, . . ., s.We estimate it as follows: , α + β = 1.Substituting these estimates back into (4.15)gives:

Figure 1 .
Figure 1.Local activation via positive feedback and depletion of the substrate in the cytosol generates an activator-enriched domain on the cortex [12].

Figure 2 .
Figure 2. The interconversions of Rho-GTPases between active and inactive forms can be modeled as an reaction-diffusion equation governing the dynamics of the slowly-diffusing activator u and the infinitely-diffusing substrate v.

Figure 3 .
Figure 3. Derivation of a pheromone profile generated from a heat equation.

Figure 4 .
Figure 4. Pheromone profile generated from a heat equation f h (x) and a similar linear pheromone profile f (x).

Figure 9 .
Figure 9. Center of mass position CM u (t) as a function of time with pheromone profile f (x).In regions away from x peak , CM u (t) changes linearly with time.

Figure 10 .
Figure 10.Center of mass movement speed dCMu(t) dt as a function of pheromone strength α with pheromone profile f (x).In regions away from x peak , dCMu(t) dt