Upper limits on the robustness of Turing models
and other multiparametric dynamical systems

Deciphering reproducible pattern formation during embryogenesis with multiparameter models calls for mechanisms demonstrating robustness in parameter space. Turing models of diffusion-driven instability [1], which serve as a pivotal category of models, only induce patterns within a highly restricted region of the parameter space [2, 3, 4, 5]. The exploration of such models for increased regions of parameter space and hence enhanced robustness is markedly limited due to their inherent complexity: they incorporate numerous parameters, including rate constants and diffusion constants of activator and inhibitor morphogens, as well as a wave number parameter as part of the analysis.

Here, we consider autonomous spatiotemporal dynamic models that include diffusion, such as the ones describing reaction-diffusion phenomena. Their partial differential equations (PDEs) are written as

where $\bm{X}\in\mathbb{R}^{n}$ are system variables, $\bm{f}$ is an $n$-valued function defined in $n$-dimensional phase space, $\bm{\theta}\in\mathbb{R}^{m}$ is the system-independent parameter vector, $\bm{D}$ is the diffusion matrix, and $\nabla^{2}$ the Laplacian.

Generally, when analyzing a dynamical system, after finding all the fixed points, one needs to study their stability and, eventually, the basins of attraction of all fixed points in the phase space. To elaborate on this procedure, for systems with diffusion-induced instability capable of producing Turing patterns, one requires the system without diffusion to be stable. Hence, one needs to solve the system of equations

for a given $\bm{\theta}$ to find fixed points $\bm{X}^{*}$, and then use the Taylor expansion of the equations around these points to linearise the model to study the stability of the system due to a small perturbation for $t\ll 1$. This procedure is equivalent to writing the $n\times n$ Jacobian matrix $\left.\bm{J}_{0}\right|_{\bm{X}^{*}}$ of the system at $\bm{X}^{*}$ and studying its eigenvalues to classify the fixed point stability. However, finding the eigenvalues of a matrix is computationally expensive $O(n^{3})$, so the above procedure is challenging for high-dimensional phase spaces. The mentioned problem becomes more acute when considering that a candidate dynamical system can have some parameters that change the system’s stability. In this case, we need to redo all the steps for every point in parameter space to study the dynamical system’s characteristics, e.g. the bifurcation diagram.

Specifically, for our reaction-diffusion system with diffusion-induced instability, the Laplacian is removed by spatial Fourier transformation at the expense of an extra parameter, the wave number. Hence, the linear stability analysis must be applied for different wave numbers to specify the dominant wavelength that finally shapes the stationary solution [3]. For a diagonal diffusion matrix $\bm{D}$, the Jacobian for a given wave number, say $k$, is

Hence, any improvement that decreases the procedure’s computational complexity is beneficial for general stability analysis, particularly for pattern formation study.

In this paper, we introduce an efficient $O(n)$ methodology based on the well-known Gershgorin’s theorem [6] in linear algebra. Our method serves as an initial screening process to rule out extensive portions of the parameter space when searching for Turing instabilities. This method, applicable numerically and analytically for simplified models, markedly accelerates the exploration process. We illustrate the efficacy of this approach using three 2-morphogen models – a sigmoidal Hill-function-based model, the Brusselator model [7], and the Lengyel-Epstein [8] model.

Gershgorin’s theorem - We start by introducing Gershgorin’s theorem and consider its geometrical interpretation. As we shall see, it is possible to construct an algorithm that checks the rows or columns of a Jacobian matrix to find those that are unstable or remain stable after introducing diffusion and consequently cannot produce a Turing pattern. The theorem [6] states that for $n\times n$ complex matrix $A=(a_{ij})$ and $r_{i}\equiv\sum_{\begin{subarray}{c}j=1\\ j\neq i\end{subarray}}^{n}|a_{ij}|$ the sum of moduli of off-diagonal elements in the $i$-th row, each union of circles $|z-a_{ii}|\leq r_{i}$ (for $i=1,2,\dots,n$) contains a number of eigenvalues of $A$ equal to the number of circles used to create the union. The analogous result holds if columns of $A$ are considered. Note that in our case, $a_{ij}=\partial_{i}f_{j}$, and for diagonal terms $a_{ii}=\partial_{i}f_{i}$ without and $a_{ii}=\partial_{i}f_{i}-k^{2}D_{i}$ with diffusion.

