Terminaison mechanisms of Turing patterns in growing systems

TL;DR

This paper investigates how Turing patterns in growing systems can naturally terminate through parameter changes, combining analytical and numerical methods to propose a controllable mechanism for pattern cessation.

Contribution

It introduces a novel mechanism for Turing pattern termination based on local concentration increases, supported by analytical stability analysis and numerical validation.

Findings

A local increase in reservoir concentration can stop Turing patterns.

Analytical conditions for pattern termination are derived.

Numerical simulations confirm the proposed termination mechanism.

Abstract

The question of the termination of a periodic spatial structure of Turing type in a growing system is addressed in a chemical engineering perspective and a biomimetic approach. The effects of the dynamical parameters on the stability and the wavelength of the structure are analytically studied and used to propose experimental conditions for which a Turing pattern stops by itself with a decreasing wavelength. The proposed mechanism is successfully checked by the numerical integration of the equations governing the dynamics of the activator and the inhibitor. We conclude that a local increase of the concentration of the reservoir which controls the injection rate of the inhibitor into the system can be used to achieve the appropriate termination of a Turing pattern.

Equations

$$\mathrm{A}$$

$$2A+\mathrm{B}$$

$$\mathrm{B}$$

$$\frac{\partial A}{\partial t}$$

$$\frac{\partial B}{\partial t}$$

$$\Delta =(k_{-3}{R}_2{)}^2-4k_1^2{k}_3/k_2\ge 0$$

$$A_T$$

$$B_T$$

$$M=\left(\begin{array}{cc}{k}_1-D_A{q}^2& k_2{A}_T^2\\ -2k_1& -\frac{k_2{k}_{-3}{R}_2}{k_1}{A}_T-D_B{q}^2\end{array}\right)$$

$$\mu$$

$$q_{\text{max}}=\sqrt{\frac{A_T(D_A+D_B)\sqrt{2k_1{k}_2{D}_A/D_B}-k_1-k_2{k}_{-3}{R}_2{A}_T/k_1}{D_B-D_A}}$$

$$\lambda =\frac{2\pi }{q_{\text{max}}}$$

$${\mu }_{\text{max}}$$

$$A(i,s=0)=A_T,B(i,s=0)=B_T,\mathrm{for}1\le i\le 10$$

$$A(i,s=0)=A_0,B(s,s=0)=B_0,\mathrm{for}11\le i\le n_0$$

$$A(1,s+1)$$

$$B(1,s+1)$$

$$R_2$$

Full text

Corresponding author: Annie Lemarchand, E-mail: anle@lptmc.jussieu.fr

1 Introduction

During embryonic development, segmented structures of the body such as the spine and the digits are formed by the production of repeated elements. Since the seminal work of Turing [1] accounting for the formation of biological pattern in the framework of reaction-diffusion models, experimental evidences of Turing structures have been given in chemistry [2, 3, 4] and biology [5, 6]. Recent years have witnessed a growing number of papers where reaction-diffusion principles are proposed to model the formation of biological periodic spatial structures [7, 8, 9, 10, 11, 12, 13]. Following Turing, a two-component chemical system composed of an autocatalytically-produced activator by consumption of an inhibitor that diffuses faster may produce periodic spatial structures such as stripes in one-dimensional (1D) systems and hexagons in 2D. In other words, a Turing pattern emerges by local self-activation and lateral inhibition [14]. The concepts developed to model living systems provide a source of inspiration in chemical engineering [15, 16, 17, 18, 19, 21, 22]. However, standard models of Turing patterns generate structures in infinite systems and the question of the termination of a striped structure in a finite system arises in a perspective of biomimetism in material science. Specifically, it is desirable to find experimentally achievable conditions creating a finite-size structure, whose growth stops by itself with decreasing oscillation amplitude and respects the decrease of the wavelength during the termination process. To this aim, we use an as simple as possible reaction-diffusion model [23] admitting a Turing structure developing behind a propagating wave front and examine the effect of all parameters on both the stability and the wavelength of the structure [5, 22]. We already used a stochastic approach to a Turing pattern [23] and showed that, contrary to intuition, the internal fluctuations may have a constructive effect and stabilize the structure outside the domain of stability predicted by a deterministic description. Here, we adopt a macroscopic approach. Our goal is to select suitable conditions from this systematic approach and to propose termination mechanisms compatible with processing in chemical engineering.

