From Toggle to Tuning: Controlling Turing Patterns in Gene Circuits

TL;DR

Analyzing network size reveals a key trade-off: small networks are easier to control but less robust, while larger networks gain robustness at the cost of tunability-suggesting a sweet spot for both evolvability and designability.

Abstract

Controlling spatial patterns in synthetic biological systems remains challenging due to poor parameter robustness and limited experimental tunability. We introduce two complementary mechanisms-the pattern switch and the pattern dial-to systematically control Turing pattern formation in gene circuits. The switch toggles pattern onset via a single parameter, while the dial enables transitions between distinct pattern types using weakly nonlinear amplitude equations. Analyzing network size reveals a key trade-off: small networks are easier to control but less robust, while larger networks gain robustness at the cost of tunability-suggesting a sweet spot for both evolvability and designability. Our results offer practical design rules for engineering programmable patterns in living systems.

Figures (5)

• Figure 1: Controlling the dispersion relation for pattern formation in genetic networks.(A-C) Representation of the genetic networks used in this work, corresponding to the 2-node (A), 3-node (B) and 4-node (C) networks, all capable of producing Turing patterns. (D) Pattern switch using the dispersion relation to turn on and off pattern formation. If the dispersion relation is below zero for all $k$, the system is stable and there is no pattern formation, corresponding to the switch being off. Through a parameter change, the dispersion relation can become positive, turning the switch for pattern formation on. Specifically for our initial parameter set, we found that $\mu_A$ can shift the dispersion relation, acting as a switch. (E) Pattern dial between different region of pattern stability. Weakly nonlinear analysis techniques is applied to investigate the stability of different patterns in parameter space, allowing us to find parameters that are capable of causing this transition. In this example, a low value of $k_A$ leads to stripes, while increasing it leads to spots.
• Figure 2: Control of the dispersion relation for pattern formation.(A) Dispersion relations for different values of the parameter $\mu_A$, identified as optimal for reducing the peak for the initial parameter set. When $\mu_A = 0.002$, the dispersion relation remains below zero, indicating that the system is stable with respect to diffusion and does not produce patterns. As $\mu_A$ increases, $\lambda_{max}$ approaches the critical value of 0 at approximately $\mu_A=0.004$. Beyond this point, the system becomes unstable to diffusion and is capable of forming patterns, in agreement with linear stability analysis. A more complete plot showing the behavior of $\lambda_{max}$ as a function of both $\mu_a$ and $k$ is shown in Fig. S1C of the SI. Parameters used: $D_A = 1, D_B = 1000, V_A = 21.544, V_B = 0.464, k_A = 21.544, k_{AB} = 21.544, k_{BA} = 0.464, b_A = 0.010, b_B = 0.010, \mu_A = 0.002, \mu_B = 0.010$. (B) Bar charts showing the number of Turing parameter sets identified for each model, with two different sample sizes: light blue represents the $10^5$ sample, and dark blue represents the $10^6$ sample. In the 2-node network, we found 22 Turing parameter sets in the $10^5$ sample and 204 in the $10^6$ sample, representing approximately 0.0002% of the parameter space. For the 3-node network, we identified 57 and 548 Turing parameter sets in the respective samples (0.00055%), and in the 4-node network, 160 and 1,864 Turing parameters were found, representing around 0.0017% of the explored parameter space. These findings align with previous results SCHOLES2019243, indicating that larger networks exhibit a greater abundance of Turing patterns in parameter space. (C-E) Histograms of the highest (blue) and lowest (orange) ranked parameters based on the sign difference metric for the three models. The best parameters for the 2-node and 3-node networks, shown in (C) and (D), achieved a sign difference of approximately 0.9, indicating that most of the distribution lies on one side of the $x$-axis. In contrast, the best value for the 4-node network, shown in (E), is around 0.7, suggesting reduced controllability of the dispersion relation, and consequently, of pattern formation, as network complexity increases.(F-H) Barplots of the sign difference score for the different parameters of each model in decreasing order. We can observe that for the 2- and 3-node models (F and G, respectively), the top parameters all have a very high sign difference score, and then experience a sudden drop. This drop is shown in the inner plots, which show the sign difference score in decreasing order for all parameters. In contrast, for the 4-node model (H), we observe a lower initial value and a continuous decrease, with a wider drop at the end.
• Figure 3: Illustration of the quadratic approximation and wavevector lattices.(A) Sketch of the projection from the full reaction-diffusion system to the Swift-Hohenberg model, with key parameter $r$ acting as the switch between pattern formation and no pattern. (B) Application of the quadratic approximation for the dispersion relation. The original dispersion relation for a given model is shown in orange, while a quadratic fit (blue) around the maximum, $k_c^2$, leads to the Swift-Hohenberg (SH) equation. (C) The effect of the negative biharmonic operator and negative Laplacian operator on a cosine wave. On the left, the negative biharmonic smooths out the cosine wave, first reducing side peaks and eventually flattening the central peak to reach a homogeneous steady state. On the right, the negative Laplacian amplifies the wave's clumps, causing them to grow tighter, eventually tending toward a Dirac delta-like function in the limit as $t\rightarrow\infty$. The combination of these two opposing effects results in a stable pattern, as shown in the lower plot. (D) Lattices generated by varying numbers of wavevectors. A single wavevector produces stripes. With two wavevectors, a square lattice forms, supporting both stripes and squares. Three wavevectors produce a hexagonal lattice, where both hexagonal spots and stripes are stable.
• Figure 4: Model approximations, stability regions and pattern transition.(A) Different patterns produced by a saturating parameter set (the last in Table S1 in the SI), where the cubic nonlinearity is sufficient to saturate the amplitude equation. We show from left to right the original model, the Taylor-approximated (TA) model, the quadratic-diagonal (QD) model and the Swift-Hohenberg (SH) model. All approximations produce the same striped pattern for this parameter set. Orientation of stripes changes depending on initial conditions. (B) Stability diagram in $c_1-c_2$ space, obtained from the amplitude equation. The original parameters correspond to a projection of saturated parameters from Eq. \ref{['network2d']} to the Swift-Hohenberg equation. The different colors show the different areas where each pattern is stable. Blue corresponds to regions where only stripes are stable, green to regions where only spots are stable and red corresponds to unstable regions. Originally our system was at point a, where only stripes are stable. Investigating the effect of varying different parameter by using the projection, we found that increasing $k_A$ leads to the spots stable region (point b), whereas decreasing it leads to the unstable region (point c). (C) Simulations of the different parameter sets at points a, b and c for the 2-node network (above) and the Swift-Hohenberg model (below), resulting in the validation of the analytical prediction of pattern change.
• Figure 5: Potential experimental applications of the pattern switch.(A) An initial spatial distribution of the parameter $\mu_A$ allows the appearance of pattern formation in specific regions of space. (i) We let $\mu_A = \mu_A^{off}$, a value for which the dispersion is below zero, for the upper region of space, and $\mu_A = \mu_A^{on}$, a value which allows for pattern formation, for the lower region. (ii) We simulate the dynamics, allowing this parameter of the PDE to have the different spatial values and observe that pattern formation indeed is only produced in the region of space where $\mu_A = \mu_A^{on}$. (B) Workflow for using the pattern switch to test for different levels of concentrations of a parameter. (i) A genetic network where anhydrotetracycline (ATC) is used to suppress the inhibition from node $A$ to $C$ ($k_{CA}$ in the 3-node network introduced before). (ii) Checking the effect of ATC. As it only has a moderate effect, it is not the best pattern switch (although it still can be used). (iii) Determining three different conditions for which ATC can act as a pattern switch. (iv) Experimentally test of those conditions. Images in (iv) were taken from Fig. 1A in martina2023 with permission. Note the concentration of ATC was not computed using our workflow. Nevertheless, the data show some similarity to our approach.