A Spatial Filtering Approach to Biological Patterning

TL;DR

This paper introduces a discrete spatial filtering framework to model how cellular interactions influence pattern formation in biological systems, providing insights into robustness and patterning mechanisms.

Contribution

It presents a novel filtering approach to predict biological patterning, linking intercellular interactions to spatial frequency modulation, with applications to developmental biology.

Findings

Patterns can form without instabilities, similar to Turing patterns.

Biological systems show robustness to environmental and cellular noise.

Spatial filtering simplifies understanding of pattern formation processes.

Abstract

Interactions between neighboring cells are essential for generating or refining patterns in a number of biological systems. We propose a discrete filtering approach to predict how networks of cells modulate spatially varying input signals to produce more complicated or precise output signals. The interconnections between cells determine the set of spatial modes that are amplified or suppressed based on the coupling and internal dynamics of each cell, analogously to the way a traditional digital filter modifies the frequency components of a discrete signal. We apply the framework to two systems in developmental biology: the Notch-Delta interaction that shapes \textit{Drosophila} wing veins and the Sox9/Bmp/Wnt network responsible for digit formation in vertebrate limbs. The latter case study demonstrates that Turing-like patterns may occur even in the absence of instabilities. Results also indicate that developmental biological systems may be inherently robust to both correlated and uncorrelated noise sources. Our work shows that a spatial frequency-based interpretation simplifies the process of predicting patterning in living organisms when both environmental influences and intercellular interactions are present.

Figures

Tables

Table 1. Table 1: Parameters used in the Notch-Delta model simulations, unless noted otherwise in the text (Figures 4 to 6 , 19 , and 20 ).

Parameter	Value	Description	Source
${\alpha }_N$	10	“leakiness” of Notch expression (RFU/hr)	[23] Table S1 (Figure S4A)
${\beta }_{D0}$	17.5	max. Delta production rate (RFU/hr)	[8] Table S3 (Figure 4C)
$x_0$	7	number of cell diameters	[8] Table S3 (Figure 4C)
${\beta }_{D_i}$	${\beta }_{D0}{e}^{-\left|i\right|/x_0}$	Delta production rate (RFU/hr) for cell $i$	[8] (Figure 4C)
${\overline{\underset{¯}{\beta }}}_D$	9.09	Delta production rate (RFU/hr) for linearization	$\frac{1}{N}{\sum }_{i=0}^{N-1}{\beta }_{D_i}$
${\beta }_N$	10	Notch production rate (RFU/hr)	[8] Table S3 (Figure 4C)
${\beta }_R$	150	reporter production rate (RFU/hr)	[8] Table S3 (Figure 4C)
$\gamma$	0.1	Notch, Delta decay rate (1/hr)	[8] Table S3 (Figure 4C)
${\gamma }_R$	0.05	reporter decay rate (1/hr)	[8] Table S3 (Figure 4C)
$k_c$	0.25	inverse cis-interaction strength	[8] Table S3 (Figure 4C)
$k_t$	5	inverse trans-interaction strength	[8] Table S3 (Figure 4C)
$n$	2	Hill coefficient for Notch-Delta activation of reporter	-
$m$	2	Hill coefficient for reporter repression of Delta	-
$k_{RS}$	300,000	affinity of reporter induction	[23] Table S1 (Figure S4A)
$k_{NS}$	$5\times {10}^7$	affinity of reporter induction	-

Table 2. Table 2: Parameters used in the simulations of digit formation, unless noted otherwise in the text (Figures 10 , 8 , 9 ). Values are from Table ST4 and text of [ 26 ] .

Parameter	Value	Description
${\alpha }_{sox9}$	0	constitutive Sox9 production rate
${\alpha }_{bmp}$	16.9	constitutive Bmp production rate
${\alpha }_{wnt}$	13.7	constitutive Wnt production rate
$k_2$	1	Bmp promotion of sox9 expression
$k_3$	1	Wnt repression of sox9 expression
$k_4$	1.59	Sox9 repression of bmp expression
$k_5$	0.1	Bmp decay rate
$k_7$	1.27	Sox9 repression of wnt expression
$k_9$	0.1	Wnt decay rate
$d_b$	2.5	diffusion coefficient for Bmp
$d_w$	1	diffusion coefficient for Wnt
$l$	1.7	distance between cells

Table 3. Table 3: Parameters used in the simulations of digit formation with a morphogen gradient, unless noted otherwise in the text (Figures 12 , 21 to 30 ). Values are from Methods in [ 27 ] (Figures 4 and 5).

Parameter	Value	Description
${\alpha }_{sox9}$	0	constitutive Sox9 production rate
${\alpha }_{bmp}$	0.1	constitutive Bmp production rate
${\alpha }_{wnt}$	1.2	constitutive Wnt production rate
${\mu }_F$	0.1	Fgf decay rate
$k_2$	1	Bmp promotion of sox9 expression
$k_3$	3	Wnt repression of sox9 expression
$k_4$	6	Sox9 repression of bmp expression
$k_5$	0.1	Bmp decay rate
$k_7$	2.4	Sox9 repression of wnt expression
$k_9$	0.1	Wnt decay rate
$k_f$	$\frac{2}{3}$	strength of Fgf influence on $k_4$, $k_7$
$d_b$	160	diffusion coefficient for Bmp
$d_w$	25	diffusion coefficient for Wnt
$d_F$	600	diffusion coefficient for Fgf
$l$	4	distance between cells

Equations

$$\begin{cases}\dot{x_{i}}(t)=f(x_{i}(t),v_{i}(t),u_{i})\\ w_{i}(t)=g(x_{i}(t))\\ y_{i}(t)=h(x_{i}(t))\\ v(t)=\left(M\otimes I_{q}\right)w(t)\ \end{cases}\$$

$$A:=\frac{\partial f}{\partial x_{i}}\bigg{|}_{(\underline{\bar{x}}^{*},\underline{\bar{v}}^{*},\underline{\bar{u}})},~{}B_{v}:=\frac{\partial f}{\partial v_{i}}\bigg{|}_{(\underline{\bar{x}}^{*},\underline{\bar{v}}^{*},\underline{\bar{u}})},~{}B_{u}:=\frac{\partial f}{\partial u_{i}}\bigg{|}_{(\underline{\bar{x}}^{*},\underline{\bar{v}}^{*},\underline{\bar{u}})},~{}C:=\frac{dh}{dy_{i}}\bigg{|}_{\underline{\bar{x}}^{*}},G:=\frac{dg}{dx_{i}}\bigg{|}_{\underline{\bar{x}}^{*}}.$$

$$\hat{u}=T^{-1}(u-\underline{u})$$

$${\left[S\right]}_{kk}=-C(A+{\lambda }_k(M)B_{v}G{)}^{-1}{B}_u$$

$$d\hat{x}_{k}(t)=\left[\left(A\otimes I_{N}\right)+\left(B_{v}G\right)\otimes\Lambda\right]\hat{x}_{k}(t)dt+\left(B_{u}\otimes I_{N}\right)d\hat{u}_{k}(t).\$$

$$\mathbb{E}\left[\left|\hat{y}_{k}\right|^{2}\right]=\mathbb{E}\left[Tr\left(\hat{y}_{k}\hat{y}_{k}^{T}\right)\right]=Tr\left(C\mathbb{E}\left[\hat{x}_{k}\hat{x}_{k}^{T}\right]C^{T}\right)=Tr\left(CQ_{k}C^{T}\right),\$$

$$\mathbb{E}\left[Tr\left({\hat{y}}_k{\hat{y}}_k^T\right)\right]$$

$$={\int }_0^{\infty }Tr\left(G_k(t)\mathbb{E}\left[d{\hat{u}}_{k}d{\hat{u}}_k^T\right]G_k(t{)}^T\right)dt={\int }_0^{\infty }Tr\left(G_k(t)G_k(t{)}^T\right)dt$$

$$\displaystyle=:||G_{k}(t)||^{2}_{\mathcal{H}_{2}},\$$

$$\left(A+\lambda_{k}B_{v}G\right)Q_{k}+Q_{k}\left(A+\lambda_{k}B_{v}G\right)^{T}+B_{u}B_{u}^{T}=0.\$$

$$\begin{cases}\dot{m}_{i}=-\gamma_{m}m+\alpha_{m}f(v_{i},u_{i},p_{i})\\ \dot{p}_{i}=-\gamma_{p}p+\alpha_{p}m\\ y=p\\ v=Mp\end{cases}\$$

$$\begin{cases}\dot{\tilde{m}}_{i}=-\gamma_{m}\tilde{m}_{i}+\alpha_{m}\left(F_{v}\tilde{v}_{i}+F_{u}\tilde{u}_{i}+F_{p}\tilde{p}_{i}\right)\\ \dot{\tilde{p}}_{i}=-\gamma_{p}\tilde{p}_{i}+\alpha_{p}\tilde{m}_{i}\\ \tilde{y}=\tilde{p}\\ \tilde{v}=M\tilde{p}\end{cases}\$$

$$\left[S\right]_{kk}=\frac{\frac{\alpha}{\gamma}F_{u}}{1-\frac{\alpha}{\gamma}\left(F_{p}+F_{v}\lambda_{k}\left(M\right)\right)}\$$

$$\frac{\alpha}{\gamma}\left(F_{v}\lambda_{k}\left(M\right)+F_{p}\right)<1.\$$

$$\tilde{u}=\left(u-\underline{u}\right)=T\hat{u}=\sum_{k=0}^{N-1}\hat{u}_{k}t_{k}.\$$

$$\tilde{y}^{*}=\sum_{k=0}^{N-1}\left[S\right]_{kk}\hat{u}_{k}t_{k},\$$

