Bifurcations and multistability in inducible three-gene toggle switch networks

This repository contains supporting code for:

• Rousseau RJ, Phillips R. Bifurcations and multistability in inducible three-gene toggle switch networks. Preprint found at arXiv:2509.24926v2 [q-bio.MN] https://doi.org/10.48550/arXiv.2509.24926v2 (2025)

The contents of the repository were used to generate the results found in figures throughout this paper, and explore the dynamics of three-gene toggle switch models for gene regulation as a function of allosteric effector concentration.

Layout

---

The repository is organized into several directories, each corresponding to the specified figure(s) in the paper:

Fig3

This directory contains the Matlab code Fig3_pact.m used to plot the probability of repressor activity as a function of inducer concentration (Figure 3).

Fig4

This directory contains the Matlab code Fig4_analyticthresh.m used to generate heatmaps tracking the analytically-derived bifurcation threshold beyond which the three-gene toggle switch can no longer have tristable dynamics, tracking change over a range of Hill coefficients and mRNA expression rates (Figure 4).

Fig5

This directory contains an input file Output_OG_Fig5.mat with the Mathematica output from solving for the fixed points of the baseline three-gene toggle switch. The Matlab code Fig5_inducer1_OG.m then uses this input to plot the dynamics of the three-gene toggle for increasing inducer concentration (Figure 5). Note that this code also generates a plot for Fig. 9(A) of the paper, and that Fig. 5(A)-(C) correspond to Fig. 13(A)-(B).

Fig6

This directory contains the input files Output_OG_Fig6_KI0_5.mat, Output_OG_Fig6_KI1_5.mat, and Output_OG_Fig6_KI5_0.mat (Mathematica output), which list fixed points derived for increasing dissociation constants within the MWC model for inducer binding. The Matlab code Fig6_inducer1diffKI_analysis.m then uses these inputs to plot the change in inducer concentration range for the intermediate tristable dynamic regime observed in the baseline three-gene toggle switch (Figure 6).

Fig7

This directory contains the input file Output_OG_Fig7.mat (Mathematica fixed point output) and the Matlab code Fig7_inducer1_nchange.m which uses this input to plot the change in the regions of multistability for a range of Hill coefficients (Figure 7).

Fig8

This directory contains the input file Output_OG_Fig8.mat (Mathematica fixed point output) and the Matlab code Fig8_inducer2_OG.m which uses this input to plot a phase diagram for the evolving stability profile of the three-gene toggle switch when tuned by two inducers (Figure 8).

Fig9

This directory contains the input files Output_OG_Fig9_v1thresh7_5.mat, Output_OG_Fig9_v1thresh5_3.mat, Output_OG_Fig9_v2thresh7_5.mat, and Output_OG_Fig9_v2thresh5_3.mat (Mathematica fixed point outputs). The Matlab code Fig9_inducer1c_perturb_v1v2.m then uses these inputs to plot the change in the two bifurcation thresholds of the single-inducer three-gene toggle under increasing perturbation away from symmetry among gene repression strengths in the circuit (Figure 9(B) and (C)).

Fig13

This directory contains the input file Output_OG_Fig13.mat (Mathematica fixed point output) and the Matlab code Fig13_inducer_selfactcomp_a2b0.m, which uses this input to plot the dynamics of the three-gene toggle switch in the presence of self-activation (Figure 13 (C) and (D)).

Fig14

This directory contains several input files Output_OG_Fig14_...mat containing fixed point output from Mathematica, evaluated for the three-gene toggle switch with self-activation at different ratios of basal expression level vs. activation strength. The Matlab code Fig14_bifurcincreaseb.m then uses these inputs to plot the dynamics of this three-gene toggle switch for increasing activation strength (Figure 14).

Fig15

This directory contains the Matlab code Fig15_inducer_selfactcompA_a2b2.m, which uses input Output_OG_Fig15A.mat to plot the dynamics of the self-activated three-gene toggle switch with competitive repressor-activator binding, where effector binding determines protein activity without altering its function (Figure 15 (A)). The directory also contains the Matlab code Fig15_inducer_selfactcompB_a2b2.m, which uses input Output_OG_Fig15B.mat to similarly plot dynamics for competitive binding where now effector binding determines protein function (Figure 15 (B)).

Fig16-18

This directory contains the Matlab code and .mat file inputs necessary to explore dynamics for the self-activated three-gene toggle switch with non-exclusive repressor-activator binding, where effector binding determines protein activity without altering its function (Figures 16-18). The workflow of this directory is:

1. Input: Output_OG_Figs_nonexcA.mat (Mathematica fixed point output)
2. Fig16_18_inducer_selfactnonexc_max2_fpstable.m: Matlab code that filters this input to isolate only stable fixed points, generating output fpstable_Figs_nonexcA.mat (included in the directory)
3. Fig16_18_inducer_selfactnonexc_max2.m: Matlab code that uses the output of the previous step to evaluate and plot bifurcations in the dynamics of the three-gene toggle switch.

Fig19-21

This directory contains the Matlab code and .mat file inputs necessary to explore dynamics for the self-activated three-gene toggle switch with non-exclusive repressor-activator binding, where effector binding determines protein function (Figures 19-21). The workflow of this directory is:

1. Input: Output_OG_Figs_nonexcB.mat (Mathematica fixed point output)
2. Fig19_21_inducer_selfactnonexcB_max2_fpstable.m: Matlab code that filters this input to isolate only stable fixed points, generating output fpstable_Figs_nonexcB.mat (included in the directory)
3. Fig19_21_inducer_selfactnonexcB_max2.m: Matlab code that uses the output of the previous step to evaluate and plot bifurcations in the dynamics of the three-gene toggle switch.