Interlinked GTPase cascades provide a motif for both robust switches and oscillators

TL;DR

This paper compares single-species GTPase switches and interlinked cascades, showing that cascades are more robust and versatile, capable of functioning as switches or oscillators depending on cofactor levels, with implications for cellular regulation.

Contribution

It introduces a thermodynamically consistent model demonstrating how interlinked GTPase cascades can function as robust switches and oscillators, expanding understanding of cellular signaling motifs.

Findings

Interlinked cascades are more robust than single-species switches.

Cascades can be tuned into oscillators by adjusting cofactor concentrations.

Both bistability and oscillations are achievable with the same feedback motif.

Abstract

GTPases regulate a wide range of cellular processes, such as intracellular vesicular transport, signal transduction, and protein translation. These hydrolase enzymes operate as biochemical switches by toggling between an active guanosine triphosphate (GTP)-bound state and an inactive guanosine diphosphate (GDP)-bound state. We compare two network motifs, a single-species switch and an interlinked cascade that consists of two species coupled through positive and negative feedback loops. We find that interlinked cascades are closer to the ideal all-or-none switch and are more robust against fluctuating signals. While the single-species switch can only achieve bistability, interlinked cascades can be converted into oscillators by tuning the cofactor concentrations, which catalyse the activity of the cascade. These regimes can only be achieved with sufficient chemical driving provided by GTP hydrolysis. In this study, we present a thermodynamically consistent model that can achieve bistability and oscillations with the same feedback motif.

Figures

Tables

Table 1. Table 1: Comparison of the switch qualities Q 𝑄 Q and hysteresis widths w hys subscript 𝑤 hys w_{\text{hys}} for different network motifs at a chemical affinity of Δ ​ μ / k B ​ T = 12 Δ 𝜇 subscript 𝑘 B 𝑇 12 \varDelta\mu/k_{\mathrm{B}}T=12 .

switch type		$Q$	$w_{\text{hys}}$
all-or-none switch		1	–
single-species switch	$x^{\ast }$	0.63	0.235
interlinked cascade	$x^{\ast }$	0.62	0.255
	$y^{\ast }$	0.84

Equations

$$X^{\ast }\underset{{\omega }_1^{-}}{\overset{{\omega }_1^{+}}{\rightleftharpoons }}X+{\mathrm{P}}_{\mathrm{i}},$$

$$I+2X^{\ast }+X+\mathrm{GTP}\underset{{\kappa }_{11}^{-}}{\overset{{\kappa }_{11}^{+}}{\rightleftharpoons }}I+3X^{\ast }+\mathrm{GDP},$$

$$\Delta \mu \equiv {\mu }_{\mathrm{GTP}}-{\mu }_{\mathrm{GDP}}-{\mu }_{{\mathrm{P}}_{\mathrm{i}}}=\Delta {\mu }_0+k_{\mathrm{B}}T\ln \frac{[\mathrm{GTP}]}{[\mathrm{GDP}][{\mathrm{P}}_{\mathrm{i}}]},$$

$$\frac{{\kappa }_{11}^{+}{\omega }_1^{+}}{{\kappa }_{11}^{-}{\omega }_1^{-}}={\mathrm{e}}^{\Delta \mu /k_{\mathrm{B}}T}.$$

$${\partial }_{t}P(n,t)$$

$$-({\alpha }_n+{\beta }_n)P(n,t).$$

$${\alpha }_n={\alpha }_n^{\kappa }+{\alpha }_n^{\omega }$$

$${\beta }_n={\beta }_n^{\kappa }+{\beta }_n^{\omega }$$

$$\frac{P_{\mathrm{ss}}(n)}{P_{\mathrm{ss}}(0)}=\prod _{i=0}^{n-1}\frac{{\alpha }_i}{{\beta }_{i+1}},$$

$$x^{\ast }\equiv \frac{1}{\Omega }\sum _{n=0}^{N}nP(n),$$

$$\frac{\mathrm{d}{x}^{\ast }}{\mathrm{d}t}$$

$$X^{\ast }$$

$$I+2X^{\ast }+X+\mathrm{GTP}$$

$$G_2^{+}+2X^{\ast }+Y+\mathrm{GTP}$$

$$Y^{\ast }$$

$$G_2^{+}+2Y^{\ast }+Y+\mathrm{GTP}$$

$$G_1^{-}+2Y^{\ast }+X^{\ast }$$

$$\frac{{\kappa }_{11}^{+}{\omega }_1^{+}}{{\kappa }_{11}^{-}{\omega }_1^{-}}$$

$$\frac{{\kappa }_{12}^{+}{\omega }_2^{+}}{{\kappa }_{12}^{-}{\omega }_2^{-}}$$

$$\frac{\mathrm{d}{x}^{\ast }}{\mathrm{d}t}$$

$$-\left({\omega }_1^{+}+[I]{\kappa }_{11}^{-}{x^{\ast }}^2+[G_1^{-}]{\omega }_{21}^{+}{y^{\ast }}^2\right)x^{\ast },$$

$$\frac{\mathrm{d}{y}^{\ast }}{\mathrm{d}t}$$

$$-\left({\omega }_2^{+}+[G_2^{+}]{\kappa }_{22}^{-}{y^{\ast }}^2+[G_2^{+}]{\kappa }_{12}^{-}{x^{\ast }}^2\right)y^{\ast }.$$

$$w_{\text{hys}}=\frac{[I_{+}]-[I_{-}]}{[I_{+}]+[I_{-}]},$$

$$0\leq Q\equiv\big{(}(1-\delta_{\text{off}})+(1-\delta_{\text{on}})+\delta_{\text{off}}^{\text{on}}\big{)}/3\leq 1,$$

$${\delta }_{\text{off}}(x^{\ast })\equiv \underset{\text{off}}{\max }{x}^{\ast }-\underset{\text{off}}{\min }{x}^{\ast }\ge 0,$$