$$M=\frac{1}{2}\left[\begin{matrix}0&1&0&{1}&0&1\\ 1&0&1&{1}&0&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 1&0&0&{1}&1&0\end{matrix}\right].\$$

$$\begin{cases}\dot{N_{i}}(t)=\beta_{N}-\gamma N_{i}(t)-\frac{N_{i}(t)v_{D_{i}}(t)}{k_{t}}-\frac{N_{i}(t)D_{i}(t)}{k_{c}}\\ \dot{D_{i}}(t)=\underline{\bar{\beta}}_{D}+u_{i}-\gamma D_{i}(t)-\frac{D_{i}(t)v_{N_{i}}(t)}{k_{t}}-\frac{N_{i}(t)D_{i}(t)}{k_{c}}\\ \dot{R}_{i}(t)=\beta_{R}\frac{\left(N_{i}(t)v_{D_{i}}(t)\right)^{n}}{k_{RS}+\left(N_{i}(t)v_{D_{i}}(t)\right)^{n}}-\gamma_{R}R_{i}(t)\\ y_{i}(t)=Cx_{i}(t)\\ w_{i}(t)=\left[\begin{matrix}N_{i}(t)\\ D_{i}(t)\end{matrix}\right]\\ v(t)=\left(M\otimes I_{2}\right)w(t)\ \end{cases}\$$

$$A$$

$$B_v$$

$$b_{1}:=\beta_{R}nk_{RS}\frac{\underline{\bar{v}}_{D}^{*n}N^{*n-1}}{\left(k_{RS}+\left(\underline{\bar{N}}^{*}\underline{\bar{v}}_{D}^{*}\right)^{n}\right)^{2}},~{}b_{2}:=\frac{\underline{\bar{N}}^{*}}{\underline{\bar{v}}_{D}^{*}}b_{1}\$$

$$B_{u}^{corr}=\left[\begin{matrix}1\\ 1\\ 0\end{matrix}\right],~{}B_{u}^{uncorr}=\left[\begin{matrix}1&0\\ 0&1\\ 0&0\end{matrix}\right].\$$

$$\begin{cases}\dot{s}_{i}(t)=\alpha_{s}+u_{i}+k_{2}b_{i}(t)-k_{3}z_{i}(t)-(s_{i}(t)-s_{0})^{3}\\ \dot{b}(t)=\alpha_{bmp}+u_{i}-k_{4}s_{i}(t)-k_{5}b_{i}(t)+\frac{d_{b}}{l^{2}}v_{b_{i}}(t)\\ \dot{z}_{i}(t)=\alpha_{wnt}+u_{i}-k_{7}s_{i}(t)-k_{9}z_{i}(t)+\frac{d_{z}}{l^{2}}v_{z_{i}}(t)\\ w_{i}(t)=\left[\begin{matrix}b_{i}(t)\\ z_{i}(t)\end{matrix}\right]\\ v(t)=(M\otimes I_{2})w(t)\ \end{cases}\$$

$$A$$

$$B_v$$

$$\begin{array}[]{c|c|c|c|c}0&N_{R}&2N_{R}&\dots&(N_{C}-1)N_{R}\\ \hline\cr 1&N_{R}+1&2N_{R}+1&\dots&(N_{C}-1)N_{R}+1\\ \hline\cr\vdots&\vdots&\vdots&\ddots&\vdots\\ \hline\cr N_{R}-2&2N_{R}-2&3R-2&\dots&N_{C}N_{R}-2\\ \hline\cr N_{R}-1&2N_{R}-1&3N_{R}-1&\dots&N_{C}N_{R}-1\end{array}.\$$

$$M:=\left(M_{R}\otimes I_{N_{R}}\right)+\left(I_{N_{C}}\otimes M_{C}\right).\$$

$$T:=\left(T_{R}\otimes I_{N_{R}}\right)\left(I_{N_{C}}\otimes T_{C}\right)=T_{R}\otimes T_{C},\$$

$$\lambda_{m+nN_{R}}\left(M\right)=\lambda_{m}\left(M_{C}\right)+\lambda_{n}\left(M_{R}\right)\$$

$$\tilde{Y}=T_{C}\hat{Y}T_{R}^{T}=T_{C}\left(\Lambda_{S}\odot\hat{U}\right)T_{R}^{T}\$$

$${\left[{\Lambda }_S\right]}_{mn}=-C{\left[A+\left({\lambda }_m(M_C)+{\lambda }_n(M_R)\right)B_{v}G\right]}^{-1}{B}_u.$$

$$\begin{cases}\dot{F}_{i}(t)=\alpha_{F}+u_{1i}-\mu_{F}F_{i}(t)+\frac{d_{f}}{l^{2}}v_{F_{i}}(t)\\ \dot{s}_{i}(t)=\alpha_{s}+u_{2i}+k_{2}b_{i}(t)-k_{3}z_{i}(t)-s_{i}(t)^{3}\\ \dot{b}(t)=\alpha_{bmp}+u_{2i}-k_{4}\left(1-k_{f}F_{i}(t)\right)s_{i}(t)-k_{5}b_{i}(t)+\frac{d_{b}}{l^{2}}v_{b_{i}}(t)\\ \dot{z}_{i}(t)=\alpha_{wnt}+u_{2i}-k_{7}k_{f}F_{i}(t)s_{i}(t)-k_{9}z_{i}(t)+\frac{d_{z}}{l^{2}}v_{z_{i}}(t)\\ \dot{w}_{i}(t)=\left[\begin{matrix}F_{i}(t)\\ b_{i}(t)\\ z_{i}(t)\end{matrix}\right]\\ v(t)=\left(M\otimes I_{3}\right)w(t)\ \end{cases}\$$

Full text

1 INTRODUCTION

Biological organisms rely on spatial variation in cell activity to coordinate diverse phenomena including contrast enhancement in the visual system [2] and body planning in developing embryos [3]. Interactions between neighboring cells play a crucial role in generating spatial patterns spontaneously from stochastic initial conditions or by refining simple inputs, such as chemical concentration gradients, into complex outputs, such as stripes in gene expression [4], [5]. Mathematical theory in developmental biology has emphasized spontaneous pattern formation through the reaction-diffusion (Turing) mechanism [6] as well as contact- or diffusion-mediated lateral inhibition [7], [3], [8]. In practice, however, the conditions necessary for spontaneous pattern formation may be prohibitively difficult to satisfy.

Prepattern processing—also known as Wolpert’s theory of positional information [9]—is an attractive and flexible alternative to spontaneous patterning, but mathematical analysis of prepattern processing has been largely limited to numerical simulations (e.g, [8]). Prepatterns may arise directly from environmental influences that differ by cell or from consistent, preinduced parameter variation across space.

We propose a discrete filtering approach to analyze how networks of interacting cells respond to prepatterns. The framework elucidates which components of spatial structure are amplified and which ones attenuated by the system to produce an output from any given input. The insights gained from this perspective challenge the conventional notion that instability is necessary for complex patterning; for example, our approach reveals that Turing-like stripes can emerge from a stable system lacking diffusion-driven instability, and furthermore that external noise reinforces rather than combats this behavior (Section 4).

In Section 2 we present the setup for the framework. We combine the internal dynamics of cell behavior with interaction between cells by modeling each cell as an input-output module coupled to other modules. We examine the steady-state gains for constant-in-time, spatially varying inputs (prepatterns) and show that the system behaves as a discrete spatial filter, where the interconnectivity between cells determines the spatial modes, while the coupling and input-output dynamics dictate how each mode is scaled to generate a readout pattern. We also examine the system response to temporal and spatially varying noise inputs, measured with the $\mathcal{H}_{2}$ norm, to determine which spatial modes are sensitive to stochastic influence. Lastly, we show how to apply the approach by considering a simple model of gene expression that exemplifies two of the most common classes of filter behaviors—highpass or lowpass—depending on the choice of parameters.

In the remaining three sections we demonstrate the utility of the filtering perspective by examining two biological case studies: the Notch-Delta system in developing fruit fly wings (Section 3) and the Sox9/Wnt/Bmp network in vertebrate digit formation (Sections 4 and 5). We conclude with a brief summary and areas for future research.

2 THE SPATIAL FILTERING APPROACH

The main contribution of this paper is a filtering perspective for analyzing prepattern processing in developmental biological systems. A central component of our approach is spatial mode decomposition, a common tool in distributed systems analysis (e.g., [10]) that has previously been applied to detect instabilities in cellular networks lacking external inputs [11]. In this section we introduce generalized notation followed by a derivation of the filter coefficients and the noise amplification factors that we will use throughout the remainder of the paper.

2.1 Notational Conventions

We use the following notational conventions (see also Supplementary Figure 14):

• •

Cells are indexed by $i$ in vector form and spatial modes are indexed by $k$ in vector form or $(m,n)$ in an array, unless noted otherwise.

• •

Inputs except white noise in the context of the $\mathcal{H}_{2}$ norm are assumed constant in time.

• •

Vectors containing strictly constant-in-space entries are designated with an underline. The entries corresponding to any fixed point in space are additionally labeled with an overbar, e.g., $\underline{u}=\underline{\bar{u}}\mathds{1}_{N}$ where $\underline{\bar{u}}\in\mathbb{R}$ and $\mathds{1}_{N}$ is the length $N$ vector of all ones.

• •

Steady-state values for time-dependent variables are designated with superscript asterisks. Constant-in-space steady-state (i.e., homogeneous) solutions to nonlinear systems are designated with both an asterisk and an underline, e.g., $\underline{y}^{*}=\underline{\bar{y}}^{*}\mathds{1}_{N}$.

• •