The paper is organized as follows. In section 2, a reaction-diffusion model involving local activation and long-range inhibition is presented. An analytical stability condition and the wavelength expression of the Turing pattern are given. The influence of the parameters of the model on the stability and the wavelength of the pattern are studied in section 3. The analysis of the results leads to the selection of a user-friendly termination mechanism in the framework of chemical engineering. The analytical predictions regarding stability and wavelength are compared to numerical results for the chosen mechanism. Section 4 contains conclusions. The possibility that the different mechanisms exhibited could be found as termination scenarios in morphogenesis is raised.

2 Model

We consider the following reaction mechanism inspired from the Schnakenberg model [24] and the Gray-Scott model [25]

[TABLE]

where R1 and R2 are reservoirs. The concentrations $R_{1}$ and $R_{2}$ of the reservoirs are assumed to be constant in time. The reaction given in Eq. (8) autocatalytically produces species A and consumes species B. Due to the ability of accelerating its own production, species A is called an activator whereas species B, which is consumed by the same process, is called an inhibitor. The macroscopic dynamics of the system is governed by two partial differential equations [9, 23]

[TABLE]

for the concentrations $A$ and $B$ of the activator and the inhibitor supposed to have different diffusion coefficients $D_{A}$ and $D_{B}$. For appropriate rate constant values, such that

[TABLE]

the system admits two steady states $(A_{0}=0,B_{0}=k_{-3}R_{2}/k_{3})$ and

[TABLE]

that are stable with respect to homogeneous perturbations. A linear stability analysis of Eqs. (13,14) reveals that the steady state $(A_{T},B_{T})$ can be destabilized by inhomogeneous perturbations [5, 3, 9, 23].

The Fourier transforms $A_{q}(t)=\int_{-\infty}^{\infty}A(x,t)e^{-iqx}dx$ and $B_{q}(t)=\int_{-\infty}^{\infty}B(x,t)e^{-iqx}dx$, where $q$ is the Fourier mode, are introduced. In the Fourier space, the linear stability operator $M$ is given by:

[TABLE]

The eigenvalues of the matrix $M$ obey the characteristic equation $\mu^{2}+\alpha\mu+\beta=0$, with $\alpha=k_{1}-\frac{k_{2}k_{-3}R_{2}}{k_{1}}A_{T}-(D_{A}+D_{B})q^{2}$ and $\beta=2k_{1}^{2}A_{T}/B_{T}-(k_{1}-D_{A}q^{2})(k_{-3}R_{2}/B_{T}+D_{B}q^{2})$. The Turing structure develops if the largest eigenvalue

[TABLE]

is real and positive [5, 3]. Indeed, a system of differential equations, linearized around a homogeneous steady state, is easily solved by diagonalizing the linear operator. Then, the solution is a linear combination of eigenmodes which exponentially depend on time according to the corresponding eigenvalues. A term associated with a real, positive eigenvalue grows in time, leading to trajectories in the concentration space that move away from the fixed point [5]. In the studied system, the destabilization of the steady state occurs in favor of a Turing pattern. Equation (19) imposes conditions on the parameter values. In particular, the diffusion coefficient $D_{B}$ of the inhibitor B must be sufficiently larger than the diffusion coefficient $D_{A}$ of the activator A: The destabilization of the homogeneous steady state $(A_{T},B_{T})$ requires local self-activation, ensured by the autocatalytic production of the activator through the reaction given in Eq. (8), as well as long-range inhibition, due to the larger diffusion coefficient of the inhibitor. The mode $q_{\textrm{max}}$, which maximizes the eigenvalue $\mu$, is the most unstable Fourier mode:

[TABLE]

In order to account for the termination of the Turing pattern in a growing system, including the fact that the structure ends with a gradually shorter spatial oscillation, we need to find conditions for which the structure tends to lose its stability while its wavelength decreases. The wavelength of the periodic structure is given by:

[TABLE]

Turing structure becomes unstable as the value of largest eigenvalue vanishes for the mode $q_{\textrm{max}}$ associated with the maximum of $\mu$:

[TABLE]

with $q_{\textrm{max}}$ given in Eq. (20). Figure 1 illustrates the behavior of $\mu_{\textrm{max}}$ for parameter values associated with a stable Turing pattern with $\mu_{\textrm{max}}>0$. It is also shown that it is sufficient to increase the value of $k_{-3}R_{2}$ to shift the curve $\mu(q^{2})$ down and lose the stability of the Turing pattern.