$${\delta }_{\text{on}}(x^{\ast })\equiv \underset{\text{on}}{\max }{x}^{\ast }-\underset{\text{on}}{\min }{x}^{\ast }\ge 0,$$

$${\delta }_{\text{off}}^{\text{on}}(x^{\ast })\equiv \underset{\text{on}}{\min }{x}^{\ast }-\underset{\text{off}}{\max }{x}^{\ast }\le 1.$$

$$[I(t)]=[I_{0}]\big{(}1+\xi(t)\big{)},$$

$$\frac{{\text{Var}}_1}{[I_0]\sqrt{2D/\mathrm{d}t}}\simeq 4.0\cdot 10^{-4},$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{f}(\mathbf{x}),$$

$$\mathbf{f}(\mathbf{x})\approx \mathbf{f}({\mathbf{x}}_0)+\mathbf{J}({\mathbf{x}}_0)(\mathbf{x}-{\mathbf{x}}_0).$$

$${\partial }_x{f}_1(x,y)$$

$${\partial }_y{f}_1(x,y)$$

$${\partial }_x{f}_2(x,y)$$

$${\partial }_y{f}_2(x,y)$$

$$U(\mathbf{x})=-\ln P_{\mathrm{ss}}(\mathbf{x}).$$

$$\frac{\mathrm{d}{x}^{\ast }}{\mathrm{d}t}$$

Full text

Introduction

Interlinked cascades are a special class of network motifs that are fundamental components of cell signalling networks. They regulate the cellular response to an input signal using complex feedback structures and involve a multitude of proteins. Ideally, these molecular switches should have an all-or-none response, thus making them ultrasensitive Tu (2008); Ferrell and Ha (2014). Although this feature is essential for signal transducers, their steady state output is by definition sensitive to input fluctuations. Molecular switches can benefit from bistability which makes them robust against input noise Ferrell and Ha (2014). Interlinked cascades often consist of two species coupled through positive and negative feedback loops. Positive feedback loops can enhance the input signal as well as the bistability region, and make the switch more robust against fluctuations. Negative feedback loops are necessary for oscillations and enable a fast input response. The combination of positive and negative feedback is a highly versatile system as it facilitates rich dynamics such as bistability or oscillations, and is robust against input noise Pfeuty and Kaneko (2009); Kim et al. (2006, 2008); Kochanczyk et al. (2017).

Bistability can already be achieved by a minimal model containing only one species and an autocatalytic reaction as in Schlögl’s model Schlögl (1972); Vellela and Qian (2009), complex feedback mechanisms are not required. A natural question that arises is why interlinked cascades should be advantageous over single-species switches. In this study, we show that the interlinked cascade has an improved switching quality, which is closer to the ideal all-or-none switch compared to its single-species counterpart. Moreover, interlinked cascades can decrease the effect of input fluctuations Hooshangi et al. (2005); Zhang et al. (2007).

In this paper, we will focus on small GTPase proteins, which are key players in many cellular processes such as gene expression, membrane trafficking or vesicular transport Wittinghofer (2014a, b); Mizuno-Yamasaki et al. (2012). Their structure comprises three domains: a nucleotide-binding domain, where a guanosine triphosphate (GTP) or guanosine diphosphate (GDP) is lodged, a catalytic-binding domain, and an effector-binding domain Wittinghofer (2014a, b). The effector-binding domain is responsible for different functions in membrane transport, such as vesicle budding, movement, tethering, and fusion Behnia and Munro (2005); Hutagalung and Novick (2011). Small GTPases have a high affinity to nucleotides and are mostly found in complexes with GTP or GDP under physiological conditions. In the GTP-bound state, small GTPases are in an open activated state and can selectively bind to a specific set of effectors with high affinity. In the GDP-bound state, small GTPases are preferably closed and the system is considered inactivated as they have a low affinity to effectors.

Small GTPases rely on cofactors to regulate their activity. Cofactors are essential catalysers since the spontaneous GTPase rate is very low without them ($10^{-4}\,\mathrm{s^{-1}}$). Small GTPases can on the one hand form a complex with GTPase-activating proteins (GAPs), which greatly accelerate the GTPase rate (deactivation) up to $10\,\mathrm{s^{-1}}$ Gideon et al. (1992). On the other hand, they can bind to guanine nucleotide exchange factors (GEFs) that catalyse the dissociation of GDP in the inactive state (up to $10^{5}$-fold Klebe et al. (1995a)), followed by an activation with GTP. Hence, cofactors are essential to realise fast and reliable switching behaviour. In this paper, we show that these cofactor concentrations must be tuned appropriately to achieve bistability or oscillations. Recently, it was proposed that small GTPases are optimised to rapidly bind to their effectors and quickly release them Nguyen et al. (2017), which is ideal for tasks like defining identities within the membrane Barr (2013). Defects or mutations in the cofactors’ functions have been experimentally linked to reduced targeting Blümer et al. (2013); Cherfils and Zeghouf (2013). Small GTPases share many similarities with kinetic proofreading Hopfield (1974); Ninio (1975); Bennett (1979), Min proteins Huang et al. (2003); Loose et al. (2008); Halatek and Frey (2012); Xiong and Lan (2015); Wu et al. (2016); Denk et al. (2018) and the $70\,\mathrm{kDa}$ heat shock protein (Hsp70), which has an analogous ATP cycle and set of cofactors Kim et al. (2013); Finka et al. (2016); De Los Rios and Barducci (2014); Nguyen et al. (2017). Recent experiments have shown that the Hsp70-Hsp90 cascade can improve client folding yields Morán Luengo et al. (2018).