“Actual” values in the standard basis are unadorned. Perturbations from constant-in-space values are designated with a tilde; time-dependent perturbations are understood to be linear approximations of “actual” nonlinear solutions, e.g., $\tilde{x}_{i}(t)\approx x_{i}(t)-\underline{\bar{x}}^{*}$. Perturbed variables in the basis $T$ are designated with a hat, e.g., $\hat{x}^{*}=T^{-1}\tilde{x}^{*}$.

2.2 System Dynamics and Filter Coefficients

We consider a generalized system of $N$ identical cells with the state variables of the $i$th cell at time $t$ given by $x_{i}(t)\in\mathbb{R}^{n}$, readout $y_{i}(t)\in\mathbb{R}$, and constant-in-time input $u_{i}\in\mathbb{R}$, which may represent an environmental stimulus or intrinsic parameter variation. Coupling occurs via $v_{i}(t)\in\mathbb{R}^{q}$ and output $w_{i}(t)\in\mathbb{R}^{q}$ where $q\leq n$. Let the vectors for the full system be the vertical concatenation $x(t)^{T}:=[x_{0}(t)^{T}x_{1}(t)^{T}...x_{N-1}(t)^{T}]$ and similarly for $u$, $y(t)$, $w(t)$, and $v(t)$. The dynamics of the $i$th cell and the full linear coupling between the $N$ cells are given by

[TABLE]

where $\otimes$ is the Kronecker product, $I_{q}$ is the $q\times q$ identity matrix, and $M\in\mathbb{R}^{N\times N}$.

The system (1) accommodates a wide range of specific deterministic models. Intercellular processes such as gene expression and protein decay are encapsulated by appropriate definition of the evolution function $f$ for chemical concentrations $x_{i}$, including the effect of environmental stimuli or parameter values $u_{i}$ as well as signals from neighbors $v_{i}$. The output $w_{i}$ is the subset of elements in $x_{i}$ that transmit signals to neighbors, with the method of transmission (e.g., diffusion, cell-to-cell contact) and the neighboring cells specified by the interconnection matrix $M$. The readout $y_{i}$ isolates a quantity of interest to the user, which may be experimentally measurable (e.g., fluorescence) or simply relevant to a particular model (see examples in Sections 3, 4, and 5).

The vectors indexed by $i$ describe patterns by the concentration of chemicals within individual cells at individual points in space. A full pattern is reconstructed from $N$ elements, each of which represents the concentration of a chemical in a single cell. However, patterns can also be thought of as combinations of spatially varying components that span multiple cells, e.g., stripes of varying thickness (frequency). When weighted and summed, these spatial modes can represent arbitrary patterns of interest. We use the term “filtering” to refer to the process by which a network of interacting cells alters the weighting of the spatial modes of the input, thereby producing a readout that is built from the same components as, but differs in appearance from, the input. A key approximation to facilitate the analysis is that coupling between spatial modes is negligible, such that the readout can be expressed as a linear sum of the same set of spatial modes used to represent the input. In analogy to traditional signal processing, the network of cells plays the role of a linear time-invariant system (filter) that modifies the frequency components of a (spatially) varying signal. The following proposition formalizes this concept mathematically.

Proposition 1**.**

If the system described by (1) satisfies

$M\mathds{1}_{N}=\mu\mathds{1}_{N}$* and $M$ is diagonalized by $T$ ($M=T\Lambda T^{-1}$),* 2. 2.

given $\underline{\bar{u}}\in\mathbb{R}$, $\exists~{}\underline{\bar{x}}^{*}\in\mathbb{R}^{n}$ such that $f\left(\underline{\bar{x}}^{*},\mu g\left(\underline{\bar{x}}^{*}\right),\underline{\bar{u}}\right)=0$ and $\underline{x}^{*}:=\mathds{1}_{N}\otimes\underline{\bar{x}}^{*}$, $\underline{u}:=\underline{\bar{u}}\mathds{1}_{N}$, 3. 3.

the homogeneous steady state $(\underline{x}^{*},\underline{u})$ is stable,

then the system may be linearized about $(\underline{x}^{*},\underline{u})$ with linearization matrices

[TABLE]

A constant-in-time, varying-in-space input $u$ to (1) equates to a perturbing input

[TABLE]

to the linearized system in the coordinate system $T$. Then the steady-state perturbed readout $\tilde{y}^{*}$ in the basis $T$ is $\hat{y}^{*}:=T^{-1}\tilde{y}^{*}=S\hat{u}$ where $S$ is a diagonal matrix and

[TABLE]

is the steady-state gain of the $k$th eigenvector of $T$.

The conditions (1) through (3) ensure that the network, when given a constant-in-space input, will admit a stable, homogeneous steady-state solution, and that the expression pattern across cells can be represented in a complete orthonormal basis other than the standard; this basis $T$ comprises the modes. The $\left[S\right]_{kk}$ collectively form the “filter coefficients”, which dictate how the corresponding $N$ spatial modes are multiplicatively scaled by the system when the input is no longer constant in space (Figure 2). In other words, the matrix $S$ “filters” the perturbed input into a perturbed readout with respect to the eigenvectors, or spatial modes, of $M$ as contained in $T$. In contrast to the conditions for spontaneous pattern formation, our approach does not require the perturbed system to be unstable; large amplification of spatial modes is possible even when the system is stable. Figure 3 shows eight examples of prototypical filter behaviors that vary with interaction type and cellular interconnectivity.

Many continuous pattern-forming and distributed dynamical systems exhibit spatial invariance of the dynamics with respect to linear transformations such as reflections, rotations, or translations [12], [10]. The discrete-space cellular network has a direct analog: If $M$ is invariant under a linear transformation, then $S$ is also invariant under the same transformation, since the system dynamics are identical within each cell and therefore the only spatial information contained within the system is contained in $M$. Formally, if we let $\Pi\in\mathbb{R}^{N\times N}$ be a linear transformation and $M$ and $\Pi$ commute, then $\Pi$ and the filter coefficient matrix $S$ also commute (i.e., the map from input to output is equivariant). Thus, $M$’s permutation $M_{\Pi}:=\Pi M\Pi^{-1}$ shares the same eigenvectors $T$ and corresponding eigenvalues $\Lambda$ as $M$, which implies that the filter coefficients for a system with interconnection matrix $M$ are the same as for that system with interconnection matrix $M_{\Pi}$.

2.3 Stochastic Influence on Patterning

The role of stochastic influences in biological patterning is a subject of ongoing theoretical and experimental interest (e.g., [13], [14]). Here, we concern ourselves with the response of spatial modes to time-varying white noise inputs, for which the $\mathcal{H}_{2}$ norm of the system quantifies the expected power of the perturbed readout. The $\mathcal{H}_{2}$ norm has previously been used to analyze energy amplification in channel flows [15], networks of cells [16], and reaction-diffusion systems [17], among others.

To begin our analysis we rewrite the linearized ordinary differential equations in the form of a nonlinear Langevin equation (Itô interpretation). Since the $N$ modes are decoupled we can write the equation for the perturbed states in the $k$th mode as

[TABLE]

Here $\hat{u}(t)$ is an $n_{u}$-dimensional independent standard Wiener process, also known as the standard Brownian motion process. Implicitly we assume that concentrations of reactants are high enough to permit us to neglect molecular-level fluctuations, which cannot accurately be described by the Langevin approach [18].

With slight abuse of notation, $d\hat{u}_{k}(t)$ is stationary, therefore the variance of the readout $y(t)=Cx(t)$ in mode $k$ does not change in time. The variance is given by

[TABLE]

where $Q_{k}:=\mathbb{E}\left[\hat{x}_{k}\hat{x}_{k}^{T}\right]$ is the covariance of the reactants in the $k$th mode.

Let $G_{k}(t)$ be the impulse response of (2) for readout $y(t)$. We could equivalently write

[TABLE]

from which we deduce that the $\mathcal{H}_{2}$ norm is equivalent to the variance of $\hat{y}_{k}$ and can be calculated as $Tr\left(CQ_{k}C^{T}\right)$ where $Q_{k}$ is the positive semi-definite solution to the Lyapunov equation

[TABLE]

The unit variance of $d\hat{u}_{k}(t)$ allows us to interpret $||G_{k}(t)||^{2}_{\mathcal{H}_{2}}$ as the ratio of the variance of the readout to the variance of the input in mode $k$. Moreover, since $d\hat{u}_{k}(t)$ is zero mean, the squared $\mathcal{H}_{2}$ norm is also equivalent to the time integral of the expected power spectral density, or the factor by which the system amplifies the average power of the readout within mode $k$. Those modes with the highest $\mathcal{H}_{2}$ norms are most strongly amplified by the external noise source.

2.4 Constructing the Interconnection Matrix

In the remainder of this paper we construct the interconnection matrix $M$ for a particular signal as follows:

The length $N$ vector of all ones $\mathds{1}_{N}$ is an eigenvector of $M$, which implies that a homogeneous steady-state solution exists. 2. 2.

The $i$th, $j$th entry $\left[M\right]_{ij}$ for $i\neq j$ is 0 if cell $i$ is not connected to cell $j$. Otherwise $0<\left[M\right]_{ij}$, where the magnitude $\left[M\right]_{ij}$ captures the “strength” of the connection. 3. 3.

The diagonal entries $\left[M\right]_{ii}$ encapsulate the “signaling cost” associated with interaction. Negative values imply the cell loses signal to transmit to its neighbors, e.g., diffusion.