In the next section we investigate the behavior of $\lambda$ and $\mu_{\textrm{max}}$ as each parameter controlling dynamics varies. Specifically, we aim at identifying diffusion coefficients or rate constants whose variation leads both to a decrease of the wavelength and a destabilization of the Turing structure, i.e. negative values for the maximum of the eigenvalue.

3 Results

The concentration $R_{2}$ of the inhibitor reservoir is first assumed to be homogeneous. Figures 2 and 3 show the variations of the wavelength $\lambda$ and the maximum value $\mu_{\textrm{max}}$ of the eigenvalue with one of the diffusion coefficients, the other parameters being constant. The results are deduced from Eqs. (20) and (21) for $\lambda$ and Eq. (22) for $\mu_{\textrm{max}}$, the expressions of the steady state $(A_{T},B_{T})$ being given in Eqs. (16) and (17). To facilitate the comparison with the numerical integration of Eqs. (13) and (14) that will be performed in the following, the wavelength is given in number of spatial cells of length $\Delta x=0.1$. As shown in Fig. 2, the decrease in the maximum value $\mu_{\textrm{max}}$ of the eigenvalue as the diffusion coefficient $D_{A}$ of species A increases is accompanied by an increase of the wavelength $\lambda$: The loss of stability of the Turing structure occurs with an increase of the spatial period. We conclude that a variation of the diffusion coefficient $D_{A}$ cannot be argued as a justification of the termination process. The behavior with respect to the diffusion coefficient $D_{B}$ of species B is different. The simultaneous loss of stability of the structure and the decrease of the wavelength are observed in Fig. 3 as $D_{B}$ decreases: The diffusion coefficient $D_{B}$ of species B can be considered as a suitable parameter in the search for a termination model.

Figures 4-7 show the variations of the wavelength $\lambda$ and the maximum value $\mu_{\textrm{max}}$ of the eigenvalue with rate constants. The variations of $\lambda$ are given in a bounded interval of rate constant values, in which the Turing pattern is stable. At one of the endpoints of the interval, the eigenvalue $\mu_{\textrm{max}}$ vanishes and at the other endpoint, the condition of existence of the steady state $(A_{T},B_{T})$ given in Eq. (15) is no longer satisfied. The two desired behaviors, i.e. the decrease of both $\lambda$ and $\mu_{\textrm{max}}$, are observed as $k_{1}$ decreases, $k_{2}$ increases, $k_{3}$ decreases, and $k_{-3}R_{2}$ increases. For an assumed homogeneous concentration $R_{2}$ of the reservoir, the variations of $\lambda$ and $\mu_{\textrm{max}}$ with $R_{2}$ are analogous to the variations with $k_{-3}R_{2}$. According to the chemical reaction given in Eq. (4), decreasing the rate constant $k_{1}$ tends to increase the concentration of species A. Following Eq. (8), increasing the rate constant $k_{2}$ of the autocatalytic step tends to increase the concentration of species A and decrease the concentration of species B. This last result seems to be inconsistent with the consequences drawn from the decrease in $k_{3}$ or the increase in $k_{-3}R_{2}$, which result in increasing the concentration of species B according to Eq. (12). However, we already stated that increasing $B$ through soliciting the reservoir $R_{2}$ results in consuming species B faster by the autocatalytic step given in Eq. (8) [9, 26]. In particular, we observed that introducing a local source of species B leads to the nonintuitive local decrease of B concentration. Hence, all the variations of the rate constants that lead to a loss of stability of the Turing pattern are eventually associated with an increase of A concentration and a decrease of B concentration.