Small GTPases from different families are linked through feedback schemes to accomplish either bistability or oscillations. Interlinked cascades are designed to achieve bistability, they are found in many parts of the cell, such as the Rab5-Rab7 in the endosome or Rab/Ypt in the secretory pathway Ortiz et al. (2002); Pusapati et al. (2012); Rivera-Molina and Novick (2009); Nottingham et al. (2011); Suda et al. (2013); Rana et al. (2015). While each species has its specific set of cofactors and effectors, all of them share the same structure. The effector of the first species is the GEF of the second species, which catalyses the activation of the second system (GEF cascade). Once activated, the second species can bind to its effector and the GAP of the first species. When bound to the latter, it inactivates the first system (GAP cascade) as illustrated in Fig. 1. Errors when deactivating the previous species have been linked to trafficking defects Rivera-Molina and Novick (2009); Nottingham et al. (2011); Suda et al. (2013); Rana et al. (2015). Interestingly, small GTPases are also tuned to achieve oscillations. For example, Rho GTPases are drivers of cellular morphogenesis and serve as pattern to guide the formation of various cellular structures Goryachev and Pokhilko (2008); Bement et al. (2015); Goryachev et al. (2016); Graessl et al. (2017). In this paper, we show that interlinked small GTPases can be tuned to act as molecular switches or oscillators with the same feedback structure, thus providing a versatile motif. Fischer-Friedrich et al. have shown experimentally that the MinDE system, which shares a similar feedback structure, displays a transition from a switch to an oscillator as the cell length increases Fischer-Friedrich et al. (2010).

Small GTPases are GTP-driven machines that rely on the GTP hydrolysis cycle, which transforms a GTP into a GDP and an inorganic phosphate $\mathrm{P}_{\mathrm{i}}$. Under physiological conditions, an excess of GTP is maintained, driving the system far from equilibrium Traut (1994). Thermodynamics states that complex behaviour such as bistability or oscillations can only be achieved far from equilibrium Qian and Reluga (2005); Ge and Qian (2009); Xiong and Lan (2015); Barato and Seifert (2017); Nguyen et al. (2018); Herpich et al. (2018). Although extensive research has been carried out on small GTPase cascades, few studies have investigated the finite energy budget of the hydrolysis cycle. Apart from thermodynamically consistent models for a phosphorylation-dephosphorylation switch Qian and Reluga (2005); Ge and Qian (2009); Bishop and Qian (2010); Zhu et al. (2010); Ge and Qian (2011), very little is known about the thermodynamic properties of interlinked cascades. Thermodynamic consistency requires mass action kinetics, thus forbidding more complex activation functions such as sigmoidal, Michaelis-Menten, and Hill functions, which are primarily used in the literature to analyse interlinked cascades Stadtman and Chock (1977); Chock and Stadtman (1977); Gardner et al. (2000); Markevich et al. (2004); Del Conte-Zerial et al. (2008); Liu et al. (2011); Tiwari and Igoshin (2012); Huang et al. (2016). Moreover, the underlying biological reactions of these activations remain unclear in contrast to the fundamental reactions considered in the present study.

Model description

Thermodynamically consistent model for single-species switches

We consider a biological switch that consists of species $X$ in a volume $\Omega$. It can toggle between an inactive GDP-bound state $X$ and an active GTP-bound state $X^{*}$. Transitions between these two states occur through a hydrolysis reaction

[TABLE]

and a nucleotide exchange reaction

[TABLE]

where $\omega_{1}^{\pm},\kappa_{11}^{\pm}$ are transition rates shown in Fig. 2(a). We assume that the input $I$ acts as a catalyst by boosting the nucleotide exchange reactions like the GEF molecules for small GTPases. Second, we suppose that two active proteins are required to activate a third one and write $\kappa_{11}^{\pm}$ with double indices to indicate this self-feedback. We choose this effective reaction in order to achieve non-linear feedback, which is a necessary condition for bistability. In reality, a GEF only binds to one active protein to activate a second one Klebe et al. (1995b); Zhu et al. (2010); Stanley and Thomas (2016). This would normally result in one $X^{*}$ term on the left side of Eq. 2. However, Goryachev et al. Goryachev and Pokhilko (2008) have shown that a model for Rho GTPases with membrane interaction containing only bimolecular reactions can be simplified to an effective trimolecular reaction similar to Eq. 2. Note that we did not consider the equilibrium binding reactions in our model, they can be coarse-grained in a thermodynamically consistent way assuming timescale separation Esposito (2012); Wachtel et al. (2018).

A transition from $X$ to $X^{*}$ along the “+” direction of reaction (2) and back to $X$ along the “+” direction of reaction (1) does not change the state of the system, however, it turns one GTP into one GDP and a Pi. Such a complete cycle consumes a chemical energy

[TABLE]

where $\mu_{S}$ is the chemical potential of species $S=\mathrm{GTP},\mathrm{GDP},\mathrm{P}_{\mathrm{i}}$ and $[S]$ its concentration. The temperature is denoted by $T$ and $k_{\mathrm{B}}$ represents Boltzmann’s constant. The last equality is the approximation for an ideal solution, where $\varDelta\mu_{0}$ is a reference value. Physiological conditions are maintained by an excess of GTP, characterised by $\varDelta\mu>0$, and equilibrium corresponds to $\varDelta\mu=0$. The transition rates must satisfy the local-detailed balance condition Seifert (2011), which reads

[TABLE]

The state of the chemical reaction system (1) and (2) is determined by only one variable, the number of active proteins $n$. The time evolution of the probability distribution $P(n,t)$ to find the system in state $n$ at time $t$ is governed by the chemical master equation (CME), which reads

[TABLE]

Here, we define rate parameters as

[TABLE]

and introduce the total number of proteins $N$ as well as the input concentration $[I]$, which is given by the number of $I$ particles divided by the volume $\Omega$. For long times $t$ and constant $[I]$, the system reaches a non-equilibrium steady state (NESS), at which the left-hand side of Eq. 5 vanishes, leading to the condition $\alpha_{n-1}P_{\mathrm{ss}}(n-1)=\beta_{n}P_{\mathrm{ss}}(n)$. The analytic solution for the steady state probability distribution reads

[TABLE]

From the CME (5), we obtain the equation for the time evolution of the concentration of active proteins

[TABLE]

which reads in the deterministic limit ($N\to\infty,\Omega\to\infty$)

