Stochastic Turing Pattern Formation in a Model with Active and Passive Transport

Abstract

We investigate Turing pattern formation in a stochastic and spatially discretized version of a reaction–diffusion–advection (RDA) equation, which was previously introduced to model synaptogenesis in C. elegans. The model describes the interactions between a passively diffusing molecular species and an advecting species that switches between anterograde and retrograde motor-driven transport (bidirectional transport). Within the context of synaptogenesis, the diffusing molecules can be identified with the protein kinase CaMKII and the advecting molecules as glutamate receptors. The stochastic dynamics evolves according to an RDA master equation, in which advection and diffusion are both modeled as hopping reactions along a one-dimensional array of chemical compartments. Carrying out a linear noise approximation of the RDA master equation leads to an effective Langevin equation, whose power spectrum provides a means of extending the definition of a Turing instability to stochastic systems, namely in terms of the existence of a peak in the power spectrum at a nonzero spatial frequency. We thus show how noise can significantly extend the range over which spontaneous patterns occur, which is consistent with previous studies of RD systems.

Exact Forms of the Macroscopic and Mesoscopic Models

1.1 Reaction Components of Macroscopic Model

Multiplying ${\mathbf {L}}$ and ${\mathbf {F}}({\mathbf {u}})$ in Sect. 2 yields

$$\begin{aligned} {\mathbf {L}} {\mathbf {F}}({\mathbf {u}}) = \left[ \begin{array}{c} \beta _1 - \mu _1 u_1 + \frac{\rho _1 u_1^2}{u_2+u_3} \\ \frac{\beta _2}{2} - \mu _2 u_2 + \frac{\rho _2 u_1^2}{2} -\alpha u_2 + \alpha u_3\\ \frac{\beta _2}{2} - \mu _2 u_3 + \frac{\rho _2 u_1^2}{2} +\alpha u_2 - \alpha u_3 \end{array}\right] . \end{aligned}$$

(A.1)

Since this multiplication is the local reaction component of the macroscopic model, we have

$$\begin{aligned} {\mathbf {G}}({\mathbf {u}}) = \left[ \begin{array}{c} g_1(u_1,u_2,u_3) \\ g_2(u_1,u_2) -\alpha u_2 + \alpha u_3\\ g_2(u_1,u_3)+\alpha u_2 - \alpha u_3 \end{array}\right] , \end{aligned}$$

(A.2)

where

$$\begin{aligned} g_1(u_1,u_2,u_3) = \beta _1 - \mu _1 u_1 + \frac{\rho _1 u_1^2}{u_2+u_3}, \quad g_2(u_1,u_2) = \frac{\beta _2}{2} - \mu _2 u_2 + \frac{\rho _2 u_1^2}{2}. \end{aligned}$$

Thus, the Jacobian of ${\mathbf {G}}$ at ${\mathbf {u}}^*$ becomes

$$\begin{aligned} \nabla {\mathbf {G}}({\mathbf {u}}^*) = \left[ \begin{array}{c c c} g_{1,1} &{} g_{1,2} &{} g_{1,3} \\ g_{2,1} &{} g_{2,2} -\alpha , &{} g_{2,3} + \alpha \\ g_{2,1} &{} g_{2,3} + \alpha &{} g_{2,2} - \alpha \end{array}\right] , \end{aligned}$$

(A.3)

where the partial derivatives $g_{1,m} = \left. \partial g_1/\partial u_m\right| _*$ are

$$\begin{aligned} g_{1,1} = -\mu _1 + \frac{2\rho _1 u_1^*}{u_2^* + u_3^*}, \quad g_{1,2} = -\rho _1 \left( \frac{u_1^*}{u_2^* + u_3^*} \right) ^2, \quad g_{1,3} = g_{1,2}, \end{aligned}$$

(A.4)

and $g_{2,m} = \left. \partial g_2/\partial u_m\right| _*$