In many biological systems, cells can be approximated to have the same distance between them and the same communication strength with each of their neighbors. In such systems, the corresponding spatial modes are sinusoidal, giving rise to stripes or spots. Lower-frequency modes correspond to longer-wavelength spatial modes, while higher-frequency modes correspond to shorter-wavelength spatial modes. The relationship between patterning wavelength and spatial mode frequency enables these systems to be interpreted from the standpoint of how the weights of the frequency components in an input are scaled to produce the readout, analogous to filtering as it is understood in discrete signal processing. In this paper we will consider basis vectors arising from a line or sheet of regularly spaced cells with periodic or no-flux boundary conditions. The modes then pertain to two common signal processing transforms: the discrete Fourier transform (DFT) for periodic boundaries or the second discrete cosine transform (DCT-2) for no-flux boundaries. The eigenvectors and eigenvalues for these transforms are well known (e.g., [19]); a review is offered in Supplementary Section 8. We will assume modes are indexed in order of increasing frequency with increasing $k$ toward $\frac{N}{2}$ for the DFT and $N$ for the DCT-2.

2.5 Minimal Model: Gene Expression with Autoregulation

The following example is a simple model that is easy to solve analytically for the filter coefficients. We begin with a brief description of gene expression for readers who may not be familiar with the biology, including terminology that will be used in later sections. We then apply the filtering approach to the example, including an expansion of the matrix notation to emphasize the role of the filter coefficients as “weights” for the spatial modes. As this model focuses on biological and filtering concepts, intercellular interaction is described only in the most general terms, leaving exploration of the underlying mechanisms to later examples.

The case studies in this paper deal with gene expression, or the process by which a gene coded in DNA is transcribed into mRNA molecules that are then translated into protein molecules (Supplementary Figure 13). The production and degradation rates for mRNA and protein may be modulated by physical or chemical factors; for example, a protein may locally interact with DNA so as to increase (promote) or decrease (inhibit or repress) the production rate for mRNA corresponding to a particular gene. In this case, the DNA-interacting protein is called a transcription factor because it directly influences whether mRNA is transcribed. The genes expressed by cells during embryonic development will determine the ultimate “identity” of the cell (e.g., a nerve or skin cell) in the adult organism.

Here, we consider a simple model of an autoregulatory process in which each cell transcribes mRNA $m$ that is translated into protein $p$ that in turn modifies the production rate of $m$. The signaling molecule $v$, generated in exact proportion to $p$, also regulates $p$ production in the self and neighbors by modulating the production rate of $m$. The system dynamics are

[TABLE]

where $\gamma_{m}$, $\gamma_{p}$ are the degradation (decay) rates of mRNA and protein respectively, and $\alpha_{m}$, $\alpha_{p}$ are the corresponding transcription or translation rates. The function $f(v_{i},u_{i},p_{i})$ captures the influence of the coupling, input, and protein on the production rate of the mRNA and therefore of the protein.

When linearized at steady state, the system becomes

[TABLE]

where F_{v}:=\frac{\partial f}{\partial v_{i}}\big{|}_{(\bar{m}^{*},\bar{p}^{*},\bar{v}^{*},\bar{u}^{*})}, F_{u}:=\frac{\partial f}{\partial u_{i}}\big{|}_{(\bar{m}^{*},\bar{p}^{*},\bar{v}^{*},\bar{u}^{*})}, and {F_{p}:=\frac{\partial f}{\partial p_{i}}\big{|}_{(\bar{m}^{*},\bar{p}^{*},\bar{v}^{*},\bar{u}^{*})}}.

Define $\alpha:=\alpha_{m}\alpha_{p}$ and $\gamma:=\gamma_{m}\gamma_{p}$. Note that the steady-state protein concentration is a linear multiple of the steady-state mRNA concentration, such that mathematically a molecule produced at rate $\alpha f(v_{i},u_{i},p_{i})$ and decayed at rate $\gamma$ would have the same steady-state concentration as the protein in (3). Indeed, it is not uncommon for the dynamics of transcription and translation to be lumped together (usually by neglecting mRNA dynamics) in mathematical models such as those presented later in this paper.

The steady-state solution to the perturbed system yields filter coefficients

[TABLE]

for $k=0,1,...,N-1$. For the homogeneous steady state to be stable—and therefore for the filtering approach to be applicable—we require

[TABLE]

We henceforth assume this condition is satisfied.

Recall that the spatial modes are given by $t_{k}$, the columns of the matrix $T$ that diagonalizes $M$. The perturbing input can be written as

[TABLE]

The coefficients $\hat{u}_{k}$ (the entries of $\hat{u}$) are the weights assigned to each of the spatial modes $t_{k}$. The steady-state perturbed readout is given by

[TABLE]

such that the readout in the $i$th cell is given by $\underline{\bar{y}}+\tilde{y}^{*}_{i}$.

The $\mathcal{H}_{2}$ norm for the $k$th spatial mode is analytically calculated to be $\frac{F_{u}}{2\alpha}\left[S\right]_{kk}$. This relationship indicates that the modes in the system respond identically to within a scaling factor to both persistent spatial disturbances and temporally varying white noise inputs.

Figure 3 exemplifies how the choice of interaction type and interconnectivity affects the filtering behavior of the system with no autoregulation. In particular, activation of neighbors tends to cause the system to amplify low spatial frequencies, while inhibition of neighbors introduces amplification at high spatial frequencies.

To investigate the effect of autoregulation, suppose we fix all parameters except $F_{p}$. As $F_{p}\rightarrow-\infty$ all filter coefficients approach 0. This attenuating behavior occurs because allowing a protein to effectively shut down its own production prevents the system from responding to signal.

For $F_{p}>0$ (autoactivation), increasing $F_{p}$ disproportionately increases the coefficients at spatial modes with low eigenvalues. For $T$ corresponding to the DFT or DCT-2, the lower eigenvalues are associated with higher-frequency spatial modes. In the case of lateral inhibition ($F_{v}<0$), the filter coefficients already amplify high-frequency spatial modes relative to intermediate ones (Figure 3), such that adding autoactivation enhances the filter’s intrinsic highpass characteristics. Indeed, mechanisms involving lateral inhibition and autoactivation have been conjectured to increase the sharpness of boundary formation in systems of patterned cells responding to exponential input [3][8].

3 1D APPLICATION: NOTCH-DELTA

The Notch-Delta patterning mechanism is a lateral inhibition system that is responsible for diverse developmental phenomena including neural and epidermal fate determination in the fruit fly Drosophila melanogaster. Cells produce both Notch and Delta, which are proteins found in the cell membrane. Delta on the surface of one cell binds Notch on the surface of neighboring cells to inhibit those neighbors’ Delta production, thereby relieving inhibition on the cell’s own Delta production by decreasing the potential for the neighbors to bind its Notch. With the appropriate interaction strengths, such mutual inhibition between neighbors will ultimately generate a checkerboard pattern in which cells expressing high Delta are adjacent to cells expressing low Delta. This has significant consequences for organismal development: Notch that is bound by Delta on an adjacent cell will cleave in two—preventing it from further signaling—and the portion left inside the cell will signal the cell to express target genes that influence the choice of cell identity. A cell whose neighbors express more Delta is more likely to have bound Notch and therefore more likely to adopt a particular fate [20], [21].

Patterning in a Notch-Delta system may arise spontaneously [22], [7] or through modification of a prepattern. In the case of Drosophila wing development, the gene veinless is expressed in an exponential gradient decreasing in either direction from what will become the center of a vein. The level of veinless expression in a cell determines the Delta production rate at that cell. Notch activity occurs in two peaks, one on either side of the center, where further vein development is restricted to occur. One model of the Notch-Delta mechanism suggests that so-called mutual inactivation, when Notch and Delta on the same cell inhibit each other’s activity, enables sharper and more robust patterning than is achieved with lateral inhibition alone [8], [23].

The authors of [8] considered a line of cells with periodic boundary conditions, corresponding to the interconnection matrix

[TABLE]

The diagonal entries are zero to reflect the fact that Notch ($N$) and Delta ($D$) interact via contact with neighbors rather than diffusion, while the factor of $\frac{1}{2}$ ensures that $v_{N}$ is the average Notch from neighbors and $v_{D}$ is the average Delta from neighbors. Because $M$ is circulant, the spatial modes correspond to the eigenvectors of the DFT matrix.

We discretize the input gradient of Delta production rate $\beta_{D}(\cdot)$ into $\beta_{D_{i}}$, $i=0,1,...,N-1$ and let $\underline{\bar{\beta}}_{D}$ be the mean of the $\beta_{D_{i}}$. We then define $u_{i}:=\beta_{D_{i}}-\underline{\bar{\beta}}_{D}$, $x_{i}^{T}:=[N_{i},D_{i},R_{i}]$, and $v_{i}:=[v_{N_{i}},v_{D_{i}}]$ where readout $R$ is a reporter for Notch activity (i.e., is expressed from a target gene for Notch activity).

The authors of [23] propose four models of the Notch-Delta patterning mechanism that involve mutual inactivation, lateral inhibition, or both. As an example we present the linearization for the mutual inactivation (MI) model; equations and linearizations for the lateral inhibition with mutual inactivation (LIMI) and simplest lateral inhibition by mutual inactivation (SLIMI) models may be found in Supplementary Section 9.

The system equations for the MI model are

[TABLE]

where $\gamma$, $\gamma_{R}$ are decay rates, $k_{t}^{-1}$ is the rate at which Delta and Notch bind each other on neighboring cells, $k_{c}^{-1}$ is the strength of mutual inactivation, and $k_{RS}$, $n$ are parameters determining how strongly bound Notch promotes reporter expression. Note that mRNA is not explicitly incorporated into the model, such that the dynamics are effectively lumped into the production and degradation terms for the proteins.

Linearization about the steady state with all $u_{i}=0$ yields

[TABLE]

where we have defined

[TABLE]