[TABLE]

As the total number of proteins $N$ is conserved, we set the sum of active and inactive concentrations to one (i.e. $x+x^{*}=1$), thus, setting $N/\Omega=1$. Note that the steady state concentration for zero input is simply given by $x^{*}([I]{=}0)=\omega_{1}^{-}/(\omega_{1}^{+}+\omega_{1}^{-})$. Having introduced the model for single-species switches, we consider interlinked cascades of different species $X,Y$ with positive self-feedback in the next section.

Thermodynamically consistent model for interlinked cascades

As in the model for the single-species switch, we assume that the input $I$ only affects species $X$ and catalyses its activation. We now introduce a positive feedback from the first $X^{*}$ to the second $Y^{*}$ component of the cascade and negative feedback from $Y^{*}$ to $X^{*}$, as shown in Fig. 3(a). The underlying biochemical reactions of the interlinked cascade are the following set of hydrolysis and nucleotide exchange reactions

[TABLE]

which are illustrated in Fig. 3(b). Note that the first two reactions of Eq. 10 are the same as Eqs. 1 and 2 of the single-species switch and that the additional feedback pathways are realised by the remaining four reactions. The molecules $G_{1}^{-}$ and $G_{2}^{+}$ can boost the hydrolysis and nucleotide exchange reactions, respectively, similar to the input $I$. In the case of small GTPases, $G_{1}^{-}$ would correspond to the GAP of species $X$ and $G_{2}^{+}$ to the GEF of species $Y$. Local-detailed balance requires

[TABLE]

The deterministic equations for the normalised active protein concentrations $x^{*}$ and $y^{*}$ are given by the following set of non-linear ordinary differential equations (ODEs)

[TABLE]

The results presented in the next section do not depend strongly on the choice of parameters. The explicit parameter values used in the simulations can be found in the Methods section.

Results

Bistability analysis of single-species switches

Deterministic analysis

We first consider the deterministic Eq. 9 and study its non-equilibrium steady state solutions for which the left-hand side vanishes. Depending on the input concentration $[I]$ and the chemical affinity $\varDelta\mu$, either one or three real solutions $x_{i}^{*}$, with $i=-,0,+$, are possible. Hence, one or three different steady states can occur as shown in Fig. 4(a). The solid lines correspond to stable steady states, the dashed lines to unstable steady states. If the chemical affinity is large enough, two stable steady states coexist and the system can exhibit hysteresis, which is a characteristic of bistability. It is interesting to note that an increase of the chemical affinity $\varDelta\mu$ enables an activation for which smaller input $[I]$ is sufficient. Bistability can only occur for non-vanishing chemical affinity, i.e., under non-equilibrium conditions. In Fig. 4(b), we show the phase diagram, where the chemical affinity and the input concentration are varied. The bistability region arises for a critical value of the chemical affinity $\varDelta\mu_{\mathrm{c}}/k_{\mathrm{B}}T\simeq 2.72$ and increases for larger affinities.

We analyse the bistability region in more detail by introducing the relative width of the hysteresis zone Qian and Reluga (2005)

[TABLE]

with upper $[I_{+}]$ and lower $[I_{-}]$ bounds of the hysteresis region for a given chemical affinity $\varDelta\mu$. As shown in Fig. 4(c), the minimal chemical affinity for hysteresis is $\varDelta\mu_{\mathrm{c}}/k_{\mathrm{B}}T\simeq 2.72$. For values $\varDelta\mu>\varDelta\mu_{\mathrm{c}}$, the relative hysteresis width growths until it saturates at its plateau value $w_{\text{hys}}\simeq 0.235$ for $\varDelta\mu/k_{\mathrm{B}}T=12$. From a theoretical perspective, a large hysteresis width has its assets and drawbacks. The major advantage is its increased robustness against input fluctuations, which is very beneficial in biological systems. However, it slows down the activation process because more input is needed for the entire activation in comparison with monostable activation functions like sigmoids. In this sense, an all-or-none activation with maximum distinction between the inactive and active state together with an appropriate hysteresis width is desirable.

Hence, as a measure of the performance of biochemical switches, we define the switch quality

[TABLE]

with the activation degrees in the “off” and “on” state as well as in the transition from “off” to “on” as

[TABLE]

The definitions enable a straightforward comparison with an all-or-none switch that would take on the optimum values $\delta_{\text{off}}=\delta_{\text{on}}=0$, $\delta_{\text{off}}^{\text{on}}=1$ and would achieve the maximum quality $Q_{\text{id}}=1$. Therefore, the switch quality measures the “distance” between the actual realisation and the ideal switch. A higher quality corresponds to a larger difference between the activations in the “off” and “on” state and is associated with an increased robustness, which are both biologically relevant. For all following investigations of the switch quality, the maximum value $\max_{\text{on}}x^{*}$ is taken from $[I]=3$ as shown in Fig. 4(a). It can be seen from the data in Fig. 4(c) that the switch quality of the single-species switch does not depend strongly on the chemical affinity and slightly increases with increasing $\varDelta\mu$ until it reaches its plateau value of $Q_{x}\simeq 0.63$ at $\varDelta\mu/k_{\mathrm{B}}T=12$.

Comparison of stochastic and deterministic description