Let us consider the radius and position of a circle in the complex plane as depicted in Fig. S1. Depending on the sign of the diagonal term $a_{ii}$, the circle’s center is on the real axis’s negative or positive side. Let us say $\lambda_{i}$ is one of the eigenvalues corresponding to the row or column $i$ of the matrix $A$. Four types of circles may appear: Type 1: As it shown in Fig. S1A, we have $a_{ii}<0$ and $h_{i}\equiv|a_{ii}|-r_{i}\geq 0$. Therefore, regardless of where the eigenvalue is inside the circle, its real part must be negative, $Re(\lambda_{i})\leq 0$. Type 2: In Fig. S1B, we can see $a_{ii}<0$, and the center of the circle is placed on the negative side of the real axis. However, since $|a_{ii}|<r_{i}$ and $h_{i}<0$, the real part of $\lambda_{i}$ can be negative or positive. In other words, the theorem is inconclusive about the sign of the real part of the corresponding eigenvalue. Type 3: The diagonal element is positive ($a_{ii}>0$), and the center of the circle is on the positive side of the real axis. Yet, $|a_{ii}|>r_{i}$ and $h_{i}>0$. Thus, the real part of the eigenvalue must be positive $Re(\lambda_{i})>0$.Type 4: Similar to Type 2, the range of Gershgorin’s circle spans from negative to positive values. Therefore we cannot conclusively decide about the real part of the eigenvalue inside this area.

We must emphasize that when two or more Gershgorin’s circles overlap, the eigenvalues lie in the union of the circles. Let us define the lower and upper bounds of each circle as $l_{i}\equiv a_{ii}-r_{i}$ and $u_{i}\equiv a_{ii}+r_{i}$, respectively. For instance, in Fig. S1C, the union of two circles shows that the real part of both eigenvalues corresponding to rows $i$ and $j$ must be in the union of their diameters $Re(\lambda_{i}),Re(\lambda_{j})\in[l_{j},u_{j}]\bigcup[l_{i},u_{i}]=[l_{j},u_{i}]$. And correspondingly, for each row (or column) of the matrix $A$, there exists intervals like

where after taking their unions and reformulating as non-overlapping, disjoint intervals result in

such that $[L_{i},U_{i}]\cup[L_{j},U_{j}]=\emptyset$, for $i\neq j$. Therefore, the $[L_{max},U_{max}]$ is the right-most disjoint interval constructed by the unions of the original intervals, and is sufficient to study it to find the sign of the real part of the largest eigenvalue. This introduces the following possibilities regarding the stability condition: 1. For $U_{max}\leq 0$, the real part of the largest eigenvalue is negative. Consequently, the system is stable. 2. For $L_{max}>0$, the real part of the largest eigenvalue is positive. Consequently, the system is unstable. 3. For $L_{max}\leq 0<U_{max}$ the situation is inconclusive.

At this stage, we can use the obtained results in two different ways. The first possibility is when the Jacobian is written in terms of the model’s parameters, $\bm{\theta}$, and we may be able to derive the right-most disjoint set parametrically. Accordingly, our inequalities define the regions in parameter space where the method can conclusively determine the stability/instability of the system. Nevertheless, the region corresponding to the inconclusive range requires different classification techniques. Note that even if the method is not always conclusive for all the regions in parameter space, it can always find a theoretical lower bound for the volume of the parameter space that the system is stable/unstable.