and $C$ is chosen depending on the readout. It can be shown that the linearized dynamical system is stable for all nonnegative and thus biologically attainable parameter values (Supplementary Section 9.1), strengthening the argument that patterning may not require instability.

To examine the effects of identical (correlated) vs. separate (uncorrelated) white noise inputs to both Delta and Notch, we first modify 5 so that a single input appears in the equations for both $\dot{N_{i}}$ and $\dot{D_{i}}$ in the correlated case and two independent inputs appear in each of these equations for the uncorrelated case. Accordingly, we then calculate the $\mathcal{H}_{2}$ norm for

[TABLE]

3.1 Comparison of Models

The MI, LIMI, and SLIMI models from reference [23] produce substantially similar readouts (Figure 6), filter characteristics, and $\mathcal{H}_{2}$ norms (Figures 19 and 20). Mutual inactivation (lower $k_{c}$) decreases the magnitude of the coefficients and therefore the final Notch and Delta concentrations, but exaggerates the intrinsic highpass characteristics of the filter, producing the sharper peaks in Notch activity predicted by [8]. Analysis of the $\mathcal{H}_{2}$ norm reveals that regardless of readout, noise that is completely uncorrelated between Delta and Notch production rates is favored by the same frequencies as the system filter, while noise that is completely correlated between the production rates is almost uniformly rejected relative to uncorrelated noise (Figures 4 and 5). Together, these observations suggest that time-varying stochastic inputs—unless they are of exceptionally large magnitude—do little to combat the intrinsic behavior of the filter, contributing to the robustness of the developmental program.

4 APPLICATION: DIGIT FORMATION

Digits in developing vertebrate embryos originate from a flat paddle-shaped layer of cells that form the limb bud. A crucial step in digit patterning involves specifying which cells in the paddle will become digits and which will die to create the space between digits [24], [25]. This periodic pattern of digit with interdigit has been proposed to originate with spatially periodic expression of the gene sox9, which produces a protein that regulates transcription of the genes wnt and bmp. In turn, these genes code transcription factors Wnt and Bmp that regulate Sox9 production [26].

Cell cultures from developing embryos grown on plates show Turing-like patterns where Sox9 is out of phase with Wnt and Bmp. Turing patterns typically arise in chemical reaction systems with at least two types of diffusible molecules produced at every point in space, where the activation/inhibition relationship between the types is such that the homogeneous solution to the resulting dynamical system is unstable owing to the difference in diffusion rates between the two molecules. Such a reaction-diffusion model has been proposed to generate the observed Sox9/Wnt/Bmp pattern from stochastic initial conditions within a particular parameter range [26]. Our discretization of the model suggests that such a pattern might be observed even if the parameters do not satisfy the conditions for diffusion-driven instability.

Consider the Sox9/Bmp/Wnt network with diffusion distance $l$ between cells. Let $s$, $b$, and $z$ represent the concentrations of Sox9, Bmp, and Wnt respectively, such that $x_{i}=[s_{i},~{}b_{i},~{}z_{i}]$ and $v=[v_{b},~{}v_{z}]$. Let the input be random cell-to-cell variation in background protein production rate, i.e., the production rate of protein in the absence of promotion or inhibition, as from cell-to-cell variability in transcription or translation rates (see Supplementary Figure 13). The dynamics within cell $i$ and the coupling are given by

[TABLE]

where $\alpha$ are background production rates, $k$ are interaction rates, and $d$ are diffusivities.

Linearization about the homogeneous steady state yields

[TABLE]

4.1 Spatial Modes in 2D

For this example we will consider a two-dimensional, rectangular $N_{R}\times N_{C}$ array of $N:=N_{C}N_{R}$ cells indexed from [math] to $N_{C}N_{R}-1$ starting in the upper lefthand corner from top to bottom and then left to right, i.e.,

[TABLE]

It is known (e.g., [11]) that if any isolated row has interconnection matrix $M_{R}\in\mathbb{R}^{N_{C}\times N_{C}}$ and any isolated column has interconnection matrix $M_{C}\in\mathbb{R}^{N_{R}\times N_{R}}$, then the full matrix $M$ for the interconnectivity of the entire array is

[TABLE]

If $T_{R}$ and $T_{C}$ diagonalize $M_{R}$ and $M_{C}$ respectively then $M$ is diagonalized by

[TABLE]

giving $N_{C}N_{R}$ eigenvalues

[TABLE]

where $m=0,1,...N_{R}-1$, $n=0,1,...,N_{C}-1$. The $(m,n)$th spatial mode is given by $T_{C}^{m}T_{R}^{nT}$.

We can explicitly relate the spatial modes for a 2D array of cells to constituent modes in the horizontal and vertical directions by recasting the vector $\hat{y}$ in matrix form. Let $U$ and $Y$ be matrices arranged as in (7) where $u_{i}$ is the input to compartment $i$. Vector form is recovered through the vectorization operation ${vec(U)=u}$. The readout matrix $Y$ is defined similarly. If the matrices $\tilde{U}$ and $\tilde{Y}$ designate perturbations from steady state in the original basis and $\hat{Y}$, $\hat{U}$ designate perturbations in the basis for the spatial modes, then

[TABLE]

where $\odot$ is the Hadamard product (element-by-element multiplication) and $\Lambda_{S}$ has $m$th, $n$th entry

[TABLE]

From this it can be seen that the full system alters the input along the $i$th vertical spatial mode and the $j$th horizontal spatial mode defined by the vertical and horizontal connectivities. For the remainder of this example, we will assume Neumann boundary conditions such that the spatial modes for the rows and columns of $M$ correspond to the DCT-2 (Figure 7).

4.2 Analysis

We pick $C$ to monitor Sox9, Bmp, or Wnt concentration and choose $s_{0}$ such that the Turing instability conditions are not satisfied, i.e., the eigenvalues of $\left(I_{N}\otimes A\right)+M\otimes\left(B_{v}G\right)$ are all negative. Nevertheless, the readout still replicates the spatially periodic patterns predicted by [26] for a range of intercellular distances (Figure 10) owing to the bandpass behavior of the filter (Figure 8). [Sox9] is out of phase with both [Bmp] and [Wnt], as indicated by the opposing signs of the coefficients in the passband.

The $\mathcal{H}_{2}$ norm measurements for the readouts qualitatively emphasize the same frequency bands as their respective filters $\left[S\right]_{kk}$ (Figure 9). For correlated or uncorrelated noise sources, readouts [Bmp] and [Wnt] experience much greater magnification than does [Sox9], suggesting that the opposing effects of Bmp and Wnt on sox9 expression may mostly cancel each other out at the level of Sox9 concentration. Relative noise amplification in the same modes favored by the $\left[S\right]_{kk}$ may ensure that stochastic influences do not counteract filter behavior, at the same time that attenuation and evenness in the response to other modes might reduce the relative influence of temporally varying inputs on the readout. The latter especially may be useful to maintain consistent behavior in a process such as digit formation that takes place over a long timespan.

While our simulations do not refute the hypothesis that a diffusion-driven instability constitutes the biological basis for digit formation, the fact that we can produce a similar pattern with an externally perturbed stable system suggests that not all apparent Turing patterns need arise from an instability. This observation could significantly ease the search for molecules and proteins that contribute to “spontaneous” stripe and spot patterning, as the parameter restrictions required for true Turing instabilities may not be biologically plausible.

5 APPLICATION: DIGIT FORMATION WITH A MORPHOGEN GRADIENT

Expanding on the work of [26], reference [27] demonstrated that changes to the parameters in the proposed Turing network for digit formation in mice can produce sox9 expression patterns matching those found in embryonic catshark fins, suggesting that the mechanism has been evolutionarily conserved. The authors augmented the model with an exponential gradient of fibroblast growth factor (Fgf), a morphogen originating at the fin edge that has been experimentally demonstrated to facilitate normal digit arrangement in mice. In their model, Fgf represses Sox9 repression of $bmp$ expression ($k_{4}$) and promotes Sox9 repression of $wnt$ expression ($k_{7}$). In simulation, the authors observed that increasing the ratio of Wnt production to Bmp production or decreasing Bmp promotion of $sox9$ expression caused the Turing pattern to transition from stripes to spots.

We implemented the model from [27] using the following evolution equations:

[TABLE]

where $\alpha_{F}$ is the Fgf production rate, $u_{1}$ represents the source of Fgf, and $u_{2}$ is random constant-in-time spatial variation in background production rate. Unlike [27], we did not normalize $F$ to $[0,1]$, but we chose $u_{1},\underline{u}_{1}$ such that $0\leq\tilde{F_{i}}^{*}+\underline{\bar{F}}^{*}\leq 1$ and $0\leq u_{1_{i}}-\underline{\bar{u}}_{1}$.

As compared to (6), the equations in (8) are rendered as perturbations to prior steady-state protein concentrations, therefore “negative” steady-state values should be interpreted as reductions in concentration relative to preexisting levels.

To handle both background production rate and localized Fgf production we use the generalization to $L$ inputs

[TABLE]

where $u_{k}\in\mathbb{R}^{N}$ is the $k$th input vector and

[TABLE]

is the linearization matrix for one subsystem with respect to the $k$th input when all inputs are held constant in time and space. To avoid ambiguity, the “filter” interpretation is defined with respect to one input, i.e., as one term in the summation (9).

The matrices for the linearization are

[TABLE]

where the steady-state concentration of Fgf is $\underline{\bar{F}}^{*}=\frac{\alpha_{F}+\underline{\bar{u}}_{1}}{\mu_{F}}$ independent of the other variables. We stabilized the homogeneous steady-state solution by setting $\alpha_{F}=0$ and using a small value of $\underline{\bar{u}}_{1}$ with the remaining parameters taken from Figures 4 and 5 in [27]. This choice of $\alpha_{F}$ completely localizes the source of Fgf to the input $u_{2}$.