In this section, we study the occupation of the two stable steady states in the bistability region and compare deterministic and stochastic descriptions. The probability distribution $P(x^{*})$ is shown in Fig. 5(a) for different input concentrations. In an experimental realisation of such switches, the protein concentration will stay in one of the bistable solutions for a very long time. In Fig. 5(b), we compare the active concentrations $x^{*}_{\text{det},i}$ for the different solutions of Eq. 9 with an average concentration $x^{*}_{\text{CME}}$ computed from the CME via Eq. 8, which is plotted for different system sizes $\Omega=25,50,100,1000$. Increasing the system size has two major advantages. First, the concentration in the active state increases, as expected. Second, the transition region from the lower stable solution to the upper one, in which apparent bistability is observed, decreases with increasing system size. For a system size $\Omega=1000$, the average of the stochastic concentration is in good agreement with the deterministic concentrations apart from the metastable region and the apparent bistability region is small. Ge and Qian Ge and Qian (2009) have shown that the chemical master equation converges to the deterministic solution with the lowest non-equilibrium potential. Note that the descriptions by the chemical master equation and the deterministic equation differ in the order in which the limits of time ($t\to\infty$) and volume ($\Omega\to\infty$) are taken. The CME considers the limits $\lim_{\Omega\to\infty}\lim_{t\to\infty}$, i.e., the time limit first and the volume limit second, whereas the deterministic equation takes the limits in opposite order, $\lim_{t\to\infty}\lim_{\Omega\to\infty}$.

Advantages of interlinked cascades over single-species switches

Multiple non-equilibrium phases

The phase space of the interlinked cascade is much richer than the phase space of the single-species switch. While only two phases (monostability and bistability) exist for the single-species switch, the interlinked cascade is additionally able to exhibit oscillation, excitability and tristability. We provide a brief overview of the most relevant nonequilibrium phases in this main text, detailed information can be found in Supplementary Note 1. Oscillation is characterised by the existence of a limit cycle, on which the trajectory circles around a source, which is the sole steady state. If a sink and a saddle arise additionally to the source, the phase is called excitability since the system can be excited to a quasi steady state. Additionally, a second type of bistability arises, which differs from the former bistability regime by an additional source and one further saddle.

The phase diagrams of the interlinked cascade for different negative feedback strengths $[G_{1}^{-}]=0.01,1.0,3.0$ are presented in Fig. 6. As a first result, there occur phase transitions from bistability to excitability and oscillation with increasing negative feedback strength from $[G_{1}^{-}]=0.01$ to $[G_{1}^{-}]=3.0$. Hence, by varying the experimentally accessible concentration of $G_{1}^{-}$ molecules, the network motif can be tuned to achieve both bistability and oscillations. Second, the diagrams reveal that the chemical affinity switches between different non-equilibrium phases.

Moreover, our model is able to achieve regions of tristability, where three stable steady states coexist. A detailed analysis in Supplementary Note 2 reveals that tristability is not an additional regime of the interlinked cascade as the probability of the third stable steady state is very low, making an experimental observation impossible. Consequently, the behaviour of the tristable system is quite similar to the bistable one.

All-or-none activation

The activating effect of the chemical affinity becomes apparent by studying the bifurcation diagrams of Fig. 7(a),(b). In these figures, the active concentrations $x^{*}$ and $y^{*}$ are plotted as functions of the input $[I]$ for chemical affinities (a) $\varDelta\mu/k_{\mathrm{B}}T=6$ and (b) $\varDelta\mu/k_{\mathrm{B}}T=12$. The activation of the second component changes drastically with increased chemical affinity and leads to a better switching behaviour.

We quantify the above statement by studying the switch quality for the interlinked cascade. The qualities $Q_{x}$ and $Q_{y}$ of both components are shown in Fig. 7(c) as function of the chemical affinity and allow for a direct comparison with the single-species switch. In contrast to the latter, the switch quality depends considerably on the chemical affinity for $\varDelta\mu/k_{\mathrm{B}}T\lesssim 10$. For larger affinities, $Q_{x}$ and $Q_{y}$ approach their plateau values, which are presented in Table 1 together with the corresponding values for the single-species and the all-or-none switch. It is apparent from these data that the second component of the interlinked cascade leads to a significantly better performance with a switch quality of $Q_{y}\simeq 0.84$ at $\varDelta\mu/k_{\mathrm{B}}T=12$. The associated hysteresis width and the quality of the first component remain almost unaffected with $w_{\text{hys}}\simeq 0.255$ and $Q_{x}\simeq 0.62$. Besides, the data reveal the importance of the chemical affinity since the switch quality is substantially improved for values $\varDelta\mu/k_{\mathrm{B}}T\gtrsim 7$ . Overall, the activation performance of the interlinked cascade is relatively close to the ideal all-or-none switch for high enough chemical affinities. The observed kinks in the switch qualities at $\varDelta\mu/k_{\mathrm{B}}T\simeq 5.5,6.5$ as well as the kink in the relative hysteresis width at $\varDelta\mu/k_{\mathrm{B}}T\simeq 8.1$ can be attributed to the onsets of the lower excitability region, bistability type 2, and the upper excitability region, respectively (see Fig. 6(b)).

Robustness against noisy input

The results presented thus far provide insights into steady state properties of biological cascade switches. In this paragraph, we study the dynamical behaviour of the system. To this end, we consider a noisy input function to obtain a realistic biological signal Tian et al. (2009); Zhang et al. (2007). The noisy input is described by

[TABLE]

where $\xi(t)$ is normally distributed white noise with zero mean $\braket{\xi(t)}=0$ and correlations $\braket{\xi(t)\xi(t^{\prime})}=2D\delta(t-t^{\prime})$. Here, the parameter $D$ denotes the noise intensity and $\braket{\cdot}$ denotes an ensemble average. This stochastic input is used in the ordinary differential Eq. 9 for the single-species switch and Eq. 12 for the interlinked cascade. In order to keep the analysis as simple as possible, we consider input fluctuations solely in the monostable regime. The response of the single-species switch and interlinked cascade to the noisy input is shown in Fig. 8 for $[I_{0}]=2.0$, $D=0.005$, and $\varDelta\mu/k_{\mathrm{B}}T=12$. The fundamental difference between the two switches is that the second component of the interlinked cascade is much less noisy than the single-species switch. This enhanced robustness against fluctuations is a striking advantage of the interlinked cascade as it increases the distinction between active and inactive state by reducing the fluctuations. Hence, the interlinked cascade acts as a low-pass filter Hooshangi et al. (2005); Zhang et al. (2007).