$$\begin{aligned} g_{2,1} = \rho _2 u_1^*, \quad g_{2,2} = -\mu _2, \quad g_{2,3} = 0. \end{aligned}$$

(A.5)

In particular, if $\beta _1 = \beta _2 = 0$, then the partial derivatives reduce to

$$\begin{aligned} g_{1,1} = \mu _1, \quad g_{1,2} = - \frac{\mu _1^2}{\rho _1}, \quad g_{2,1} = \frac{\mu _2}{\mu _1} \rho _1, \quad g_{2,2} = -\mu _2. \end{aligned}$$

(A.6)

1.2 Characteristic Equation for the Linearized Macroscopic Equations

The linear operator ${\mathcal {L}}(k)$ in the characteristic equation (22) can be reduced to the form

$$\begin{aligned} {\mathcal {L}}(k) = \left[ \begin{array}{c c c} 2\kappa _1 q &{} 0 &{} 0 \\ 0 &{} \kappa _2(q - i \sqrt{1-(1+q)^2}) &{} 0 \\ 0 &{} 0 &{} \kappa _2(q + i \sqrt{1-(1+q)^2}) \end{array}\right] + \nabla {\mathbf {G}}({\mathbf {u}}^*) \end{aligned}$$

(A.7)

in accordance with the substitution $q = \cos (k) -1$. The corresponding characteristic equation is

$$\begin{aligned} p(q,\lambda ) = -\lambda ^3 + p_2(q) \lambda ^2 + p_1(q) \lambda + p_0(q). \end{aligned}$$

where

$$\begin{aligned} p_j(q) = \sum _{m=0}^2 p_{j,m}q^m. \end{aligned}$$

The coefficients $p_{j,m}$ are given by

$$\begin{aligned} p_{2,0}&= g_{1,1} + 2g_{2,2} - 2\alpha , \nonumber \\ p_{2,1}&= 2 (\kappa _1 + \kappa _2), \nonumber \\ p_{2,2}&= 0, \end{aligned}$$

(A.8)

$$\begin{aligned} p_{1,0}&= 2\alpha ( g_{1,1} + g_{2,2}) -g_{2,2}^2 -2(g_{1,1}g_{2,2} - g_{1,2} g_{2,1}), \nonumber \\ p_{1,1}&= 2\left[ (\alpha - g_{2,2})(2\kappa _1 + \kappa _2) -g_{1,1}\kappa _2 + \kappa _2^2\right] , \nonumber \\ p_{1,2}&= - 4 \kappa _1 \kappa _2, \end{aligned}$$

(A.9)

and

$$\begin{aligned} p_{0,0}&= (2\alpha - g_{2,2})(2g_{1,2}g_{2,1} - g_{1,1}g_{2,2}), \nonumber \\ p_{0,1}&= -2 \left[ (2\alpha - g_{2,2})g_{2,2} \kappa _1 +(\alpha g_{1,1} - g_{1,1}g_{2,2} + g_{1,2}g_{2,1})\kappa _2 + g_{1,1} \kappa _2^2\right] ,\nonumber \\ p_{0,2}&= -4\kappa _1 \kappa _2(\alpha -g_{2,2}+\kappa _2). \end{aligned}$$

(A.10)

Note that if $\beta _1 = \beta _2 = 0$, then the coefficients $p_j(q)$ are independent of $\rho _1$ and $\rho _2$. Moreover, Eq. (A.6) implies that the partial derivatives are independent of $\rho _2$, and only $g_{1,2}$ and $g_{2,1}$ depend on $\rho _1$. However, the latter coefficients always appear in $p_{j,m}$ as the product $g_{1,2}g_{2,1} = -\mu _1\mu _2$, which is independent of $\rho _1$. This establishes the above claim.

The coefficient $p_j(q)$ also determines the Turing–Hopf bifurcation condition (28)

$$\begin{aligned} p_\text {TH}(q) = p_0(q) + p_1(q)p_2(q) = \sum _{m=0}^3 p_{\text {TH},m} q^m. \end{aligned}$$