The diffusion coefficients and the rate constants characterize dynamics and are intrinsic to the reaction-diffusion system. Nevertheless, it is always possible to imagine spatial variations of the dynamical parameters. Well-chosen variations of the diffusion coefficient $D_{B}$ of the inhibitor and each of the four rate constants of the chemical mechanism could be a priori used to build a termination model. In the framework of the application to developmental biology, steric hindrance and molecular crowding in the tail of an embryo may be invoked to justify the decrease of the diffusion coefficients. In chemical engineering, a local increase of temperature could be used to induce a local increase of the rate constants. However, local increase of confinement or temperature is susceptible to simultaneously affect several dynamical parameters [13, 22, 27, 28, 29, 30, 31, 32]. Whereas a decrease of $D_{B}$ is desired to destabilize the Turing pattern while decreasing its wavelength, a simultaneous decrease of $D_{A}$ would be detrimental. Similarly, an increase of $k_{2}$ and $k_{-3}$ due to temperature increase could be satisfying but the joint decrease of $k_{1}$ and $k_{3}$ could blur the effect on the Turing structure. The simplest way to imagine the control of a targeted parameter leading to the desired behavior is to impose well-chosen spatial variations of the reservoir concentration $R_{2}$. Indeed, the product $k_{-3}R_{2}$ plays the role of an apparent rate constant for the backward reaction given in Eq. (12) that can be fixed at will in chemical engineering in the case of a single dynamical system with uniquely defined intrinsic parameters.

According to Fig. 7, increasing $R_{2}$ tends to destabilize the Turing pattern and decrease its wavelength. We examine if the results deduced from a stability analysis can be used in a dynamical approach. The results of the numerical integration of Eqs. (13) and (14) for a homogeneous concentration $R_{2}$ and a piecewise linear profile are given in Fig. 8. The initial condition is a step function between the steady state $(A_{T},B_{T})$ in the first $10$ cells on the left and the steady state $(A_{0},B_{0})$ in the remaining cells. The initial total number of cells is set at $n_{0}=610$. Introducing the cell label $i=x/\Delta x$, where $\Delta x$ is the cell length, and the discrete time $s=t/\Delta t$, where $\Delta t$ is the time step, we choose:

[TABLE]

To account for the growth of the system and simultaneously avoid boundary effects that may alter the wavelength of the structure [16], we impose a fixed boundary on the left and a free growing end on the right [9, 23, 26]. For parameter values for which the steady state $(A_{T},B_{T})$ is unstable with respect to inhomogeneous perturbations, a Turing pattern develops after the passage of a chemical wave front. More precisely, according to Eqs. (13) and (14) and due to the no-flux boundary conditions applied on the left boundary, the concentrations in the first cell obey

[TABLE]

so that both $A(1,s)$ and $B(1,s)$ are extremum of the Turing pattern in the first spatial cell $i=1$.

Spatial cells are added to the right end of the system at the front speed to counterbalance the progression of the wave front and mimic system growth: At all the discrete times $s$ for which the concentration $B(n-600,s)$ of species B in the $n-600$ cell becomes smaller than $0.99B_{0}$, the total number $n$ of cells is increased by 1. Provided that the front propagates at a speed smaller than $\Delta x/\Delta t$, this protocol ensures that a layer of about $600$ cells remains in the stationary state $(A_{0},B_{0})$ on the right of the system, so that the propagation of the front is not significantly affected by the finite size of the system. To draw Fig. 8b, we have chosen the parameter values given in the caption of Figs. 1 and imposed $k_{-3}=8.76$ for the following spatial profile for the concentration $R_{2}$ of the inhibitor reservoir:

[TABLE]

The simulation is stopped at time $t_{end}$ for which the wave front has passed cell $i=1000$. It is worth noting that the Turing pattern is unchanged for larger values of the final integration time. Then, only the position of the concentration gradients associated with the traveling wave evolve in time but the Turing structure has stopped growing and remains in a steady state with a fixed number of wavelengths. As desired, the increase of the concentration $R_{2}$ leads to the termination of the Turing structure.

As illustrated in Fig. 7, the Turing structure is expected to be stable in the range $1\leq i<500$ for which $k_{-3}R_{2}=8.76$ and unstable in the range $i\geq 1000$ for which $k_{-3}R_{2}=10$. More precisely, according to Eq. (22), the maximum of the eigenvalue $\mu_{\textrm{max}}$ vanishes for $k_{-3}R_{2}=9.75$, i.e. $R_{2}=1.11$ for $k_{-3}=8.76$, which occurs in spatial cell $i=900$. Hence, the Turing pattern is predicted to be stable in the range $0\leq i<900$ and unstable beyond this domain. The results shown in Fig. 8b confirm the analytical predictions. The amplitude of the spatial oscillations decreases between $i\simeq 500$ and $i\simeq 1000$. The system is in a steady state in the range $1000\leq i<1500$.