The second possibility is when studying a system’s linear stability numerically. We propose an algorithm that classifies a given Jacobian matrix into “stable”, “unstable” and “inconclusive” stability groups (see Algorithm (1) in Supplementary Materials - IV A).

Reaction-diffusion models - To study a pattern-forming system given by Eq. (1) with stationary solution $\bm{X}^{*}=(X^{*}_{1},\dots,X^{*}_{n})$, the linear stability of the Jacobian $\left.\bm{J}\right|_{\bm{X}^{*}}$ in Eq. (3) is studied by our proposed method. We can write the lower and upper bounds corresponding to row (or column) $i$ as $l_{i}=\partial_{i}f_{i}-r_{i}$ and $u_{i}=\partial_{i}f_{i}+r_{i}$ for $r_{i}\equiv\sum_{j\neq i}\partial_{j}f_{i}$. Then, the stability/instability criteria of the Jacobian determine the Turing pattern conditions.

It is interesting to see the effect of introducing diffusion. For a given wave number $k$, the inclusion of diffusion shifts all the diagonal terms by $-k^{2}D_{i}$ (see Eq. (3)). Effectively, since diagonal terms correspond to the location of the centers of Gershgorin’s circles, it is geometrically equivalent to saying all the circles shift to the left as shown in Fig. S1D. Note, since the off-diagonal terms have not changed, the circles’ radii remain unchanged.

Indeed, for any given circle partially on the positive side of the real axis, there exists a maximum shift by a wave number defined by $k^{*}_{i}=\sqrt{(r_{i}+\partial_{i}f_{i})/D_{i}}$ that transfers the circle entirely to the negative side of the real axis by $k^{*2}_{i}D_{i}$ (for details see Supplementary Materials - II). Furthermore, the real part of the eigenvalue corresponding to that circle must be negative for all the higher wave numbers $k_{i}>k^{*}_{i}$. By finding $k^{*}_{max}\equiv\max\{k^{*}_{1},\dots,k^{*}_{n}\}$ for all rows (or columns), the linear stability of different wave numbers can be restricted to the range $k$, such that for $k=0$, or the case with no diffusion, $\left.\bm{J}_{0}\right|_{\bm{X}^{*}}$ specifies the stability condition and $\bm{J}(k)$ for $k\in(0,k^{*}_{max}]$ determines the evolution of the dominant wave number in a perturbed system.

Thus, we must have three different regimes: (1) For $\left.\bm{J}_{0}\right|_{\bm{X}^{*}}$, when $U_{max}\leq 0$, all the circles are on the negative side of the real line and including diffusion for any given wave number shifts the circles further to the left. Hence, all real parts of eigenvalues are negative, and diffusion cannot excite any wave number. Consequently, the system is incapable of producing a Turing pattern (“super-stable”). (2) On the contrary, when $L_{max}>0$ for $\left.\bm{J}_{0}\right|_{\bm{X}^{*}}$, the initial stationary state is unstable and incapable of producing a Turing pattern (“unstable”). (3) Finally, when $[L_{max},U_{max}]$ is inconclusive ($L_{max}<0$ and $U_{max}>0$), $\left.\bm{J}_{0}\right|_{\bm{X}^{*}}$ must be studied by finding its eigenvalues. And if one finds that it is stable, the maximum of the dispersion relation $\lambda(k)$ finds the dominant wave number for pattern formation by restricting $k$ to $[0,k^{*}_{max}]$. These regimes are included in our algorithm to speed up the process for checking the possibility of Turing-pattern formation for a given parameter set. Only parameters for which the classification is “inconclusive” need further study of their eigenvalues and are in principle able to form patterns.