In order to investigate its behavior over $q \in [-2,0]$, we first compute

$$\begin{aligned} p_{\text {TH},0}&= p_{0,0} + p_{1,0}p_{2,0} = -2 (g_{1,1} + g_{2,2}) \nonumber \\&\quad \left[ 2\alpha ^2 - \alpha (g_{1,1}+g_{2,2}) + (g_{1,1} g_{2,2} - g_{1,2} g_{2,1}) + g_{2,2}^2 - 2 \alpha g_{2,2}\right] . \end{aligned}$$

(A.11)

Using the fact that $g_{2,2} = - \mu _2$ and the stability condition (19)

$$\begin{aligned} g_{1,1} + g_{2,2} <0, \quad g_{1,1}g_{2,2} -g_{1,2}g_{2,1} >0, \end{aligned}$$

one can show $p_{\text {TH},0}> 0$. Similarly, we find the sign of higher-order coefficients. The leading coefficient is

$$\begin{aligned} p_{\text {TH},3} = -8\kappa _1 \kappa _2 (\kappa _1 + \kappa _2), \end{aligned}$$

(A.12)

which is negative. The next leading coefficient can be written as

$$\begin{aligned} p_{\text {TH},2}&= 4\left[ \kappa _2^3 - (g_{1,1}+g_{2,2})\kappa _2 (2\kappa _1 + \kappa _2) \right. \nonumber \\&\quad + \left. \alpha ( 2\kappa _1^2 + 4\kappa _1\kappa _2 + \kappa _2^2) + 2 \mu _2 \kappa _1 (1 + \kappa _2)\right] , \end{aligned}$$

(A.13)

which is positive. Introducing variables

$$\begin{aligned} \varphi _1 = - g_{1,1} + \mu _2, \quad \varphi _2 = g_{1,1}g_{2,2} - g_{1,2} g_{2,1} \end{aligned}$$

which are all positive, the first-order coefficient takes the form of

$$\begin{aligned} - \frac{p_{\text {TH},1}}{2}&= 4\alpha ^2(2\kappa _1 + \kappa _2) + 2\kappa _2^2 \mu _2 + \alpha [ 2\kappa _2^2 + 6 \kappa _1 ( \mu _2 + \varphi _1) + \kappa _2 ( 3 \mu _2 + 5 \varphi _1)] \nonumber \\&\quad + 2 \kappa _1 (\mu _2^2 + \mu _2 \varphi _1 + \varphi _2 ) + \kappa _2 ( \mu _2^2 + \mu _2 \varphi _1 + \varphi _1^2 + \varphi _2), \end{aligned}$$

(A.14)

which right-hand side becomes positive. Thus, we have $p_{\text {TH},1}<0$. From the above results, we conclude that $p_\text {TH}(q) >0$ for $q \in [-2,0]$.

1.3 Components of Correlation Matrices in the Mesoscopic Model

The nonzero components of ${\mathbf {B}}_{1,n}({\mathbf {u}})$ are

$$\begin{aligned}{}[{\mathbf {B}}_{1,n}]_{1,1}&= \beta _1 + \mu _1 + \rho _1 \frac{u_{1,n}^2}{u_{2,n}+u_{3,n}}, \nonumber \\ [{\mathbf {B}}_{1,n}]_{2,2}&= \alpha (u_{2,n} + u_{3,n}) + \beta _2 + \mu _2 u_{2,n} + \rho _2 u_{1,n}^2, \nonumber \\ [{\mathbf {B}}_{1,n}]_{2,3}&= B_{32}^{1,n} = - \alpha (u_{2,n} + u_{3,n}), \nonumber \\ [{\mathbf {B}}_{1,n}]_{3,3}&= \alpha (u_{2,n} + u_{3,n}) + \beta _2 + \mu _2 u_{3,n} + \rho _2 u_{1,n}^2, \nonumber \\ [{\mathbf {B}}_{1,n}]_{1,2}&= B_{21}^{1,n} = B_{13}^{1,n} = B_{31}^{1,n} =0, \end{aligned}$$