The increase of $R_{2}$ not only destabilizes the Turing structure but also modifies the steady state values and the propagation speed of the wave front. The comparison between Figs. 8a and 8b shows that, as $R_{2}$ increases, the wave front propagates faster, $A_{T}$ increases, $B_{T}$ decreases and $B_{0}$ increases. As a consequence of the variation of $A_{T}$ and $B_{T}$, the oscillations of A and B concentrations are not symmetrical in the range $500\leq i<900$. The decrease of the wavelength predicted in Fig. 7 is more difficult to check by a qualitative analysis. Using the numerical results illustrated in Fig. 8b, we evaluate the local wavelength by computing the number of cells between two minima of the A concentration profile. The results are given in Fig. 9 and compared to the analytical prediction deduced from Eqs. (20) and (21). The agreement between the numerical and analytical results is very satisfying in the range $600\leq i<900$. Oscillations of very small amplitude are observed in Fig. 8b in the range $900\leq i<1000$, proving that a very damped Turing structure remains in a small area where instability was predicted. The wavelength of the structure in the range $1\leq i<500$ is slightly affected by the increase of $R_{2}$ from cell $i=500$ but the deviation from the analytical prediction is only $2.5$ percent. This small difference is related to the linear approximation used in wavelength evaluation that neglects non-linear terms that may be more important for large structures. Interestingly, the wavelength is sensitively decreased in the expected area in which the concentration of the reservoir R2 increases: As shown in Fig. 9, the wavelength is reduced from 72 spatial cells to less than 61, before the structure disappears. We conclude that an increase in the concentration of the reservoir R2 related to the inhibitor $B$ is sufficient to account for the destabilization of the Turing pattern associated with a decrease of the wavelength. As anticipated by the results given in Fig. 7, according to which an increase of $R_{2}$ decreases the wavelength $\lambda$ and leads to a negative eigenvalue $\mu_{\max}$ around $(A_{T},B_{T})$, we suggest that an appropriate spatial variation of $R_{2}$ can be used in chemical engineering to stabilize the homogeneous steady state and induce a termination of the Turing pattern in a growing system.

4 Conclusion

In a biomimetic approach, we have addressed the question of the termination of a Turing structure in a growing system. A free boundary is imposed at the growing part, which ensures that the wavelength of the pattern is not perturbed by fixed boundary conditions. After deriving analytical expressions for the stability condition and the wavelength of the structure, we perform a systematic analysis of the effect of all dynamical parameters on the pattern. Apart from the variation of the diffusion coefficient of the activator, a well-chosen variation of the dynamical parameters leads to the desired behavior, i.e. the simultaneous loss of stability and the decrease of the wavelength. In particular, an increase of the effective rate constant $k_{-3}R_{2}$, where $k_{-3}$ is the rate constant of the reaction injecting the inhibitor from the reservoir at the concentration $R_{2}$, is associated with a destabilization of the Turing pattern accompanied by a decrease of the wavelength.

Imposing a spatial variation of the concentration of the reservoir R2 turns out to be an appropriate protocol for chemical engineering. However, the proposed procedure imposes the total length of the structure but not its number of wavelengths. In the framework of developmental biology, for example in the case of the growth of the digits or the spine of the vertebrates, the termination process has to respect the total number of segments for a possible variation in the length of the global structure. Therefore, it is necessary to imagine that the system itself is able to count the number of already formed segments and to trigger the variation of a parameter leading to smaller subsequently formed segments. If the concept of the Turing structure is kept in the formation of biological patterns, the presented results could be used to suggest such relevant parameters. The local increase of the rate constant $k_{-3}$ that would be activated when a given number of segments has already been formed can be straightforwardly proposed. Similarly, the local increase of the rate constant $k_{2}$ controlling the autocatalytic production of the activator or the local decrease of the rate constant $k_{1}$ or $k_{3}$, associated with the absorption of the activator or the inhibitor by reservoirs, would lead to the desired phenomenon. The local decrease of the diffusion coefficient of the inhibitor offers an alternative. The nature of the mechanism that would trigger such a response of the system when a given number of segments has been created remains an open question.

Acknowledgements

This publication is part of a project that has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 711859 and has benefited from financial resources for science awarded by the Polish Ministry of Science and Higher Education in the years 2017-2021 for the implementation of an international cofinanced project.