5.1 Spatial Modes on a Hexagonal Lattice

The hexagonal lattice is the tightest 2D packing arrangement for cells of fixed area and is found in a number of natural systems such as the wing epithelial cells in Drosophila [28]. For this example we will derive the lattice from a rectangular array where the columns are offset by 30∘ from vertical and assume periodic boundary conditions as well as identical spacing between all neighbors. With the cells numbered as shown in Figure 11, the $(m,n)$th spatial mode corresponds to the $m$th mode horizontally and the $n$th mode on a line at a 60∘ angle from each row. Cells in the hexagonal lattice interact with each of their six nearest neighbors such that there are “diagonal interconnections” between rows of cells. We account for the diagonal connections as follows: Define $C_{D}:=min\left(N_{R},N_{C}\right)$ and $R_{D}:=max\left(N_{R},N_{C}\right)$, and let $P\in\mathbb{R}^{C_{D}\times C_{D}}$ be the permutation matrix with lower diagonal ones and the last entry of the first column also one. Define $P_{F}:=diag\left(P^{0},P^{1},...,P^{R_{D}}\right)$. Then the interconnection matrix for a hexagonal lattice with periodic boundary conditions is given by

[TABLE]

If we let $M_{R}=M_{C}=M_{F}=M_{0}\in\mathbb{R}^{N_{0}\times N_{0}}$ be the circulant diffusion matrix, then $M$ has eigenvalues

[TABLE]

From this we see that $\Lambda=\Lambda^{T}$ and $\left[\Lambda\right]_{mn}=\left[\Lambda\right]_{(N_{0}-m)(N_{0}-n)}$.

Further discussion of diagonal interconnectivites and planar lattices more generally is available in Supplementary Sections 8.3 and 8.4.

5.2 Analysis

As in the original Sox9-Bmp-Wnt Turing model in [26], the filter coefficients in the Fgf-augmented model form a bandpass at mid-range frequencies, resulting in the roughly periodic output patterning that alternates between Sox9 and Wnt (Figure 12). Increasing $\alpha_{wnt}$ decreases the magnitude of the bandpass, while decreasing $k_{2}$ concentrates amplification at a small range of frequencies inside the bandpass. Either of these parameter changes tends to shrink contiguous regions of high [Sox9], consistent with the transition from stripes to spots observed in [27]. Parameters yielding more spotlike patterns also tend to suppress the influence of both correlated and uncorrelated noise for readout [Sox9], though the effect on readout [Wnt] is negligible (Supplementary Figures 27, 28, 29, and 30). The distal edge where the source of Fgf is localized exhibits relatively higher Wnt than Sox9 expression, as observed in vivo [27]; in our normalized images, the effect is most visible at higher values of $\alpha_{wnt}$.

If we assume cells are approximately $12$ to $15$ $\mu$m in diameter [29], then for $\alpha_{wnt}=1.2,k_{2}=1$ the filter and $\mathcal{H}_{2}$ norm analysis indicate that wavelengths of about 84 to 105 $\mu$m will be most strongly amplified in the result. The prediction is in decent agreement with the experimental images in Figure 2 of [27], which exhibit periodicity on the order of 80 to 100 $\mu$m. Some of the error may be accounted for by the difference in domain shape between filter simulations and actual limb paddles (rhomboidal vs. elliptical) as well as the presence of growth in the living animal. Nevertheless, this observation suggests that the framework correctly identifies the range of spatial modes that will be most influential in forming the “actual” biological pattern.

6 CONCLUSIONS

In this paper we have presented a framework to analyze how networks of interacting cells modify spatially varying inputs, either from environmental factors or intrinsic parameter variation, to produce patterned outputs. Three biologically relevant examples indicate that qualitatively similar patterns may arise from different physical implementations (Section 3), from both stable and unstable fixed points (Section 4), as well as from variable filter behaviors when certain postprocessing steps are applied (Section 5). Furthermore, these biological models appear robust to correlated and uncorrelated space-and-time-varying white noise inputs, a critical feature for maintaining consistency during embryonic development.

We have demonstrated in a theoretical context how a filtering approach can offer insight into system behavior at an intermediate level between the exact physical implementation and the measured result. In an experimental context, evaluating systems at the filter level may clarify when alterations to the input are capable of distinguishing between alternative explanations for an observed behavior. For example, systems with near-identical filter coefficients are predicted to respond equivalently to inputs of all kinds (e.g., the three Notch-Delta models in Section 3), suggesting that pure input-output probing is unlikely to illuminate the underlying mechanism. Conversely, model systems with disparate filter coefficients may not vary much in their response to certain inputs but differ drastically in reponse to others, such that experiments in which inputs to the real system can be finely controlled may suffice to differentiate more accurate models from less accurate ones. Of interest in both cases is the extent to which a particular system may impose structure upon an output pattern as compared to how much structure must be present in the prepattern.

A critical assumption in our development of the filtering framework is that linearization about a homogeneous steady state is sufficient to capture relevant system behavior. Future work should focus on incorporating nonlinear dynamics as well as investigating the influence of external inputs on spatially distributed, networked systems in the vicinity of unstable or nonhomogeneous steady states. Additional areas for further research include patterning in time-varying or perturbed networks and system response to non-white noise inputs. Lastly, although we have focused our applications on models in developmental biology, the generality of our framework suggests possible applications to synthetically engineered biological circuits as well.

Overall, we believe a spatial frequency-based interpretation simplifies the process of predicting how intermolecular and intercellular interactions affect patterning mechanisms in living organisms. It is our hope that the viewpoint developed here will help us to elucidate—and elaborate upon—nature’s designs.

ACKNOWLEDGMENTS

The authors would like to thank Andy Packard for providing incisive feedback on an earlier version of this work, and an anonymous reviewer for suggestions to strengthen the presentation of the manuscript.

7 MORE ON FILTER COEFFICIENTS

7.1 Derivation of Filter Coefficients

Here we provide a derivation of the filter coefficients from Proposition 1.

Let $\underline{u}:=\underline{\bar{u}}\mathds{1}_{N}$ be a spatially homogeneous input and assume $\exists~{}\underline{\bar{x}}^{*}\in\mathbb{R}^{n}$ such that $f\left(\underline{\bar{x}}^{*},\mu g\left(\underline{\bar{x}}^{*}\right),\underline{\bar{u}}\right)=0$. Then $\underline{x}^{*}:=\mathds{1}_{N}\otimes\underline{\bar{x}}^{*}$ is a homogeneous steady state. The remaining steady-state quantities are similarly designated $\underline{y}^{*}=\underline{\bar{y}}^{*}\mathds{1}_{N}$, $\underline{\bar{w}}^{*}=g(\underline{\bar{x}}^{*})$, and $\underline{v}^{*}=M_{q}\underline{w}^{*}=\mathds{1}_{N}\otimes\underline{\bar{v}}^{*}$ (where $M_{q}:=M\otimes I_{q}$). Let $\tilde{x}_{i}(t),\tilde{u}_{i},\tilde{y}_{i}(t),\tilde{w}_{i}(t),\tilde{v}_{i}(t)$ denote perturbations about that steady state. The full system linearized about $(\underline{x}^{*},\underline{u})$ yields perturbed dynamics

[TABLE]

where A:=\frac{\partial f}{\partial x_{i}}\big{|}_{(\underline{\bar{x}}^{*},\underline{\bar{v}}^{*},\underline{\bar{u}})}, B_{v}:=\frac{\partial f}{\partial v_{i}}\big{|}_{(\underline{\bar{x}}^{*},\underline{\bar{v}}^{*},\underline{\bar{u}})}, B_{u}:=\frac{\partial f}{\partial u_{i}}\big{|}_{(\underline{\bar{x}}^{*},\underline{\bar{v}}^{*},\underline{\bar{u}})}, C:=\frac{dh}{dy_{i}}\big{|}_{\underline{\bar{x}}^{*}}, and G:=\frac{dg}{dx_{i}}\big{|}_{\underline{\bar{x}}^{*}} are the linearization matrices.

Assume that the interconnection matrix $M\in\mathbb{R}^{N\times N}$ is diagonalizable and let $M=T\Lambda T^{-1}$ be the diagonalization (so that $M_{q}$ is diagonalized by $T\otimes I_{q}$). Define $\hat{x}(t):=\left(T^{-1}\otimes I_{n}\right)\tilde{x}(t)$, $\hat{u}:=T^{-1}\tilde{u}$, and $\hat{y}(t):=T^{-1}\tilde{y}(t)$. Recasting (10) in the coordinate system $T$, we obtain the dynamical system

[TABLE]

In contrast to the conditions for spontaneous pattern formation, we will not require this system to be unstable; large amplification of spatial modes is possible even when the system is stable. The steady-state perturbed readout in basis $T$ is

[TABLE]

where $S$ is a diagonal matrix with entries

[TABLE]

for $k=0,1,...,N-1$ and $\left[S\right]_{kk}$ collectively form the “filter coefficients” for the corresponding $N$ spatial modes. The matrix $S$ is thus analogous to a digital filter that processes the input $\tilde{u}$ into readout $\tilde{y}^{*}$ with respect to the eigenvectors, or spatial modes, of $M$ as contained in $T$.

7.2 Multiple Orthogonal Signals Sharing Same Spatial Modes

Let $M_{0},M_{1},...,M_{q-1}\in\mathbb{R}^{N\times N}$ be the interconnection matrices for $q$ orthogonal signals and assume they all commute (share the same basis). Let $\Delta_{i}\in\mathbb{R}^{q\times q}$ signify the matrix with $(i,i)$th entry one and all other entries zero, such that the full interconnectivity is