(A.15)

Similarly, the elements of the diagonal matrices ${\mathbf {B}}_{2,n}({\mathbf {u}})$ and ${\mathbf {B}}_{3,n}({\mathbf {u}})$ are

$$\begin{aligned}{}[{\mathbf {B}}_{2,n}]_{1,1}&= \kappa _1 (u_{1,n-1} + 2 u_{1,n} + u_{1,n+1}), \nonumber \\ [{\mathbf {B}}_{2,n}]_{2,2}&= \kappa _2 (u_{2,n-1} + u_{2,n}), \nonumber \\ [{\mathbf {B}}_{2,n}]_{3,3}&= \kappa _2 (u_{3,n} + u_{3,n+1}), \end{aligned}$$

(A.16)

and

$$\begin{aligned}{}[{\mathbf {B}}_{3,n}]_{11}&= -\kappa _1 (u_{1,n-1} + u_{1,n}), \nonumber \\ [{\mathbf {B}}_{3,n}]_{22}&= -\kappa _2 u_{2,n-1}, \nonumber \\ [{\mathbf {B}}_{3,n}]_{33}&= -\kappa _2 u_{3,n}.\end{aligned}$$

(A.17)

1.4 Coefficients in Stochastic Turing Pattern Condition

The first element of the power spectral density is given by Eq. (56), and the coefficients have the following explicit forms:

$$\begin{aligned} \theta _{1,1}&= 2[{\mathbf {B}}_3^*]_{1,1},\nonumber \\ \theta _{1,0}&= [{\mathbf {B}}_1^*]_{1,1} + [{\mathbf {B}}_2^*]_{1,1} + 2[{\mathbf {B}}_3^*]_{1,1}, \end{aligned}$$

(A.18)

$$\begin{aligned} \theta _{2,1}&= - 2\alpha \kappa _2 + 2 g_{2,2} \kappa _2 - 2\kappa _2^2,\nonumber \\ \theta _{2,0}&= g_{2,2}^2 - 2\alpha g_{2,2}, \end{aligned}$$

(A.19)

and

$$\begin{aligned} \theta _{3,2}&= -4\left( \alpha \kappa _1 \kappa _2 - g_{2,2} \kappa _1 \kappa _2 + \kappa _1 \kappa _2^2\right) ,\nonumber \\ \theta _{3,1}&= -2\left( \alpha g_{1,1} \kappa _2 + 2 \alpha g_{2,2} \kappa _1 - g_{1,1} g_{2,2} \kappa _2 +g_{1,1} \kappa _2^2 + g_{1,2} g_{2,1} \kappa _2 - g_{2,2}^2 \kappa _1\right) ,\nonumber \\ \theta _{3,0}&= -2 \alpha g_{1,1} g_{2,2} + 4 \alpha g_{1,2} g_{2,1} + g_{1,1} g_{2,2}^2 - 2 g_{1,2} g_{2,1} g_{2,2}. \end{aligned}$$

(A.20)

Taking the limit $\alpha \rightarrow \infty$ in the stochastic Turing condition (57), we have a quadratic Eq. (60) and its coefficients are the following:

$$\begin{aligned} \zeta _2&= 2 g_{2,2}^2 u_1^*, \nonumber \\ \zeta _1&= \mu _1(1+ u_1^*)g_{2,2}^2 + 2 g_{1,2} g_{2,1} \kappa _2 u_1^* + g_{1,1} g_{2,2}^2 u_1^* - 2 g_{1,2}g_{2,1}g_{2,2}u_1^*, \nonumber \\ \zeta _0&= \mu _1(1+u_1^*)g_{1,2}g_{2,1} \kappa _2. \end{aligned}$$

(A.21)