Having discussed the enhanced robustness qualitatively, we will now quantify our statements by computing the variances of the single-species switch ($\text{Var}_{1}$) and of the second component of the interlinked cascade ($\text{Var}_{2}^{y}$). The values for $[I_{0}]=2.0$, $D=0.005$, $t_{\mathrm{f}}=10^{6}$, $\mathrm{d}t=0.01$, and $\varDelta\mu/k_{\mathrm{B}}T=12$ are

[TABLE]

which leads to a variance ratio of $\text{Var}_{2}^{y}/\text{Var}_{1}\simeq 0.046$. This small number reveals the enhanced robustness against noisy input of the second component which makes it easier for such biological switches to react to input signals.

We now turn to the behaviour of the variance ratio of the interlinked cascade $\text{Var}_{2}^{y}/\text{Var}_{2}^{x}$ as a function of the chemical affinity $\varDelta\mu$. Figure 9 provides the results obtained from simulations for varied input concentrations $[I_{0}]$ (top row) and negative feedback strengths $[G_{1}^{-}]$ (bottom row). Major differences between the input concentrations $[I_{0}]\leq 1.4$ and $[I_{0}]\geq 2.0$ can be seen from the data in Fig. 9(a) and (b). There is a clear trend of increasing variance ratio up to $[I_{0}]\simeq 1.4$ along with almost no activation within this regime. For $[I_{0}]\simeq 2.0$, two peaks emerge and the associated activation functions are remarkably stronger. Further increasing the input leads to a decrease of the peak heights and the variance ratio. The peak positions appear to correlate with the maximum increase of $x^{*}$ and $y^{*}$, which occur at different chemical affinity $\varDelta\mu$. The negative feedback strength $[G_{1}^{-}]$ has a strong impact on the variance ratio for $\varDelta\mu/k_{\mathrm{B}}T\gtrsim 3$ (Fig. 9(c)). This is the most striking aspect of the data as it discloses a trade-off between robustness and oscillations. Increasing the negative feedback strength enables the interlinked cascade to perform oscillations, however, at the expense of a large variance ratio. Finally, increasing the chemical affinity above $\varDelta\mu/k_{\mathrm{B}}T\gtrsim 6$ decreases the fluctuations, as recently observed in biochemical oscillators by del Junco and Vaikuntanathan (2018).

Discussion and summary

In this study, we develop a thermodynamically consistent model for biochemical cascade switches that solely considers fundamental biochemical reactions. The use of trimolecular reactions may at first glance seem artificial. Nonetheless, our choice can be justified by the following reasons. First, we choose effective reactions as we are interested in a simple thermodynamically consistent network achieving both bistability and oscillations. It is well known that for an interlinked model as simple as ours, where the species have two states each, bistability or oscillations can only be achieved using trimolecular reactions if we assume mass action kinetics Schlögl (1972); Novák and Tyson (2008). In that sense, our interlinked cascade can be seen as a combination of Schlögl’s model Schlögl (1972) and the Brusselator Nicolis and Prigogine (1977). Second, it turns out that trimolecular reactions can be the result of several bimolecular reactions and additional species states. These states can then be coarse-grained into an active and inactive state using timescale separation, as shown for the Rho GTPase Goryachev and Pokhilko (2008). In addition, minimal models with only bimolecular reactions that achieve bistability or oscillations have been identified Wilhelm (2009). However, they require additional species compared to their usual trimolecular counterparts, which is equivalent to having additional protein states.

Generally, switches must operate far from equilibrium ($\varDelta\mu\neq 0$) in order to observe interesting regimes such as bistability, excitability, and oscillation. Our model for a single-species switch is able to exhibit bistability if the chemical affinity is high enough. Bistability makes the switch robust against input fluctuations and is accompanied by a hysteresis behaviour. As a measure of the performance of switches, we define the switch quality using the activation degrees in the ”off“ and ”on“ states as well as in the transition between them. An ideal all-or-none switch corresponds to the maximum switch quality of 1 and exhibits hysteresis. We show that the quality of the single-species switch is considerably lower than the optimum and slightly increases with higher chemical affinity.

Adding a second species and connecting them with positive and negative feedback relations leads to a more sophisticated model. The phase space of the interlinked cascade is richer and involves additional regimes, such as excitability, oscillation, and tristability. The interlinked cascade has several advantages over the single-species switch. First, the switch quality of the interlinked cascade is much higher than the one of the single-species switch and is closer to an ideal all-or-none switch. If we consider an input signal with noise, the interlinked cascade is shown to be much less noisy. Hence, it has an increased robustness against fluctuating signals and the distinction between the ”off“ and ”on“ state becomes clearer. Berg et al. have also shown that fluctuations decrease the sharpness of the activation curve, providing further evidence that the interlinked cascade is a better motif for biochemical switches Berg et al. (2000). The evidence presented in this study suggests that there exists a trade-off between switch quality, hysteresis width, and chemical affinity, which has to be balanced in biochemical switches.

The performance of the cascade could be further increased by adding more species to it. However, additional species may lead to a decrease of robustness as there is an optimal cascade length Detwiler et al. (2000); Thattai and van Oudenaarden (2002). For small GTPases, interlinked cascades have only been identified between two species so far and regarding three species, merely GEF cascades (without negative feedback) have been discovered Mizuno-Yamasaki et al. (2012). There is abundant room for further progress in examining switches with different feedback structures, see for example Kim et al. (2006, 2008); Tsai et al. (2008). Since these studies do not assume mass action kinetics, they lack thermodynamic consistency. We have explored variations of the interlinked cascade Ehrmann (2018). By replacing the negative feedback loop with a positive one, we obtain another thermodynamically consistent switch. We can still achieve bistability and tristability, but as a drawback of this motif, the regimes of excitability and oscillation are no longer accessible. In case of a double negative feedback motif, it is very difficult to find a biologically relevant parameter range and most of the settings lead to an irreversible activation curve. These feedback structures could be relevant for systems with complex activation functions like gene expression Tian et al. (2009). Nonetheless, our simple model presented here is sufficient to explain the versatility of interlinked small GTPases, which solely comprise positive-negative feedback structures.