[TABLE]

If $\Lambda_{i}=T^{-1}M_{i}T$ for $i=0,1,...,q-1$, then the basis $T\otimes I_{q}$ diagonalizes $M$ such that the filter coefficients are given by

[TABLE]

7.3 Derivation of Equivariance

Here we derive the equivariance property of the input-output map comprising the filter coefficients.

Let $\tilde{y}^{*}=S\tilde{u}$ and $\tilde{y}_{\Pi}^{*}=S\Pi\tilde{u}$. An equivalent statement to “$S$ is equivariant” is then $\tilde{y}_{\Pi}^{*}=\Pi S\tilde{u}\implies\tilde{y}_{\Pi}^{*}=\Pi\tilde{y}^{*}$. To see this, take

[TABLE]

Since

[TABLE]

then

[TABLE]

such that (13) becomes

[TABLE]

If $\Lambda=T^{-1}MT$, then since $M=\Pi M\Pi^{-1}\implies\Pi^{-1}M\Pi$, $\Lambda=T^{-1}\Pi^{-1}M\Pi T$ (i.e., $T$ diagonalizes the permuted version of $M$ with the same resultant eigenvalues). Note that the immutability of $M$ under permutation $\Pi$ confers immutability of $\Lambda$ under permutation $T^{-1}\Pi T$, which is just the permutation in the basis of $M$; i.e.,

[TABLE]

or equivalently, $\Lambda$ and $\left(T^{-1}\Pi T\right)$ commute.

8 TUTORIAL ON SPATIAL MODES

Here we present a brief introduction to two sets of 1D spatial modes corresponding to common signal processing transforms. These spatial modes have direct interpretations as spatial frequencies. We then provide two useful observations for calculating the spatial modes and eigenvalues for 2D arrays with diagonal interconnections, as well as an interpretation of spatial frequencies for cells arranged in arbitrary planar lattices.

8.1 Discrete Fourier Transform (DFT)

If the $N$ cells form a ring indexed clockwise or counterclockwise, then $M$ is circulant. The eigenvectors of a circulant matrix form the discrete Fourier basis such that the spatial modes of $T$ correspond exactly to the frequencies of sinusoids.

We can choose $T$ to be the discrete Fourier transform matrix (DFT) where the $j$th entry of the $k$th eigenvector, $j,k=0,1,...,N-1$, is given by

[TABLE]

with $i:=\sqrt{-1}$. $T$ is conjugate symmetric. If we let $m_{0},~{}m_{1},~{}...,m_{N-1}$ denote the entries in the first row of $M$, then the eigenvalues of $M$ are given by

[TABLE]

which corresponds to the coefficients of the discrete Fourier transform (DFT) of the first row of $M$.

If $M$ is symmetric in addition to circulant, then we can alternatively select the eigenvectors such that all entries are real.

Observation 1**.**

Let $M\in\mathbb{R}^{N\times N}$ be a symmetric circulant matrix where $m_{0}\in\mathbb{R}^{N}$ is the first row and define $m$ as the periodization of $m_{0}$. Let the matrix $T$ have entries

[TABLE]

Then $T$ is a basis for $M$ with eigenvalues

[TABLE]

$k=0,1,...,N-1$.

Proof.

Let $W$ be the unitary DFT matrix, i.e., the $j$th entry of the $k$th column is

[TABLE]

We can express $T$ as

[TABLE]

Then

[TABLE]

Since $M$ is real and even (symmetric), the DFT is also real. The symmetry of $M$ together with the symmetry of $W$ and the fact that diagonal matrices are symmetric also imply that $W^{H}MW=(W^{H}MW)^{T}=W^{T}M^{T}W^{*T}=WMW^{H}$, so we can somewhat simplify (14) to

[TABLE]

For $k=1,2,...,\frac{N}{2}$ (N odd) or $k=1,2,...,\frac{N-1}{2}$ (N even), the $(N-k)$th row or column of $W$ is equal to the $k$th row or column of $W^{H}$. Therefore $\left[WW\right]_{k,N-k}=1$. Because of the orthogonality of complex exponentials, the remaining entries are [math]. By the same logic we deduce an identical structure for $W^{H}W^{H}$ such that $WW=W^{H}W^{H}$. Now (15) becomes

[TABLE]

which is just $M$ diagonalized by the complex exponential DFT matrices, as desired. From this we derive that the eigenvalues are the same as the DFT coefficients of $h_{0}$, which owing to symmetry may be calculated as

[TABLE]

∎

Owing to the periodicity of cosine, the eigenvalues of symmetric circulant $M$ (and hence the corresponding $\left[S\right]_{kk}$) are symmetric about the highest-frequency eigenvector associated with $k=\frac{N}{2}$.

A situation of particular interest occurs when $v_{i}$ is a diffusible molecule and the connection strength is equal between cells. Then $M$ is a scaled version of the circulant finite differences (Laplacian) matrix

[TABLE]

with eigenvalues

[TABLE]

This form of $M$ corresponds to the second differences matrix for a system with periodic boundary conditions [19].

8.2 Second Discrete Cosine Transform (DCT-2)

If the $N$ cells are organized in a line, then the two cells on the end each communicate with only one neighbor. If the mode of communication is a diffusible molecule and the cells are indexed from one end of the line to the other, then the connectivity takes the form of a second differences matrix with Neumann boundary conditions centered at the midpoint:

[TABLE]

The spatial modes of $T$ form the basis for the second discrete cosine transform (DCT-2). The $j$th entry of the $k$th eigenvector, $j,k=0,1,...,N-1$, is given by

[TABLE]

(for $k=0$, divide by additional factor of $\sqrt{2}$) with corresponding eigenvalue

[TABLE]

The highest frequency is $k=N-1$ [19]. Unlike the case of circulant $M$, there are no guarantees of symmetry in the filter $S$ and $T$ itself is not conjugate symmetric.

8.3 Eigenvectors and Eigenvalues for Arrays with Diagonal Interconnections

Observation 2**.**

Consider an $N_{R}\times N_{C}$ array of cells. Let $M_{F}$ describe the interconnectivity of $C_{D}:=min(N_{R},N_{C})$ elements in the forward diagonal direction and $M_{B}$ the interconnectivity of $C_{D}$ elements in the backward diagonal direction. Let $P\in\mathbb{R}^{C_{D}\times C_{D}}$ be the permutation matrix with lower diagonal ones and the last entry of the first column also one, such that $P_{F}:=P_{B}^{T}:=diag\left(P^{0},P^{1},P^{2},...,P^{R_{D}}\right)$ where $R_{D}:=max(N_{R},N_{C})$. In total, the interconnectivity of an array with horizontal, vertical, and diagonal components is described by

[TABLE]

where $M_{R}$, $M_{C}$, $M_{F}$, and $M_{B}$ are circulant. Furthermore, if $N_{R}=N_{C}=N_{0}$ and the real or complex DFT basis $T_{0}$ diagonalizes each of $M_{R}$, $M_{C}$, $M_{F}$, and $M_{B}$ individually, then $\left(T_{0}\otimes T_{0}\right)$ diagonalizes $M$.

Proof.

Let $T_{0}$ diagonalize $M_{0}$ such that $\left(T_{0}\otimes T_{0}\right)$ diagonalizes $\left(M_{R}\otimes I_{N_{R}}\right)+\left(I_{C}\otimes M_{C}\right)$. Since $M_{0}$ is circulant,

[TABLE]

and therefore

[TABLE]

If we expand $P_{F}$ and use the fact that $P^{-1}=P^{N_{R}-1}=P^{T}$, then the summation simplifies to

[TABLE]

which is diagonalized by the complex exponential DFT vectors $T_{0}^{-1}$ and $T_{0}$ as follows:

[TABLE]

which is a sum of diagonal matrices because permutation matrices are circulant and therefore diagonalized by the DFT matrices $T_{0}$, and the Kronecker product of two diagonal matrices is diagonal. Since $P_{B}=P_{F}^{T}=P_{F}^{-1}$, the same derivation for diagonal connectivity in the backward direction gives

[TABLE]

which is also diagonalized by DFT matrices as

[TABLE]

where each summand is diagonal. This implies that the complex exponential form of the DFT is the basis for an array that is diagonally connected in either or both directions.

If $M_{0}$ is symmetric, then in addition to the complex exponential basis $T_{0}$ we might also choose the basis $T^{\prime}_{0}$ with the $j$th entry of $k$th eigenvector given by

[TABLE]

in which case each term in the summation (19) is no longer diagonal because $T^{\prime}_{0}$ does not diagonalize the nonsymmetric permutation matrices. Hence for the real-valued basis $T^{\prime}_{0}$ we require that the array is diagonally connected in both forward and backward directions such that the full connectivity matrix is given by the sum of (20) and (18):

[TABLE]

The matrix $\left(P^{-q}+P^{q}\right)$ is circulant and symmetric and hence diagonalized by $T^{\prime}_{0}$. To complete the argument we appeal to the symmetry of $m_{q}$. Specifically, if $N_{R}$ is odd,

[TABLE]

where each individual term is diagonalized by $T^{\prime}_{0}$, and hence the whole summation is diagonal. If $N_{R}$ is even, the summation (22) breaks into

[TABLE]

When $N_{R}$ is even, $P^{R/2}$ alone is circulant symmetric, hence the additional term is also diagonalized by $T^{\prime}_{0}$. Therefore when $M_{0}$ is circulant symmetric, the matrix

[TABLE]

is diagonalized by $T^{\prime}:=T^{\prime}_{0}\otimes T^{\prime}_{0}$. Note that this implies