Role of diffusing species - The dispersion relation of systems in which all species are diffusing always satisfies $Re(\lambda(k))<0$ for $k>k^{*}_{max}$. In other words, asymptotically, as long as all species are diffusers, we have $\lim_{k\rightarrow\infty}Re(\lambda(k))\rightarrow-\infty$. When some but not all species in a reaction-diffusion model diffuse, introducing the diffusion coefficients into the Jacobian matrix shifts some circles to the left while the others remain in the same place – see Fig. S2. Although one can calculate the $k^{*}_{max}$ values for the diffuser rows in the matrix, a special situation can arise for $k>k^{*}_{max}$ in the dispersion relation. For instance, consider two among three species are diffusers as shown in Fig. S2. As we can see, after shifting the diffusers circle by an amount $-k^{*2}_{max}D_{1}$ and $-k^{*2}_{max}D_{2}$ respectively, the third eigenvalue corresponds to the non-diffuser element remains positive for all $k>k^{*}_{max}$. Consequently, the real part of the dispersion relation can remain positive with no upper bounds. As a result, no dominant wave number exists to create a stationary pattern. Fortunately, it is easy to find these cases algorithmically without calculating the eigenvalues.

Tightening the bounds - Given an invertible matrix $\bm{D}$, $\bm{B}=\bm{D}\bm{A}\bm{D}^{-1}$ introduces an equivalence relation between square matrices $\bm{A}$ and $\bm{B}$ such that matrix $\bm{B}$ has the same eigenvalues as $\bm{A}$ (see Supplementary Material - III for details). Defining an invertible $n\times n$ diagonal matrix $\bm{D}$ as the identity matrix with exception of matrix element $D_{ii}=1/d_{i}$, the transformed matrix $\bm{D}A\bm{D}^{-1}$ has the form

The transformation’s effect is similar to dividing all the elements of the row $i$ by $d_{i}$ and multiplying the elements of the column $i$ by $d_{i}$. Consequently, the diagonal term $a_{ii}$ remains the same.

Consider the rows of the resulting matrix. The radii of all circles corresponding to rows other than $i$ expand by the amount $|a_{ji}|(d_{i}-1)$ (for $d_{i}>1$), and the $i$-th radius shrinks by the factor of $1/d_{i}$, while the centers of all circles stay the same. Since the eigenvalues of the transformed matrix are the same as the original one, one can hope the shrunk circle becomes isolated from the rest since the expansions of the other radii are smaller than the shrinking of the single radius. In practice, one can find $d_{i}$ that isolates the circle with the largest center from the rest. The interval of all rows except the $i$-th is $[l_{j},u_{j}]=\left[a_{jj}-|a_{ji}|(d_{i}-1)-r_{j},\quad a_{jj}+|a_{ji}|(d_{i}-1)+r_{j}\right]$, and the interval for row $i$ is $[l_{i},u_{i}]=\left[a_{ii}-\frac{r_{i}}{d_{i}},\quad a_{ii}+\frac{r_{i}}{d_{i}}\right]$.

To tighten the bounds, we have two distinct cases that can be studied separately. Case (1): When the diagonal term is positive, or $a_{ii}>0$, to isolate the circle corresponding to row $i$, its leftmost point, or $l_{i}=a_{ii}-\frac{r_{i}}{d_{i}}$, must be larger than every other circle’s rightmost point, or $u_{j}=a_{jj}+|a_{ji}|(d_{i}-1)+r_{j}$. As explained in detail in Supplementary Materials - III, if all the rows $j$ and the largest one, $i$, satisfy the inequalities

simultaneously, there exists a $d_{i}$ that isolates the rightmost eigenvalue, and the Jacobian is conclusively unstable. Case (2): When the largest diagonal term is negative, or $a_{ii}<0$, the latter implies all the other diagonal terms are negative too. In this case, we search for a possible shrinkage value $d_{i}$, such that while the circle of the row $i$ shrinks with its upper bound on the negative side of the real axis, the growth of all the other circles keeps them at the negative side. Combined, this leads to two conditioned bounds as

Hence, if a non-empty interval can be found that satisfies the above inequalities, the Jacobian is conclusively “super-stable” (see Algorithm 2 in Supplementary Materials - IV B).