Small GTPases have cofactors that are essential to achieve bistability or oscillations. Our model indicates that increasing the negative feedback, i.e., the GAP concentration, converts biochemical switches into oscillators. It would be very interesting to verify this prediction experimentally. Furthermore, a future study investigating how small GTPases respond to a heat shock, where the concentration ratio of GTP over GDP decreases, would be interesting Traut (1994). Besides, we expect the catalytic power of the cofactors to vary during stress to maintain the function of the cascade, as observed in the Hsp70 system Nguyen et al. (2017). An issue that was not addressed in this study is the membrane-cytoplasmic reaction and a third cofactor, namely the guanosine nucleotide dissociation inhibitor (GDI). It binds to the inactive GDP-bound state and transports it across the membrane Mizuno-Yamasaki et al. (2012); Wittinghofer (2014a, b). From a dynamical perspective, the GDIs can further inactivate the system in the absence of feedback and can control the total concentration of active species. Assuming fast kinetics Del Conte-Zerial et al. (2008), we do not expect major qualitative differences when including GDIs to our model.

We emphasise that our model is, to the best of our knowledge, a simple thermodynamically consistent model that demonstrates the biological functions of interlinked cascades, such as bistability and oscillations. Although the focus of this study is on GTPase cascades, our results are general and can be applied to a broader class of systems. In a future step, we could incorporate our cascade into a larger signalling network to investigate how the switch quality affects biological processes.

Methods

The explicit values for the parameters of the single-species switch and interlinked cascade are chosen in such a way to satisfy biological requirements as well as to get reasonable bifurcation and phase diagrams. While most of the parameters are given in natural units ($1.0\,\mathrm{s^{-1}}$), the rate $\kappa_{22}^{+}$ must be larger to observe interesting non-equilibrium regions. However, it must be smaller than $4.0\,\mathrm{s^{-1}}$ to achieve reversible activation curves. We choose $\kappa_{22}^{+}=3.9\,\mathrm{s^{-1}}$ to enlarge the nonmonostable region and we take $\kappa_{11}^{+}=\kappa_{22}^{+}$ to keep the system as simple as possible. The rates $\omega_{1}^{-}$ and $\omega_{2}^{-}$ must be smaller than their counterparts $\omega_{1}^{+}$ and $\omega_{2}^{+}$ in order to achieve a dominant inactivation in the hydrolysis reactions. The transition rates are in the same order of magnitude as in Pfeuty and Kaneko (2009). We set the parameters to $\omega_{1}^{+}=1.0\,\mathrm{s^{-1}}$, $\omega_{1}^{-}=0.05\,\mathrm{s^{-1}}$, $\omega_{2}^{+}=1.0\,\mathrm{s^{-1}}$, $\omega_{2}^{-}=0.01\,\mathrm{s^{-1}}$, $\kappa_{11}^{+}=3.9\,\mathrm{s^{-1}}$, $\kappa_{22}^{+}=3.9\,\mathrm{s^{-1}}$, $\kappa_{12}^{+}=1.0\,\mathrm{s^{-1}}$, $\omega_{21}^{+}=1.0\,\mathrm{s^{-1}}$, $[G_{1}^{-}]=1.0$, $[G_{2}^{+}]=1.0$, $D=0.005$, $\mathrm{d}t=0.01$. The remaining four reaction rates $\kappa_{11}^{-}$, $\kappa_{22}^{-}$, $\kappa_{12}^{-}$, $\omega_{21}^{-}$ are fixed by the local-detailed balance conditions of Eq. 11. Note that the negative feedback strength $[G_{1}^{-}]$ is changed in Figs. 6 and 9.

The statics and dynamics of the system are analysed by deterministic mass action kinetics as well as chemical master equations and chemical Langevin equations. A fourth order Runge-Kutta method Strogatz (2000) is used to solve the ordinary differential equations numerically with an integration time step $\mathrm{d}t=0.01$ and a final integration time $t_{\mathrm{f}}$. We used Mathematica (Wolfram Research, Inc., Champaign, Illinois), Python (Python Software Foundation), and C++ (Standard C++ Foundation) for the simulations.

Data accessibility

All data essential to evaluate and establish the conclusions of this study are presented in the paper itself. Computer codes can be obtained from the author upon request.

Author contributions

A.E. performed research. A.E., B.N., and U.S. designed research, discussed the results and wrote the paper.

Competing interests

The authors declare no conflict of interest.

Supplementary Information

**Interlinked GTPase cascades provide a motif for both robust switches and oscillators

** Andreas Ehrmann, Basile Nguyen, and Udo Seifert

II. Institut für Theoretische Physik, Universität Stuttgart, 70550 Stuttgart, Germany

(Dated: March 16, 2024)

Supplementary Note 1. Linear stability and phase plane analysis

Linear stability analysis

The stability of the fixed points of the system is determined by a linear stability analysis in which we consider the following set of non-linear differential equations

[TABLE]

with the concentration vector $\bm{x}=(x_{1},\dots,x_{M})^{T}$ of $M$ different protein species and the equation vector $\bm{f}(\bm{x})=(f_{1}(\bm{x}),\dots,f_{M}(\bm{x}))^{T}$. In linear stability analysis, such non-linear systems are approximated around each fixed point $\bm{x}_{0}$ by a linear system via

[TABLE]

To this end, one determines the fixed points $\bm{x}_{0}$ by solving $\bm{f}(\bm{x}_{0})=0$, computes the Jacobian $\mathbf{J}(\bm{x})\equiv\partial_{\bm{x}}\bm{f}(\bm{x})$ and its eigenvalues $\lambda_{i}$ at each fixed point $\bm{x}_{0}$. The stability of each fixed point is based on the real parts of its eigenvalues.