[TABLE]

is also diagonalized by $T^{\prime}$. ∎

Remark 1**.**

For a system with diagonal connections only ($M_{R}=M_{C}=0$), then if $N_{R}=N_{C}$ odd, the diagonal transformations are identical to a 2D DFT rotated 45∘. For $N_{R}=N_{C}$ even, the array becomes divided into two separate classes that are transformed separately; i.e., the underlying network graph is no longer connected. This is because for $N_{R}=N_{C}$ odd, the array has a compartment at the center, while for $N_{R}=N_{C}$ even, the center would (in physical space) represent a crossing of intersections.

Observation 3**.**

Let $N_{C}=N_{R}=N_{0}$, $M_{0}$ circulant with complex exponential basis vectors $T_{0}$ and consider the full forward diagonal interconnection matrix $M=P_{F}^{T}\left(M_{0}\otimes I_{N_{0}}\right)P_{F}$. The $(m,n)$th eigenvalue of $M$ is

[TABLE]

where $m_{k}$ is the $k$th entry of the first row or column of $M_{0}$.

Proof.

For $N_{C}=N_{R}=N_{0}$ and $M_{F}$ as defined above, the $k$th entry of the $N_{0}$-point DFT of the first row of $P^{q}$ is $e^{\frac{2\pi ikq}{N_{0}}}$, $k=0,1,...,N_{0}$ and the diagonalization $T_{0}^{-1}P^{q}T_{0}$ is the matrix with the DFT entries on the diagonal. This implies that we can write (19) as

[TABLE]

where we have defined

[TABLE]

For $N_{0}$ odd, we use the symmetry $m_{q}=m_{-q}$ to rewrite the summation (23) as

[TABLE]

Conveniently,

[TABLE]

from which we infer

[TABLE]

is the eigenvalue for the $\left(m,n\right)$th spatial mode owing to diagonal connectivity, $N_{0}$ odd. If $N_{0}$ is even, we write (23) as

[TABLE]

and note that

[TABLE]

to write the $(m,n)$th eigenvalue as

[TABLE]

∎

8.4 Spatial Modes on Planar Lattices

It is straightforward to generalize a frequency-based interpretation to cells arranged in periodic planar lattices, which are well described mathematically. Throughout the following discussion we will refer to coordinates in physical space as $\hat{e}$, the unit vector pointing “east,” and $\hat{s}$, the unit vector pointing “south”. This choice of vector orientations mimics the numbering scheme in an array, whereby indices increase horizontally left to right (with $\hat{e}$) and vertically top to bottom (with $\hat{s}$). We will assume a system of cells indexed in an $N_{R}\times N_{C}$ array with periodic boundary conditions such that $M$ is diagonalized by $T_{R}\otimes T_{C}$ where both $T_{C}$ and $T_{R}$ are real or complex DFT bases of appropriate dimension.

Let cells in physical space be arranged in a planar lattice described by vectors $a_{R}$ and $a_{C}$ corresponding respectively to the rows and columns of the indexed array. Without loss of generality we orient $a_{R}$ along $\hat{e}$ (such that $a_{R}\cdot\hat{e}=|a_{R}|$). We define the unit vectors $\hat{a}_{R}:=\frac{a_{R}}{|a_{R}|}=\hat{e}$ and $\hat{a}_{C}:=\frac{a_{C}}{|a_{C}|}$. Letting $\theta$ be the angle between $a_{R}$ and $a_{C}$, we can write $a_{R}=|a_{R}|\hat{e}$ and $a_{C}=|a_{C}|\cos\theta\hat{e}+|a_{C}|\sin\theta\hat{s}$. Note that the eigenfunctions are periodic in $n$ with period $N_{C}|a_{R}|$ along $a_{R}$ and periodic in $m$ with period $N_{R}|a_{C}|$ along $a_{C}$.

Observation 4**.**

For an $N_{R}\otimes N_{C}$ cellular lattice with lattice vectors $a_{R}$, $a_{C}$ and periodic boundary conditions, the $(m,n)$th spatial mode corresponds to a plane wave of frequency

[TABLE]

in physical space, with an “absolute” frequency of

[TABLE]

pointing at an angle

[TABLE]

from the $\hat{x}$ axis.

Remark 2**.**

One may liken the translation from physical space into matrix space to “sampling” in space from an underlying pattern with “spatial sampling frequency” $\frac{1}{|a_{R}|}$ along $a_{R}$ and $\frac{1}{|a_{C}|}$ along $a_{C}$. The translation into spatial modes, or the DFT, recovers normalized frequency components from the discrete samples. Maintaining a constant surface area but increasing the number of cells occupying that surface area (i.e., $N_{R}^{\prime}=c_{R}N_{R}$, $N_{C}^{\prime}=c_{C}N_{C}$, $a_{R}^{\prime}=\frac{1}{C_{R}}a_{R}$, $a_{C}^{\prime}=\frac{1}{c_{C}}a_{C}$) does not change the physical range of space over which the modes are described but does increase the “resolution” or “sampling rate” of the system by a factor of $c_{R}$ along $a_{R}$ and $c_{C}$ along $a_{C}$, enabling the system to modify higher frequencies than before and therefore permitting finer filtering of a continuous-in-space input gradient.

9 NOTCH-DELTA MODELS

9.1 Mutual Inactivation (MI)

For this system we can explicitly calculate the steady-state values $\underline{\bar{N}}^{*},\underline{\bar{D}}^{*},\underline{\bar{v}}_{N}^{*},\underline{\bar{v}}_{D}^{*}$ for $u=0$. First we note that for our choice of $M$ the homogeneous solution satisfies $\underline{\bar{v}}_{N}^{*}=\underline{\bar{N}}^{*}$ and $\underline{\bar{v}}_{D}^{*}=\underline{\bar{D}}^{*}$. After algebra, we find that $\underline{\bar{N}}^{*}$ is the positive root of a quadratic, $\underline{\bar{D}}^{*}$ is found in terms of $\underline{\bar{N}}^{*}$, and $\underline{\bar{R}}^{*}$ is expressed in terms of $\underline{\bar{N}}^{*}$ and $\underline{\bar{D}}^{*}$:

[TABLE]

where $K:=\frac{k_{c}k_{t}}{k_{c}+k_{t}}$.

The filter coefficients $[S]_{kk}$ are given by $-C\left(A+\lambda_{k}(M)B_{v}G\right)^{-1}B_{u}$, $k=0,1,...,N-1$. To find them we can exploit the structure of $C$ and $B_{u}$. For the sake of demonstration we will take the readout to be the reporter protein such that $C=[0~{}0~{}1]$, although the procedure applies equally well to arbitrary choices of $C$.

First we notate

[TABLE]

and apply the matrix inversion lemma to obtain

[TABLE]

Observe that

[TABLE]

Premultiplying (25) by $C$ extracts the bottom row, while postmultiplying by $B_{u}$ extracts the middle entry of that row, which is given by

[TABLE]

Substituting $b_{1}$ and $b_{2}$, we simplify the expression to

[TABLE]

where

[TABLE]

The dynamical system corresponding to these filter coefficients is analytically stable for our chosen $M$ with any biologically relevant parameter values (i.e., when the parameters in 1 are positive, as they must be in a living system). Since $A+\lambda_{k}\left(M_{0}\right)B_{v}G$ is a block triangular matrix, its eigenvalues are the eigenvalues of the diagonal blocks, i.e., $-\gamma_{R}$ along with the eigenvalues of $A_{1}$. Since $-\gamma_{R}$ is always negative, checking for stability amounts to checking the sign of the eigenvalues of $A_{1}$.

Using the fact that $\underline{\bar{v}}_{D}^{*}=\underline{\bar{D}}^{*}$ and $\underline{\bar{v}}_{N}^{*}=\underline{\bar{N}}^{*}$ for our choice of $M$ yields

[TABLE]

with $K$ as and $A_{1}$ as defined earlier. The eigenvalues are given by the zeros of the characteristic polynomial, found by solving for $s$ in

[TABLE]

The first term multiplies out to

[TABLE]

and the second contributes the following terms, independent of $s$:

[TABLE]

To be biologically attainable the parameters and steady-state values must all be positive, such that the quadratic in $s$ has positive coefficients for the first- and second-order terms. If the roots are complex then assuming nonzero decay and nontrivial solutions, the real part is given by

[TABLE]

guaranteeing stability.

If the roots are real, then they will be negative if the zeroth-order term is positive. However, if the zeroth-order term is negative, then one root will be positive and the system will not be stable. Neglecting the expressions in $\gamma$, which by observation must be positive, the contributions to the zeroth-order term are

[TABLE]

By our choice of $M_{0}$, $\lambda_{k}\left(M_{0}\right)\in\left[-1,1\right]$, therefore the bottom expression is minimized to [math] by $\lambda_{k}\left(M_{0}\right)=1$. Since this expression contains all the possible negative contributions to the zeroth-order term, the overall zeroth-order term cannot be negative, and hence the roots of the characteristic polynomial must be negative, implying stability of the system with given $M_{0}$ for all biologically relevant parameter choices.

9.2 Lateral Inhibition with Mutual Inactivation (LIMI)

The system equations are the same as for the mutual inactivation model (5), except that now Delta production is repressed by reporter protein:

[TABLE]

When linearized at steady state, the relevant matrices are

[TABLE]

where $b_{1},b_{2}$ are defined as before and

[TABLE]

9.3 Simplest Lateral Inhibition by Mutual Inactivation (SLIMI)

The system equations are

[TABLE]

When linearized at steady state, the relevant matrices are

[TABLE]

where now

[TABLE]

10 DIGIT FORMATION