Conclusion - We have introduced an efficient $O(n)$ method that uses Gershgorin’s theorem to define robustness bounds in dynamical systems, particularly beneficial for reaction-diffusion models such as diffusion-driven Turing systems [3]. This approach can eliminate unstable non-diffusive cases or solutions that remain stable post-diffusion while also setting an upper limit for wave numbers capable of pattern formation. When applied to specific models in the appendix, our method does not only enhance numerical algorithms’ speed, but it also facilitates an analytical study of parameter space. The former capability is exemplified by a Hill-function-based Turing model which leads to a rejection of $99.3\%$ of the parameter combinations in table 1), while the latter is demonstrated for the classic Brusselator model in Fig. S3 (and by the Lengyel-Epstein model in Supplementary Materials - V B). The method’s utility increases significantly when accounting for parameters that alter the behavior of a potential dynamical system. For example, generating bifurcation profiles of dynamical systems necessitates repeated linear stability analysis for many different parameter values, a challenging endeavor in high dimensional phase spaces.

Appendix on applications to Turing models — To test the fraction of rejections and, consequently, the speed up due to the Algorithms (1) and (2), we use a biologically inspired model capable of producing Turing patterns [3]. The reaction-diffusion PDE with nine free parameters is written as

Note that the nonlinear terms are Hill functions that regulate activation and inhibition of molecule production, e.g., in gene expression. The Jacobian of the linearized form of the above equations is a two-by-two matrix, and in practice, the computational cost of calculating its eigenvalues is not much different than our algorithm. However, we selected this model since the algorithm’s correctness can be easily checked by comparing the determinant and trace of the Jacobian.

In this simulation, we selected one billion parameter combinations and applied Algorithm (1) to classify them into “unstable”, “super-stable”, “inconclusive” and “no fixed point”. Note that the case “no fixed point” refers to the parameter combinations for which the root-finding algorithm could not find any stationary solutions. We first used our Algorithm (1) for a row-wise comparison, and after that, by using the inconclusive results from the first run, we used the algorithm again for a column-wise calculation. And finally, we classified the remaining inconclusive cases using our Algorithm (2). These results are presented in table 1, showing that more than $99.3\%$ of the parameter combinations were rejected. This provides an upper limit on the robustness of Turing patterns given a certain sampling of parameter space [3].

Brusselator model - Next, we look for inequalities that separate the parameter space of the Brusselator model into “inconclusive“ or otherwise. The Brusselator is a two-species reaction-diffusion model with a set of PDEs given by

for two parameters $a,b>0$. Using the stability analysis, we can derive the Jacobian of the model at its fixed point $(u^{*},v^{*})=(a,b/a)$ and the corresponding rows and columns intervals as

Note that since the intervals depend on the parameters, it is easier to check the sign of the circles’ centers and classify the conditions than directly taking the union.

For the first row, the circle’s center is at $b-1$; therefore, for $0<b\leq 1$, it is on the negative side of the real axis. At the same time, the center of the second row’s circle is always on the negative side ($-a^{2}$). Consequently, the region of parameter space in which the Jacobian is stable before and after the inclusion of diffusion derives as

for $a,b>0$. At the same time, for $b>1$, the center of the first row’s circle is on the positive side of the real axis. Thus, the Jacobian is unstable when the first interval’s lower bound becomes positive, leading to

The conditions from both Eqs. (12) and (13) are shown in Fig. S3A. Similarly, for the column’s intervals, if $b<1/2$, both centers are on the negative sides of the real axis, and the condition of stability writes as $b\leq 1/2$ for $a,b>0$. However, for $b>1/2$, the interval of the first column stays between $-1$ and $2b-1$, which means the case is inconclusive. Fig. S3B depicts the columns’ results.

Supplementary Materials:
Upper limits on the robustness of Turing models
and other multiparametric dynamical systems

In these supplementary materials, we provide details on the results presented in the main text. We also present our Algorithms (1) and (2), and additional results on the Lengyel-Epstein model.