In one dimension, as for the single-species switch in Eq. 9 of the main text, the analysis is relatively simple and the stability of the fixed points is represented by the potential $U(x)=-\partial_{x}f_{1}(x)$. For two-dimensional systems as in Eq. 12 of the main text for the interlinked cascade, stable fixed points arise if both eigenvalues have negative real parts $\operatorname{Re}(\lambda_{1})<\operatorname{Re}(\lambda_{2})<0$ and are called sinks. In all other cases, the fixed point is unstable. For one positive and one negative value $\operatorname{Re}(\lambda_{1})<0<\operatorname{Re}(\lambda_{2})$, it is called saddle. For a source, both are positive $\operatorname{Re}(\lambda_{1})>\operatorname{Re}(\lambda_{2})>0$. The entries of the Jacobian for the interlinked cascade are given by

[TABLE]

Phase plane analysis of the interlinked cascade

As the coupled non-linear ODEs (Eq. 12 of the main text) are not analytically solvable, we perform a phase plane analysis to examine the fixed points of the system and determine all possible classes of dynamical behaviour through an exhaustive search over parameters. In Fig. S1, a typical phase plane is presented with parameter values $[I]=1.5$ and $\varDelta\mu/k_{\mathrm{B}}T=12$. The two nullclines $\mathrm{d}x^{*}/\mathrm{d}t=0$ and $\mathrm{d}y^{*}/\mathrm{d}t=0$ are represented by the blue lines and the circles represent its non-equilibrium steady states. The two stable fixed points are marked as black filled circles, while the three unstable fixed points are marked as black open circles for the saddles and a grey open circle for the source. When the input $[I]$ is varied, the $x^{*}$-nullcline primarily changes, whereas the $y^{*}$-nullcline is nearly unaffected.

The number of sinks $l$ (stable), saddles $m$ (unstable) and sources $n$ (unstable) characterise the different non-equilibrium phases of the system, which are denoted by a triplet $(lmn)$. Seven different phases can be identified by an exhaustive search over parameters: monostability (100), bistability type 1 (210), bistability type 2 (221), excitability (111), oscillation (001), tristability type 1 (320), and tristability type 2 (331). These phases are connected by three types of bifurcations, namely saddle-node, Hopf, and saddle-node infinite period (SNIPER) bifurcations. A saddle-node bifurcation occurs when the system goes through a bistability region. A Hopf bifurcation occurs when the system performs a transition between a monostability (100) and oscillation (001) region and a SNIPER bifurcation occurs in case of a transition between an excitability (111) and oscillation (001) region.

Oscillation is characterised by the existence of a limit cycle, on which the trajectory circles around the source, which is the sole steady state. This is shown in Fig. S2(a),(b). If a sink and a saddle arise additionally to the source, the phase is called excitability due to its activation properties. If the initial concentrations $(x_{0}^{*},y_{0}^{*})$ are beyond the sink’s attraction region, the system can be excited to a quasi-stationary manifold. As presented in Fig. S2(c),(d) the trajectory follows the unstable manifold of the saddle and then returns along the stable manifold of the sink. The trajectory in the excitable regime looks very similar to one oscillation phase of the trajectory in the oscillatory regime. The driving mechanism of this behaviour could be a heteroclinic orbit which becomes a limit cycle after the SNIPER bifurcation. Additionally, a second type of bistability arises, which differs from the former bistability regime by an additional source and one further saddle. The same holds for the two different tristability regimes and the phases with a source and an additional saddle are called type 2, while the others are called type 1.

Supplementary Note 2. Examining tristability with the effective potential

For two-dimensional systems, it is not possible to calculate an explicit potential function $U(x,y)$ because our interlinked model is not a gradient system, meaning that $\partial_{y}f_{1}(x,y)\neq\partial_{x}f_{2}(x,y)$ Verd et al. (2014); Strogatz (2000). However, an effective potential can be calculated according to the steady state probability distribution

[TABLE]

We consider a stochastic version of the system in order to determine the probability distribution by studying the chemical Langevin equations of the interlinked cascade, which are given by Gillespie (2000)

[TABLE]

with independent Gaussian white noise

[TABLE]

where $i,j\in\{1,\dots,12\}$. We omit the explicit time dependence of $x^{*}(t)$, $y^{*}(t)$, and $\xi_{i}(t)$ in the ODEs for simplicity reasons. In the simulations, we set the noise intensity to $D=0.005$. Note that in the thermodynamic limit $\Omega\to\infty$, the deterministic equation 12 of the main text is recovered.

By increasing the reaction rate from $\kappa_{12}^{+}=1.0\,\mathrm{s^{-1}}$ in Fig. 5(a) of the main text to $\kappa_{12}^{+}=3.9\,\mathrm{s^{-1}}$, the new regimes of tristability type 1 and type 2 are observable, which is shown in the phase diagram of Fig. S3(a). The associated activation function for $\varDelta\mu/k_{\mathrm{B}}T=12$ is shown in Fig. S3(b). Note that the switch quality $Q_{y}$ is even better as for the regime studied in the main text (if we neglect the improbable middle steady state). This supports our statements about the improved performance of the interlinked cascade. Theoretically, the system is able to exhibit tristability for this parameter set. In order to examine the occupation probabilities of these three stable steady states, the effective potential $U(x^{*},y^{*})=-\ln P_{\mathrm{ss}}(x^{*},y^{*})$ determined by the chemical Langevin Eq. S5 is studied. As presented in the contour plot in Fig. S3(c), the potential of the third NESS in the top left corner is four levels higher than the potential of the other two. Correspondingly, the probability is four orders of magnitude smaller, making an experimental observation impossible. Hence, the behaviour of the tristable system is quite similar to the bistable one. Note that the first component in Fig. S3(b) exhibits tristability in the same range as the second component although there is no evidence in the figure as the additional stable and unstable steady states have the same $x^{*}$-value at the lower stable steady state as shown in Fig. S3(c).

We emphasise that the statements in this paragraph hold for all parameter regimes studied in this paper. However, the existence of another parameter range, in which genuine tristability may occur is conceivable but is not an aim of this research.