Stability analysis of a signaling circuit with dual species of GTPase switches
SSTABILITY ANALYSIS OF A SIGNALING CIRCUIT WITH DUAL SPECIESOF GTPASE SWITCHES
LUCAS M. STOLERMAN , PRADIPTA GHOSH , , ∗ , AND PADMINI RANGAMANI ∗ A BSTRACT . GTPases are molecular switches that regulate a wide range of cellular pro-cesses, such as organelle biogenesis, position, shape, and function, vesicular transportbetween organelles, and signal transduction. These hydrolase enzymes operate by tog-gling between an active (”ON”) guanosine triphosphate (GTP)-bound state and an in-active (”OFF”) guanosine diphosphate (GDP)-bound state; such a toggle is regulatedby GEFs (guanine nucleotide exchange factors) and GAPs (GTPase activating proteins).Here we dissect a network motif between monomeric (m) and trimeric (t) GTPases as-sembled exclusively in eukaryotic cells of multicellular organisms. To this end, we de-velop a system of ordinary differential equations in which these two classes of GTPasesare interlinked conditional to their ON/OFF states within a motif through feedforwardand feedback loops. We provide explicit formulas for the steady states of the system andperform classical local stability analysis to systematically investigate the role of the dif-ferent connections between the GTPase switches. Interestingly, a feedforward from theactive mGTPase to the GEF of the tGTPase was sufficient to provide two locally stablestates: one where both active/inactive forms of the mGTPase can be interpreted as havinglow concentrations and the other where both m- and tGTPase have high concentrations.Moreover, when a feedback loop from the GEF of the tGTPase to the GAP of the mGT-Pase was added to the feedforward system, two other locally stable states emerged, bothhaving the tGTPase inactivated and being interpreted as having low active tGTPase con-centrations. Finally, the addition of a second feedback loop, from the active tGTPase tothe GAP of the mGTPase, gives rise to a family of steady states that can be parametrizedby a range of inactive tGTPase concentrations. Our findings reveal that the coupling ofthese two different GTPase motifs can dramatically change their steady state behaviorsand shed light on how such coupling may impact information processing in eukaryoticcells.
Date : September 18, 2020.
Key words and phrases.
GTPases, Biochemical switches, Network motifs, Stability analysis. Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla CA 92093. Department of Medicine,University of California, San Diego, La Jolla, CA 92093. Department of Cellular and Molecular Medicine, University of California, San Diego, La Jolla, CA 92093. Moores Comprehensive Cancer Center, University of California, San Diego, La Jolla, CA 92093. ∗ To whom correspondence should be addressed. e-mail: [email protected]; [email protected] a r X i v : . [ q - b i o . S C ] S e p L.M. STOLERMAN, P. GHOSH, AND P. RANGAMANI1. Introduction 32. Model Development 52.1. Assumptions. 52.2. Governing Equations. 53. Mathematical analysis and results 83.1. Feedforward connection: Recruitment of tGEF by active mGTPases ( mG ∗ → tGEF ) . 93.2. Feedforward connection mG ∗ → tGEF with feedback loop tGEF → mGAP : Recruitment of tGEFby active mGTPases and tGEF colocalization with mGAP. 133.3. Feedforward connection mG ∗ → tGEF with feedback loops tGEF → mGAP and tG ∗ → mGAP : Recruitment of tGEF by active mGTPases, tGEF colocalization with mGAP, and activationof mGAP by active tGTPases. 173.4. Numerical Simulations. 214. Discussion 265. Acknowledgments 29References 29Appendix A. Proof of Proposition 3.1 33Appendix B. Proof of Theorem 3.2 35Steady states. 35Local Stability Analysis. 38Appendix C. Proof of Theorem 3.3 40ContentsTABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 3
1. I
NTRODUCTION
Each eukaryotic cell has a large number of GTP-binding proteins (also called
GTPases or G-proteins); they are thought to be intermediates in an extended cellular signalingand transport network that touches on nearly every aspect of cell function [1, 2, 3]. Oneunique feature of GTPases is that they serve as biochemical switches that exist in an‘OFF’ state when bound to a guanosine diphosphate (GDP), and can be turned ‘ON’when that GDP is exchanged for a guanosine triphosphate (GTP) nucleotide. [1, 4].Turning the GTPase ‘ON’ is the key rate limiting step in the activation-inactivation pro-cess, requires an external stimulus, and is catalyzed by a class of enzymes called guaninenucleotide exchange factors (GEFs) [5]. G proteins return to their ‘OFF’ state when thebound GTP is hydrolyzed to guanosine diphosphate (GDP) via an intrinsic hydrolaseactivity of the GTPase; this step is catalyzed by GTPase-activating proteins (GAPs) [6].Thus, GEFs and GAPs play a crucial role in controlling the dynamics of the GTPaseswitch and the finiteness of signaling that it transduces [7, 8, 9, 10]. Dysregulation ofGTPase switches has been implicated in cellular malfunctioning and is commonly en-countered in diverse diseases [11, 12, 13, 14]. For example, hyperactivation of GTPases[15, 16] is known to support a myriad of cellular phenotypes that contribute to aggres-sive tumor traits [17, 18]. Such traits have also been associated with aberrant activity ofGAPs [15] or GEFs [19, 20, 21, 22, 23, 24]. These works underscore the importance ofGTPases as vital regulators of high fidelity cellular communication.There are two distinct types of GTPases that gate signals: small or monomeric (m) andheterotrimeric (t) GTPases. mGTPases are mostly believed to function within the cell’sinterior and are primarily concerned with organelle function and cytoskeletal remodeling[25, 26, 27]. tGTPases, on the other hand, were believed to primarily function at the cell’ssurface from where they gate the duration, type and extent of signals that are initiated byreceptors on the cell’s surface [28, 29]. These two classes of switches were believed tofunction largely independently, until early 1990’s when tGTPases were detected on in-tracellular membranes, e.g., the Golgi [29, 30], and studies alluded to the possibility thatthey, alongside mGTPases, may co-regulate organelle function and structure [31]. But itwas not until 2016 that the first evidence of functional coupling between the two switches– m- and tGTPases– emerged [32]. Using a combination of biochemical, biophysical,structural modeling, live cell imaging, and numerous readouts of Golgi functions, it wasshown that m- and tGTPase co-regulate each other on the Golgi [32]. mGTPase (Arffamily) must be turned ‘ON’ for it to engage with a GEF (GIV/Girdin) for tGTPase,and the latter subsequently triggers the activation of a tGTPase (of the Gi sub-family).Upon activation, the tGTPase enhances a key GAP for mGTPases (i.e., ArfGAP2/3) thatturns ‘OFF’ the mGTPase, thereby terminating mGTPase signaling. This phenomenonof co-regulation between the two classes of GTPases was shown to be critical for main-taining the finiteness of signaling via both species of GTPase on the membrane, which
L.M. STOLERMAN, P. GHOSH, AND P. RANGAMANI was critical for maintaining Golgi shape and function (trafficking within the secretorypathway). In doing so, this dual GTPase circuit was converting simple chemical signalsinto complex mechanical outputs (membrane trafficking). Because emerging evidencefrom protein-protein interaction networks and decades of work on both species of GT-Pases suggest that co-regulation through coupling between the GTPases is possible andlikely on multiple organellar membranes, it begs the question – what are the advantagesof two species of GTPase switches coupled through a closed feedback loop within a net-work motif as opposed to single switches? The answer to this question has not yet beenexperimentally dissected or intuitively theorized.Mathematical models of signaling networks have contributed significantly to our under-standing of how network motifs might function [33, 34, 35, 36, 37, 38]. Continuous-timedynamical systems, commonly represented by systems of ordinary differential equations(ODEs), are powerful tools for building rich and insightful mathematical models [39].For example, a comprehensive steady state analysis of ODE system helped frame theconcept of “zero-order ultrasensitivity” where large responses in the active fraction of aprotein of interest are driven by small changes in the saturated ratio of the enzymes [40].Similarly, modeling biochemical networks with dynamical systems also has revealed theexistence of bistable switches and biological oscillators within a feedback network archi-tecture [41, 42, 43]. The simple system proposed by Ferrell and Xiong [43] served as abasis for modeling cellular all-or-none responses, and hence, crucial for decision-makingwithin several signaling processes. Dynamical systems have also be used for mappingchemical reactions into differential equations; when numerically integrated (clustered)these systems can be used to predict the evolution of large reaction networks [44, 45, 46].For example, clustering methods have revealed the existence of recurrent structures, theso-called “network motifs’, the dynamics of which have been inferred with a Booleankinetic system of differential equations [46]. These models produced various dynamicfeatures corresponding to the motif structure, which allowed the understanding of under-lying biology (i.e., gene expression patters). Furthermore, studies using network ODEmodels have revealed that information is processed in cells through intricate connectionsbetween signaling pathways rather than individual motifs [47, 48, 49]. These findings ledto the so-called “learned behavior” hypothesis, on comprised information within intracel-lular biochemical reactions, was verified by analyzing extended signal duration, feedbackloops, and multiple thresholds of stimulation and signal outputs. Finally, from a systemsbiology modeling perspective, large systems of differential equations are usually hardto analyze, but when combined with experiments, they can give rise to quantitative pre-dictions [50, 51, 52], as exemplified by the work of Yen et al. [53]. Finally, dynamicalsystems modeling of biochemical networks can help understand how stochastic noiseis attenuated or amplified within cells [54, 55, 56, 57, 58], e.g., the trade-off betweensensitivity and noise reduction in any given network topology, and how some networktopologies could perform a dual function of noise attenuation and adaptation [58].
TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 5
Here we build mathematical models to investigate the dynamic properties of a coupledbiochemical GTPase switches. Beginning with the uncoupled GTPase switches (Fig.1A)as our starting point, we specifically sought to understand the stability features of a newnetwork motif (Fig.1B) to gain insights into any relative advantages of coupling twodistinct classes of GTPase switches. To this end, we obtained the steady states of thisnew network motif, to understand the input-output relationships. Then we studied thedynamic behavior under small perturbation around these steady states using local sta-bility analysis [59, 60]. We investigated the different feedforward and feedback loopsbetween these two motifs, all representing the observed biochemical and biophysicalevents during signal transduction (Fig.1C). To gain biological intuition, we investigatedthe appropriate space of initial conditions that guarantee the existence of the differentsteady states, thus obtaining insight to the different regimes of the GTPase concentrationlevels. In what follows, we present the model assumptions and derivation in §Section 2,the local stability analysis and numerical simulations in §Section 3, and discuss our find-ings in the context of GTPase signalling networks in §Section 4.2. M
ODEL D EVELOPMENT
In this section, we introduce our mathematical model for the GTPase coupled circuit(Fig.1B). We begin with outlining the model assumptions in §Section 2.1 and describethe reactions and governing equations in detail in §Section 2.2.2.1.
Assumptions.
Our model describes the time evolution of the concentrations forthe different system components. Table 1 contains the set of reactions in our system.In this work, we only consider the toggling of GTPases that are mediated by GEFs andGAPs that activate and inactivate them, respectively. Moreover, and guided directly byexperimental findings [32], we explicitly include the recruitment/engagement of tGEFby mG ∗ (feedforward represented by arrow 1 in Fig.1B) and two feedback loops (arrows2 and 3 in Fig.1B) representing the activation of mGAP by tGEF ∗ and tG ∗ , respectively.To develop the model equations, we consider a well-mixed regime to investigate how thekinetic reactions affects the dynamics of this GTPase circuit. We also assume that theconcentrations of the species are in large enough amounts that deterministic kinetics hold[64, 65]. Finally, for mathematical tractability, all reactions in the system are governedby mass-action and, therefore, nonlinear kinetics such as Hill functions or Michaelis-Menten are not considered [66].2.2. Governing Equations.
We developed a system of ODEs that describe coupled tog-gling of two switches, i.e., cyclical activation and inactivation of monomeric and trimericGTPases within the network motif shown in Fig.1B and described in Fig.1C. In what fol-lows, the brackets represent concentrations, which are nonnegative real numbers. Thesystem equations are given by
L.M. STOLERMAN, P. GHOSH, AND P. RANGAMANI (C) F IGURE A network motif in which two species of GTPases are in-terlinked. (A) Uncoupled monomeric and trimeric GTPase switches arerepresented by mGTPase and tGTPase, respectively. The black star de-notes the active forms. Activation and inactivation are regulated by GEFsand GAPs, where the first letter (m or t) indicates the associated GTPase.(B) Our proposed mathematical model describes the interaction betweenthe two GTPase switches. Arrows 1, 2, and 3 show the feedforward andfeedback loops that were found experimentally. (C) Description and bi-ological meaning of each arrow connecting the GTPase switches. Refer-ences: [32] for arrow 1 (*), [32] for arrow 2 (**), [32, 61] for arrow 3(***), and [61, 30, 62, 63, 31] for evidence of cooperativity between mand tGTPases (****).
TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 7 T ABLE
1. GTPase circuit reactions and rates used in the model
List of Reactions Reaction Rate mG ∗ activation mG + mGEF ∗ k mGon −−→ mG ∗ k mGon [ mGEF ∗ ][ mG ] mG ∗ inactivation mG ∗ + mGAP ∗ k mGoff −−→ mG k mGoff [ mGAP ∗ ][ mG ∗ ] Feedforward from mG* to tGEF (arrow 1) mG ∗ + tGEF k Ion −−→ tGEF ∗ k Ion [ tGEF ][ mG ∗ ] tG ∗ activation tG + tGEF ∗ k tGon −−→ tG ∗ k tGon [ tGEF ∗ ][ tG ] tG ∗ inactivation tG ∗ + tGAP ∗ k tGoff −−→ tG k tGoff [ tGAP ∗ ][ tG ∗ ] Feedback loop from tGEF to mGAP (arrow 2) tGEF ∗ + mGAP k IIon −−→ mGAP ∗ k IIon [ tGEF ∗ ][ mGAP ] Feedback loop tG* to mGAP (arrow 3) tG ∗ + mGAP k IIIon −−→ mGAP ∗ k IIIon [ mGAP ][ tG ∗ ] d [ mG ] dt = − k mGon [ mGEF ∗ ][ mG ] + k mGoff [ mGAP ∗ ][ mG ∗ ] (2.1) d [ mG ∗ ] dt = k mGon [ mGEF ∗ ][ mG ] − k mGoff [ mGAP ∗ ][ mG ∗ ] − k Ion [ tGEF ][ mG ∗ ] (2.2) d [ tG ] dt = − k tGon [ tGEF ∗ ][ tG ] + k tGoff [ tGAP ∗ ][ tG ∗ ] (2.3) d [ tG ∗ ] dt = k tGon [ tGEF ∗ ][ tG ] − k tGoff [ tGAP ∗ ][ tG ∗ ] − k IIIon [ mGAP ][ tG ∗ ] (2.4) d [ tGEF ] dt = − k Ion [ tGEF ][ mG ∗ ] (2.5) d [ tGEF ∗ ] dt = k Ion [ tGEF ][ mG ∗ ] − k IIon [ tGEF ∗ ][ mGAP ] (2.6) d [ mGAP ] dt = − k IIIon [ mGAP ][ tG ∗ ] − k IIon [ tGEF ∗ ][ mGAP ] (2.7) d [ mGAP ∗ ] dt = k IIIon [ mGAP ][ tG ∗ ] + k IIon [ mGAP ][ tGEF ∗ ] (2.8) where k mGon and k tGon are the activation rates of mG and tG, respectively. The inactivationrates for mG and tG are given by k mGoff and k tGoff . We denote k Ion as the rate of tGEFactivation by mG ∗ through the feedforward connection. Finally, the mGAP activationrates through the tGEF ∗ and tG ∗ feedback loops are given by k IIon and k IIIon , respectively.To complete the system definition, all model components must have nonnegative initialconditions. In particular, if k Ion = k IIon = k IIIon = 0 , then our system describes twouncoupled GTPase switches (Fig.1A) that each have the same dynamics of the singleGTPase model proposed in [4].
L.M. STOLERMAN, P. GHOSH, AND P. RANGAMANI
3. M
ATHEMATICAL ANALYSIS AND RESULTS
In this section, we explore the role of the feedforward and feedback loops on the systemdynamics. To simplify our mathematical analysis and investigate the isolated contribu-tions of the different network architectures, we assume that all activation rates, whennon-zero, are equal and given by k on , as well the inactivation rates for both m- and t-GTPase switches ( k mGoff = k tGoff = k off ). We also assume that the concentrations of [ mGEF ∗ ] and [ tGAP ∗ ] are constant in our model.It is convenient to rewrite our ODE system in the form d x dt = S. v(x) , where x represents the vector of concentrations for the different components, S is the sto-ichiometric matrix and v(x) is a vector with the different reaction rates [67, 68]. Thus wedefine the components x (1) = [ mG ] , x (2) = [ mG ∗ ] , x (3) = [ tG ] , x (4) = [ tG ∗ ] , x (5) =[ tGEF ] , x (6) = [ tGEF ∗ ] , x (7) = [ mGAP ] , and x (8) = [ mGAP ∗ ] . We also write thereaction velocities as v = k on [ mGEF ∗ ] x (1) , v = k off x (2) x (8) , v = k on x (2) x (5) , v = k on x (3) x (6) v = k off [ tGAP ∗ ] x (4) , v = k on x (6) x (7) , v = k on x (4) x (7) . The × stoichiometric matrix for the system given by Eqs. 2.1 – 2.8 is then given by S = arrow 1 arrow 2 arrow 3 mG − mG* − − tG − tG* − − tGEF − tGEF* − mGAP − − mGAP* (3.1)where the rows and columns of S (Eq. 3.1) represent the components and reactions,respectively. The right null space of S comprises the steady state flux solutions, and theleft null space contains the conservation laws of the system [67]. On the other hand,the column space contains the dynamics of the time-derivatives, and the rank of S isthe actual dimension of the system in which the dynamics take place. In the followingsubsections, we analyze Eqs. 2.1 – 2.8 when the mGTPase and tGTPase switches arecoupled through: (i) A feedforward connection mG ∗ → tGEF only (arrow 1), (ii)feedforward connection mG ∗ → tGEF and feedback loop tGEF → mGAP (arrows 1and 2) and (iii) feedforward connection mG ∗ → tGEF and feedback loops tGEF → mGAP and tG ∗ → mGAP (arrows 1, 2, and 3). TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 9
Feedforward connection: Recruitment of tGEF by active mGTPases ( mG ∗ → tGEF ) . To analyze Eqs. 2.1 – 2.8 with the feedforward connection only (arrow 1 inFig.1B), we assume k IIon = k IIIon = 0 , which means that the feedback loop arrows 2 and3 are not considered this first analysis. This represents the simple feedforward connec-tion of the two GTPase switches, which, in cells appears to be mediated via activation–dependent coupling of mG* to tGEF [32]. In this case, the stoichiometric matrix (Eq.3.1) is × . Conservation laws.
For this particular system, the total concentrations [ tG tot ] := [ tG ] +[ tG ∗ ] and [ tGEF tot ] := [ tGEF ] + [ tGEF ∗ ] are constant over time and are strictlypositive. For this reason, it is convenient to introduce the fractions T := [ tG ][ tG tot ] and G := [ tGEF ][ tGEF tot ] of inactive tGTPase and tGEF in the system, respectively, and let T ∗ and G ∗ denote the fraction of their active forms. We then use T + T ∗ = 1 and G + G ∗ = 1 torewrite the system in the form d [ mG ] dt = − k on [ mGEF ∗ ][ mG ] + k off [ mGAP ∗ ][ mG ∗ ] (3.2) d [ mG ∗ ] dt = k on [ mGEF ∗ ][ mG ] − k off [ mGAP ∗ ][ mG ∗ ] − k on [ tGEF tot ](1 − G ∗ )[ mG ∗ ] (3.3) d T ∗ dt = k on [ tGEF tot ] G ∗ (1 − T ∗ ) − k off [ tGAP ∗ ] T ∗ (3.4) d G ∗ dt = k on (1 − G ∗ )[ mG ∗ ] (3.5) From the stoichiometric matrix (Eq. 3.1), we observe that [ mG ] + [ mG ∗ ] + [ tGEF tot ] G ∗ = C (3.6)where C > is constant over time. We compute the left null space of the stoichiometricmatrix and confirm the total of three conservation laws in this case. The conservationlaw given by Eq. 3.6 reduces the system to three unknowns, which eases the steady stateand stability analysis. Steady states.
To find biologically plausible (nonnegative) steady states of the systemgiven by Eqs. 3.2 – 3.6, we must find (cid:91) [ mG ] , (cid:92) [ mG ∗ ] , (cid:99) G ∗ , and (cid:99) T ∗ such that the time-derivatives in Eqs. 3.2–3.5 are zero and the conservation law given by Eq. 3.6 is satisfied.Defining k = k off k on , we must solve the following system: [ mGEF ∗ ] (cid:91) [ mG ] − k [ mGAP ∗ ] (cid:92) [ mG ∗ ] − [ tGEF tot ](1 − (cid:99) G ∗ ) (cid:92) [ mG ∗ ] = 0[ tGEF tot ] (cid:99) G ∗ (1 − (cid:99) T ∗ ) − k [ tGAP ∗ ] (cid:99) T ∗ = 0(1 − (cid:99) G ∗ ) (cid:92) [ mG ∗ ] = 0 (cid:91) [ mG ] + (cid:92) [ mG ∗ ] + [ tGEF tot ] (cid:99) G ∗ = C From the third equation above, we must have (cid:92) [ mG ∗ ] = 0 or (cid:99) G ∗ = 1 . Thus we dividethe steady state analysis in these two cases and summarize our results in the followingproposition, whose proof can be found in the appendix A. Proposition 3.1.
Let k = k off k on . The steady states (cid:98) x = (cid:16) (cid:91) [ mG ] , (cid:92) [ mG ∗ ] , (cid:99) T ∗ , (cid:99) G ∗ (cid:17) of thesystem given by Eqs. 3.2 – 3.6 are given by • Steady state 1: (cid:98) x = (cid:32) , ,
11 + k [ tGAP ∗ ] C , C [ tGEF tot ] (cid:33) (3.7) if and only if C ≤ [ tGEF tot ] and • Steady state 2: (cid:98) x = (cid:18) k [ mGAP ∗ ][ mGEF ∗ ] + k [ mGAP ∗ ] ( C − [ tGEF tot ]) , [ mGEF ∗ ][ mGEF ∗ ] + k [ mGAP ∗ ] ( C − [ tGEF tot ]) , [ tGEF tot ][ tGEF tot ] + k [ tGAP ∗ ] , (cid:19) . (3.8) if and only if C ≥ [ tGEF tot ] . Given the explicit expressions for the steady states and the parameter range in whichthey exist, we perform a local stability analysis to determine if these states are stable orunstable under small perturbations. We adopt the classical linearization procedure basedon the powerful Hartman-Grobman theorem [59, 60]. We show that the steady statesare locally asymptotically stable , which means that any trajectory will be attracted to thesteady state provided the initial condition is sufficiently close.
Local Stability Analysis.
Using that [ mG ] = C − [ mG ∗ ] + [ tGEF tot ] G ∗ (from Eq. 3.6),we obtain the following three-dimensional system: d [ mG ∗ ] dt = f ([ mG ∗ ] , G ∗ , T ∗ ) d T ∗ dt = f ([ mG ∗ ] , G ∗ , T ∗ ) d G ∗ dt = f ([ mG ∗ ] , G ∗ , T ∗ ) where f ([ mG ∗ ] , G ∗ , T ∗ ) = k on [ mGEF ∗ ] ( C − [ mG ∗ ] − [ tGEF tot ]) − k off [ mGAP ∗ ][ mG ∗ ] − k on [ tGEF tot ](1 − G ∗ )[ mG ∗ ] , TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 11 f ([ mG ∗ ] , G ∗ , T ∗ ) = k on [ tGEF tot ] G ∗ (1 − T ∗ ) − k off [ tGAP ∗ ] T ∗ , and f ([ mG ∗ ] , G ∗ , T ∗ ) = k on (1 − G ∗ )[ mG ∗ ] . To perform the local stability analysis, we calculate the Jacobian matrix evaluated at thesteady state J (cid:104) (cid:92) [ mG ∗ ] , (cid:99) T ∗ , (cid:99) G ∗ (cid:105) = ∂f ∂ [ mG ∗ ] ∂f ∂ T ∗ ∂f ∂ G ∗ ∂f ∂ [ mG ∗ ] ∂f ∂ T ∗ ∂f ∂ G ∗ ∂f ∂ [ mG ∗ ] ∂f ∂ T ∗ ∂f ∂ G ∗ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:16) (cid:92) [ mG ∗ ] , (cid:99) T ∗ , (cid:99) G ∗ (cid:17) (3.9)and by showing that all its eigenvalues have a negative real part, we can prove that thesteady state is LAS [59], provided we further assume that the strict inequalities fromProposition 3.1 hold. This is the content of the following theorem. Theorem 3.1.
Let C be the conservation quantity from Eq. 3.6. Then,(1) If C < [ tGEF tot ] , the steady state 1 (Eq. 3.7) is LAS.(2) If C > [ tGEF tot ] , the steady state 2 (Eq. 3.8) is LAS.Proof. To simplify our notation, we introduce k mGEF = k on [ mGEF ∗ ] , k mGAP = k off [ mGAP ] ,and k tGAP = k off [ tGAP ∗ ] . All calculations were done with MATLAB’s R2019b sym-bolic toolbox and we proceed with the analysis of each case separately.(1) Suppose C < [ tGEF tot ] . As we have seen in the previous section, in this casethe steady state is given by Eq. 3.7. The Jacobian matrix (Eq. 3.9) is given by k on ( C − [ tGEF tot ]) − k mGEF − k mGAP − k mGEF [ tGEF tot ]0 − Ck on − k tGAP k tGAP k on [ tGEF tot ] Ck on + k tGAP − k on (cid:16) C [ tGEF tot ] − (cid:17) . The first eigenvalue in this case is given by λ = − Ck on − k tGAP and the othertwo ( λ and λ ) are such that λ + λ = k on ( C − [ tGEF tot ]) − k mGAP − k mGEF < and λ λ = − k on k mGEF ( C − [ tGEF tot ]) > from which we conclude that λ and λ are both negative and therefore the steadystate is LAS. (2) Suppose now that C > [ tGEF tot ] . Following our previous analysis, the steadystate is given by Eq. 3.8. The Jacobian matrix in this case is given by − k mGAP − k mGEF k mGEF [ tGEF tot ] (cid:16) k on ( C − tGEF tot ) k mGAP + k mGEF − (cid:17) − k tGAP − k on [ tGEF tot ] k tGAP k on [ tGEF tot ] k tGAP + k on [ tGEF tot ] − k on k mGEF ( C − [ tGEF tot ]) k mGAP + k mGEF and the eigenvalues are given by λ = − k tGAP − k on [ tGEF tot ] , λ = − k mGEF − k mGAP and λ = − k on k mGEF ( C − [ tGEF tot ]) k mGAP + k mGEF , which are all negative and this com-pletes the proof. (cid:3) Biological interpretation of the stability features of the feedforward connection mG ∗ → tGEF . The feedforward connection from mG* to tGEF allows for the emergence of thesteady state 1 (Eq. 3.7) with zero mG and mG* values. This zero concentration can beinterpreted as a scenario in which nearly all the available mG proteins are activated tomG*, and that nearly all the mG* species have successfully engaged with the availabletGEFs, thereby maximally recruiting tGEF on the Golgi membranes. The inequality
C < [ tGEF tot ] must hold for existence and local asymptotic stability to steady state1. Recalling the definition of G ∗ and that Eq. 3.6 holds for all times, including t = 0 ,this relationship between C and [ tGEF tot ] can be rewritten as [ mG ](0) + [ mG ∗ ](0) < [ tGEF ](0) where [ tGEF ](0) = [ tGEF tot ] − [ tGEF ∗ ](0) is initial concentration ofcytosolic tGEF that is yet to be recruited by mG* to the membranes. Thus, if the totalamount of mG protein is initially less than the concentration of tGEF in cells, Theorem3.1 ensures that the steady state 1 with no mG and mG* will emerge. Moreover, thereduced system (Eqs. 3.2–3.6) will converge to the steady state 1, provided the initialand steady state concentration values are sufficiently close. Similarly, the steady state 2will exist when [ mG ](0) + [ mG ∗ ](0) > [ tGEF ](0) . This can be interpreted as a scenariowhere the total amount of mG protein concentration is higher than the total concentrationof tGEF in cells. In this case, the reduced system will converge to steady state 2 wheresome distribution of mG, mG*, tG, tG* are present (Eq. 3.8), given sufficiently closeinitial and steady state concentration values.Fig.2 illustrates the two possible steady states (gray-colored “1” in the × table) pro-moted by the feedforward connection. steady state 1 can be interpreted as a configurationwhere the copy numbers of both active and inactive mGTPase are low, while the tGTPasecopy numbers remain high. On the other hand, in steady state 2 both m- and tGTPaseshave high copy numbers in both their active and inactive forms. Our results suggestthat the feedforward from mG to tGEF, which initiates the coupling between the two G TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 13 protein switches, can drive the system to two possible configurations depending on thecellular concentrations of total mG and tGEF. If the initial tGEF is larger than the totalmG, the feedforward connection will result in a significant decrease of the total mG andresult in the activation of a fraction of the tGEF ( (cid:99) G ∗ = C [ tGEF tot ] in Eq. 3.7). On the con-trary, if the initial tGEF is less than the total mG, then the available tGEF will be fullyengaged ( (cid:99) G ∗ = 1 in Eq. 3.8) and there will be a residual mG concentration in the system.We conclude that the initial difference between the copy numbers of total mG and tGEF(a cytosolic protein that is recruited to the membrane by mG*) is the main factor that willdetermine the steady state of the coupled GTPase switches.3.2. Feedforward connection mG ∗ → tGEF with feedback loop tGEF → mGAP :Recruitment of tGEF by active mGTPases and tGEF colocalization with mGAP.. We analyze the case where the feedback loop tGEF → mGAP (arrow 2 in Fig.1B)is added to the system with the feedforward connection. In cells, this feedback looprepresents a tGEF* colocalization with mGAP on Golgi membranes that facilitates therecruitment of GAP proteins [32]. In this case, the model equations are given by thefollowing system: d [ mG ] dt = − k on [ mGEF ∗ ][ mG ] + k off [ mGAP ∗ ][ mG ∗ ] (3.10) d [ mG ∗ ] dt = k on [ mGEF ∗ ][ mG ] − k off [ mGAP ∗ ][ mG ∗ ] − k on [ tGEF ][ mG ∗ ] (3.11) d [ tG ] dt = − k on [ tGEF ∗ ][ tG ] + k off [ tGAP ∗ ][ tG ∗ ] (3.12) d [ tG ∗ ] dt = k on [ tGEF ∗ ][ tG ] − k off [ tGAP ∗ ][ tG ∗ ] (3.13) d [ tGEF ] dt = − k on [ tGEF ][ mG ∗ ] (3.14) d [ tGEF ∗ ] dt = k on [ tGEF ][ mG ∗ ] − k on [ tGEF ∗ ][ mGAP ] (3.15) d [ mGAP ] dt = − k on [ mGAP ][ tGEF ∗ ] (3.16) d [ mGAP ∗ ] dt = k on [ mGAP ][ tGEF ∗ ] (3.17) where we keep all activation and inactivation rates at the same value ( k on and k off re-spectively), as done in Section 3.1. As we did in the previous section, we first analyzethe conservation laws of this particular system. In this case, the stoichiometric matrix(Eq. 3.1) is × . Conservation Laws.
We begin by observing that the total amount of tGTPase is con-served in this system. Thus we may use the fraction T ∗ as in Section 3.1 and that is thefirst conservarion law. The total amount of mGAP is also conserved, as we sum Eqs. [ mGAP ] = [ mGAP tot ] − [ mGAP ∗ ] (3.18)and substitute the above expression for [ mGAP ] in Eqs. 3.15 and 3.17. We choose tokeep the concentrations of mGAP as a variable for notational simplicity and do not defineits fraction. Summing Eqs. 3.10, 3.11, 3.15, and 3.17, and integrating over time, we get [ mG ] + [ mG ∗ ] + [ tGEF ∗ ] + [ mGAP ∗ ] = C (3.19)where C ≥ is constant over time. Moreover, Eqs. 3.14, 3.15, and 3.17 when summedand integrated give [ tGEF ] + [ tGEF ∗ ] + [ mGAP ∗ ] = C (3.20)for C ≥ also constant. We compute the left null space of the stoichiometric matrix(Eq. 3.1) and confirm the total of four conservation laws in this case. The reduced systemis given by the following equations: d [ mG ] dt = − k on [ mGEF ∗ ][ mG ] + k off [ mGAP ∗ ][ mG ∗ ] (3.21) d [ mG ∗ ] dt = k on [ mGEF ∗ ][ mG ] − k off [ mGAP ∗ ][ mG ∗ ] − k on [ tGEF ][ mG ∗ ] (3.22) d T ∗ dt = k on [ tGEF ∗ ](1 − T ∗ ) − k off [ tGAP ∗ ] T ∗ (3.23) d [ tGEF ] dt = − k on [ tGEF ][ mG ∗ ] (3.24) d [ tGEF ∗ ] dt = k on [ tGEF ][ mG ∗ ] − k on [ tGEF ∗ ] ([ mGAP tot ] − [ mGAP ∗ ]) (3.25) d [ mGAP ∗ ] dt = k on ([ mGAP tot ] − [ mGAP ∗ ]) [ tGEF ∗ ] (3.26) with the conservation laws given by Eqs. 3.19 and 3.20. In what follows, we calculatethe steady states of the system. Steady states and local stability analysis.
We begin by introducing k = k off k on to sim-plify our notation. To find the steady states, we must find nonnegative solutions of thefollowing system: − [ mGEF ∗ ] (cid:91) [ mG ] + k (cid:92) [ mGAP ∗ ] (cid:92) [ mG ∗ ] = 0 (3.27) (cid:92) [ tGEF ∗ ](1 − (cid:99) T ∗ ) − k [ tGAP ∗ ] (cid:99) T ∗ = 0 (3.28) (cid:92) [ tGEF ] (cid:92) [ mG ∗ ] = 0 (3.29) (cid:16) [ mGAP tot ] − (cid:92) [ mGAP ∗ ] (cid:17) (cid:92) [ tGEF ∗ ] = 0 (3.30) (cid:91) [ mG ] + (cid:92) [ mG ∗ ] + (cid:92) [ tGEF ∗ ] + (cid:92) [ mGAP ∗ ] = C (3.31) (cid:92) [ tGEF ] + (cid:92) [ tGEF ∗ ] + (cid:92) [ mGAP ∗ ] = C (3.32) TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 15
From Eq. 3.29, we must have (cid:92) [ tGEF ] = 0 or (cid:92) [ mG ∗ ] = 0 . Moreover, from Eq. 3.30, (cid:92) [ tGEF ∗ ] = 0 or (cid:92) [ mGAP ∗ ] = [ mGAP tot ] and thus we have four possible combinationsto analyze.We study each case separately and obtain the necessary and sufficient inequalities involv-ing the parameters C , C , and [ mGAP tot ] that ensure the existence of each steady state.As we did in the previous section, we also show that the steady states are LAS providedthe strict inequalities are satisfied. We summarize our analysis in the following theorem,whose proof can be found in Appendix B. Theorem 3.2.
The steady states (cid:98) x = (cid:16) (cid:91) [ mG ] , (cid:92) [ mG ∗ ] , (cid:99) T ∗ , (cid:92) [ tGEF ] , (cid:92) [ tGEF ∗ ] , (cid:92) [ mGAP ∗ ] (cid:17) of the system given by Eqs. 3.19 - 3.26 are given by • Steady state 1: (cid:98) x = (cid:18) , , C − [ mGAP tot ]( C − [ mGAP tot ]) + k [ tGAP ∗ ] , C − C , C − [ mGAP tot ] , [ mGAP tot ] (cid:19) (3.33) if and only if C ≥ C and C ≥ [ mGAP tot ] . The steady state is LAS if C > C and C > [ mGAP tot ] . • Steady state 2: (cid:98) x = (cid:18) k [ mGAP tot ] ( C − C )[ mGEF ∗ ] + k [ mGAP tot ] , [ mGEF ∗ ] ( C − C )[ mGEF ∗ ] + k [ mGAP tot ] ,C − [ mGAP tot ]( C − [ mGAP tot ]) + k [ tGAP ∗ ] , , C − [ mGAP tot ] , [ mGAP tot ] (cid:19) (3.34) if and only if C ≥ C and C ≥ [ mGAP tot ] . The steady state is LAS if C > C and C > [ mGAP tot ] . • Steady state 3: (cid:98) x = (0 , , , C − C , , C ) (3.35) if and only if C ≥ C and C ≤ [ mGAP tot ] . The steady state is LAS if C > C and C < [ mGAP tot ] . • Steady state 4: (cid:98) x = (cid:18) kC [ mGEF ∗ ] + kC ( C − C ) , [ mGEF ∗ ][ mGEF ∗ ] + kC ( C − C ) , , , , C (cid:19) (3.36) if and only if C ≥ C and C ≤ [ mGAP tot ] . The steady state is LAS if C > C and C < [ mGAP tot ] .Biological interpretation of the stability features of the feedforward connection mG ∗ → tGEF with single feedback loop tGEF → mGAP . The feedforward connection mG ∗ → tGEF together with the feedback loop tGEF → mGAP (arrows 1 and 2 in Fig.1B, re-spectively) allows for the emergence of four steady states. steady states 1 and 2 (Eqs.3.33 and 3.34) are similar to the two steady states obtained in Section 3.1, although with different concentration values. On the other hand, steady states 3 and 4 (Eqs. 3.35 andEqs. 3.36) newly emerge in the system, both with tG* attaining zero concentration. Thiszero concentration can be interpreted as a scenario in which nearly all the available tGT-Pase is inactivated. Recalling the definitions of C and C and the fact that Eqs. 3.19 and3.20 hold at all times, including at t = 0 , we can write C = [ mG ](0) + [ mG ∗ ](0) +[ tGEF ∗ ](0) + [ mGAP ∗ ](0) and C = [ tGEF ](0) + [ tGEF ∗ ](0) + [ mGAP ∗ ](0) . Inthis way, from the inequalities obtained in Theorem 3.2 for C and C , we obtain rela-tionships among the initial conditions of the original system (Eqs. 3.10 – 3.17) that areassociated with each one of the four steady states.For the existence and local asymptotic stability of steady state 1 (Eq. 3.33), where mGand mG* have zero concentration values, the inequalities C > C and C > [ mGAP tot ] must hold. The first inequality can be written as [ mG ](0) + [ mG ∗ ](0) < [ tGEF ](0) ,which was obtained in Section 3.1 as the existence condition for the steady state withno mG and mG* (Eq. 3.7). On the other hand, the inequality C > [ mGAP tot ] can bewritten as [ mG ](0) + [ mG ∗ ](0) + [ tGEF ∗ ](0) > [ mGAP ](0) , where [ mGAP ](0) isthe initial concentration of cytosolic mGAP that is yet to be recruited by tGEF* to themembranes. Therefore, two conditions guarantee the existence of steady state 1: (1) Thetotal amount of mG protein must be initially less than the concentration of tGEF and (2)The sum of the concentrations of total mG protein and tGEF* must be initially higherthan the concentration of mGAP. If both conditions hold, then Theorem 3.2 ensures thatsteady state 1 will emerge and the reduced system (Eqs. 3.21 – 3.26 along with Eqs.3.19 and 3.20) will converge to the steady state 1, provided the initial and steady stateconcentration values are sufficiently close.A similar analysis holds for steady states 2, 3 and 4. For simplicity, we present the re-quired initial conditions for each steady state without repeating the conclusions that fol-lows from Theorem 3.2. For steady state 2 ((Eq. 3.34)), where mG, mG*, tG, and tG* arepresent, the inequalities C > C and C > [ mGAP tot ] become [ mG ](0) + [ mG ∗ ](0) > [ tGEF ](0) and [ tGEF ](0) + [ tGEF ∗ ](0) > [ mGAP ](0) , respectively. Hence, the to-tal amount of mG protein must be initially higher than the concentration of tGEF andthe total amount of tGEF must be initially higher than concentration of mGAP. Forsteady state 3 (Eq. 3.35), where mG and mG* have zero concentration values andthe tGTPase is fully inactivated, the inequalities C > C and C < [ mGAP tot ] be-come [ mG ](0) + [ mG ∗ ](0) < [ tGEF ](0) and [ mG ](0) + [ mG ∗ ](0) + [ tGEF ∗ ](0) < [ mGAP ](0) , respectively. Hence, the total amount of mG protein must be initially lessthan the concentration of tGEF and the sum of the concentrations of total mG proteinand tGEF* must be initially less than the concentration of mGAP. For steady state 4(Eq. 3.55), where mG and mG* are present and tG* concentration is zero, the inequal-ities C > C and C < [ mGAP tot ] become [ mG ](0) + [ mG ∗ ](0) > [ tGEF ](0) and [ tGEF ](0) + [ tGEF ∗ ](0) < [ mGAP ](0) , respectively. Hence, the total amount of mG TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 17 protein must be initially higher than the concentration of tGEF and the total amount oftGEF must be initially less than the concentration of mGAP.Fig.2 illustrates the four possible steady states (gray-colored “1+2” in the × ta-ble) promoted by the feedforward connection mG ∗ → tGEF and the feedback loop tGEF → mGAP . Steady states 1 and 2 have the same interpretation of the two steadystates obtained in Section 3.1. On the other hand, steady states 3 and 4 were obtainedthrough the sole contribution of the feedback loop tGEF → mGAP . These states sharethe common feature of having tGTPase fully inactivated. However, steady state 3 can beinterpreted as a configuration where the copy numbers of mG and mG* are low, while insteady state 4, these copy numbers are high.3.3. Feedforward connection mG ∗ → tGEF with feedback loops tGEF → mGAP and tG ∗ → mGAP : Recruitment of tGEF by active mGTPases, tGEF colocaliza-tion with mGAP, and activation of mGAP by active tGTPases. We analyze the casewhere the feedback loop tG ∗ → mGAP (arrow 3 in Fig.1B) is added to the system withthe feedforward connection and feedback loop tGEF → mGAP . This connection rep-resents the release of free Gβγ promoting mGAP activation. We analyze the full systemgiven by Eqs. 2.1 – 2.8 in the case where where all rates are equal. The model equationsare thus given by the following system: d [ mG ] dt = − k on [ mGEF ∗ ][ mG ] + k off [ mGAP ∗ ][ mG ∗ ] (3.37) d [ mG ∗ ] dt = k on [ mGEF ∗ ][ mG ] − k off [ mGAP ∗ ][ mG ∗ ] − k on [ tGEF ][ mG ∗ ] (3.38) d [ tG ] dt = − k on [ tGEF ∗ ][ tG ] + k off [ tGAP ∗ ][ tG ∗ ] (3.39) d [ tG ∗ ] dt = k on [ tGEF ∗ ][ tG ] − k off [ tGAP ∗ ][ tG ∗ ] − k on [ mGAP ][ tG ∗ ] (3.40) d [ tGEF ] dt = − k on [ tGEF ][ mG ∗ ] (3.41) d [ tGEF ∗ ] dt = k on [ tGEF ][ mG ∗ ] − k on [ mGAP ][ tGEF ∗ ] (3.42) d [ mGAP ] dt = − k on [ mGAP ][ tG ∗ ] − k on [ tGEF ∗ ][ mGAP ] (3.43) d [ mGAP ∗ ] dt = k on [ mGAP ][ tG ∗ ] + k on [ tGEF ∗ ][ mGAP ] . (3.44) In what follows, we compute the conservation laws and four 1-parameter steady statefamilies for the system given by Eqs. 3.37 – 3.44. We also obtain the necessary condi-tions for the conserved quantities that guarantee the existence of each steady state family.
Conservation Laws.
As in 3.2, we also observe that the total amount of mGAP is con-stant over time, so Eq. 3.18 still holds. On the other hand, the total tGTPase follows a new conservation law that we derive here. Summming Eqs. 3.37 – 3.40, 3.42, and 3.44,we have [ mG ] + [ mG ∗ ] + [ tG ] + [ tG ∗ ] + [ tGEF ∗ ] + [ mGAP ∗ ] = ˜ C . (3.45)and summing Eqs. 3.39 – 3.42 and Eq. 3.44 and integrating over time, we obtain [ tG ] + [ tG ∗ ] + [ tGEF ] + [ tGEF ∗ ] + [ mGAP ∗ ] = ˜ C (3.46)where ˜ C and ˜ C must be nonnegative constants. We compute the left null space of thestoichiometric matrix (Eq. 3.1) and confirm the total of three conservation laws, whichare given by Eqs. 3.18, 3.45, and 3.46. These equations reduce Eqs. 3.37 – 3.44 to afive-dimensional system, whose steady states can be obtained. Steady states.
We compute the steady states of the system when the time derivatives inEqs. 3.37 – 3.44 are equal to zero. Denoting k = k off k on as in the previous section and re-moving the linearly dependent equations, the problem reduces to finding the nonnegativesolutions of the following system: − [ mGEF ∗ ] (cid:91) [ mG ] + k (cid:92) [ mGAP ∗ ] (cid:92) [ mG ∗ ] = 0 (3.47) − (cid:92) [ tGEF ∗ ] (cid:100) [ tG ] + k [ tGAP ∗ ] (cid:91) [ tG ∗ ] = 0 (3.48) (cid:16) [ mGAP tot ] − (cid:92) [ mGAP ∗ ] (cid:17) (cid:91) [ tG ∗ ] = 0 (3.49) (cid:92) [ tGEF ] (cid:92) [ mG ∗ ] = 0 (3.50) (cid:92) [ tGEF ∗ ] (cid:16) [ mGAP tot ] − (cid:92) [ mGAP ∗ ] (cid:17) = 0 (3.51) along with the conservation laws given by Eqs. 3.18, 3.45, and 3.46.Eq. 3.48 gives tG ∗ = (cid:92) [ tGEF ∗ ] (cid:100) [ tG ] k [ tGAP ∗ ] and Eq. 3.49 then becomes (cid:16) [ mGAP tot ] − (cid:92) [ mGAP ∗ ] (cid:17) (cid:92) [ tGEF ∗ ] (cid:100) [ tG ] = 0 . From Eq. 3.51, we conclude that (cid:100) [ tG ] can be any nonnegative real number satisfying Eqs.3.45 and 3.46. We define ξ = (cid:100) [ tG ] and characterize four ξ -dependent families of steadystates similarly as we did in Section 3.2. We summarize our results in the followingtheorem, whose proof can be found in the appendix C. Theorem 3.3.
The ξ -dependent families of steady states (cid:98) x ξ = (cid:16) (cid:91) [ mG ] , (cid:92) [ mG ∗ ] , (cid:100) [ tG ] , (cid:91) [ tG ∗ ] , (cid:92) [ tGEF ] , (cid:92) [ tGEF ∗ ] , (cid:92) [ mGAP ∗ ] (cid:17) of the system given by Eqs. 3.47 – 3.51 with the conservation laws given by Eqs. 3.18,3.45, and 3.46 are given by TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 19 • Family 1: (cid:98) x ξ = , , ξ, (cid:16) ˜ C − [ mGAP tot ] − ξ (cid:17) ξk [ tGAP ∗ ] + ξ , ˜ C − ˜ C , (cid:16) ˜ C − [ mGAP tot ] − ξ (cid:17) k [ tGAP ∗ ] k [ tGAP ∗ ] + ξ , [ mGAP tot ] (cid:33) (3.52) only if ≤ ξ + [ mGAP tot ] ≤ ˜ C ≤ ˜ C . • Family 2: (cid:98) x ξ = (cid:18) ( ˜ C − ˜ C ) k [ mGAP tot ][ mGEF ∗ ] + k [ mGAP tot ] , ( ˜ C − ˜ C )[ mGEF ∗ ][ mGEF ∗ ] + k [ mGAP tot ] ,ξ, (cid:16) ˜ C − [ mGAP tot ] − ξ ] (cid:17) ξk [ tGAP ∗ ] + ξ , , (cid:16) ˜ C − [ mGAP tot ] − ξ ] (cid:17) k [ tGAP ∗ ] k [ tGAP ∗ ] + ξ , [ mGAP tot ] (cid:33) (3.53) only if ≤ ξ + [ mGAP tot ] ≤ ˜ C ≤ ˜ C . • Family 3: (cid:98) x ξ = (cid:16) , , ξ, , ˜ C − ˜ C , , ˜ C − ξ (cid:17) (3.54) only if max(0 , ˜ C − [ mGAP tot ]) ≤ ξ ≤ ˜ C ≤ ˜ C . • Family 4: (cid:98) x ξ = (cid:32) ( ˜ C − ˜ C ) k ( ˜ C − ξ )[ mGEF ∗ ] + k ( ˜ C − ξ ) , ( ˜ C − ˜ C )[ mGEF ∗ ][ mGEF ∗ ] + k ( ˜ C − ξ ) , ξ, , , , ˜ C − ξ (cid:33) . (3.55) only if max(0 , ˜ C − [ mGAP tot ]) ≤ ξ ≤ ˜ C ≤ ˜ C .Biological interpretation of the stability features of the feedforward connection mG ∗ → tGEF with two feedback loops tGEF → mGAP and tG ∗ → mGAP . The feedbackloop tG ∗ → mGAP , when added to the feedforward connection mG ∗ → tGEF and thefeedback loop tGEF → mGAP , allows for the emergence of four steady state families.Those families have the inactive tGTPase with different range values that can be obtainedin the steady state analysis. Interestingly, the Families 1 – 4 resemble the steady states1 – 4 from Section 3.2. Family 1 has no mG and mG* at steady state, and both tG andtG* have non-zero steady state values (similarly to steady state 1). Moreover, Family 2has both m- and tGTPases with non-zero steady states (similarly to steady state 2). ForFamily 3, mG and mG* steady state values are zero, and the tGTPase is fully inactivated(similarly to steady state 3). Finally, Family 4 has tGTPase is fully inactivated, and bothmG and mG* have non-zero steady states (similarly to steady state 4). Recalling thedefinitions of ˜ C and ˜ C and the fact that Eqs. 3.45 and 3.46 hold at all times, including t = 0 , we can infer necessary relationships among the initial conditions for each steadystate family.The inequality ˜ C ≤ ˜ C can be rewritten as [ mG ](0) + [ mG ∗ ](0) ≤ [ tGEF ](0) is nec-essary for the emergence of Family 1 (Eq. 3.52) whith zero mG and mG* values, which can be interpreted as a scenario in which nearly all the available mG proteins are acti-vated to mG*, and that nearly all the mG* species have successfully engaged with theavailable tGEFs, thereby maximally recruiting tGEF on the membranes. For Family 1, [ mGAP tot ] ≤ ˜ C also holds and can be written as [ mGAP ](0) ≤ [ mG ](0) + [ mG ∗ ](0) +[ tG ](0) + [ tG ∗ ](0) + [ tGEF ∗ ](0) , where [ mGAP ](0) is the initial concentrations of cy-tosolic mGAP that is yet to be recruited by tGEF* and tG* to the membranes. Therefore,two initial conditions are necessary for the existence of Family 1: The total amount ofmG protein must be initially less than the concentration of tGEF and (2) The summedconcentrations of total mG, total tG and tGEF* must be initially higher than the con-centration of mGAP. Finally, the inequality ≤ ξ + [ mGAP tot ] ≤ ˜ C can be writtenas ≤ ξ ≤ [ mG ](0) + [ mG ∗ ](0) + [ tG ](0) + [ tG ∗ ](0) + [ tGEF ∗ ](0) − [ mGAP ](0) .Remarkably, we conclude that the initial balance between the summed concentrations oftotal mG, total tG, tGEF* and the available mGAP is the upper bound for the tG con-centration, which completely characterizes the necessary conditions for the emergenceof Family 1.A similar analysis can be done for Families 2, 3, and 4. For the existence of Family 2(Eq. 3.53), where mG, mG* tG, tG* are present (when ξ > ), the inequalities ˜ C ≤ ˜ C and [ mGAP tot ] ≤ ˜ C must hold and can be rewritten as [ mG ](0) + [ mG ∗ ](0) ≥ [ tGEF ](0) and [ mGAP ](0) ≤ [ tG ](0) + [ tG ∗ ](0) + [ tGEF ](0)[ tGEF ∗ ](0) . Hence, thetotal amount of mG protein must be initially higher than the concentration of tGEF andthe summed concentrations of total tG and total tGEF proteins must be initially higherthan the concentration of mGAP. Finally, the inequality ≤ ξ + [ mGAP tot ] ≤ ˜ C indicates that initial balance between the summed concentrations of total tG, total tGEFand the available mGAP is the upper bound for the tG concentrations. For Family 3(Eq. 3.54), where mG and mG* have zero concentration values and the tGTPase is fullyinactivated, the inequality ˜ C ≤ ˜ C becomes [ mG ](0) + [ mG ∗ ](0) ≤ [ tGEF ](0) . Asfor Family 1, the total amount of mG protein must be initially less than the concentrationof tGEF. Moreover, from ˜ C − [ mGAP tot ] ≤ ξ , the initial balance between the summedconcentrations of total mG, total tG, tGEF* and the available mGAP is the lower boundfor the tG concentration. For Family 4 (Eq. 3.55), where mG and mG* are present andtG* concentration is zero, ˜ C ≤ ˜ C becomes [ mG ](0) + [ mG ∗ ](0) ≥ [ tGEF ](0) . As forFamily, 2 the total amount of mG protein must be initially higher than the concentrationof tGEF. Moreover, from ˜ C − [ mGAP tot ] ≤ ξ , the initial balance between the summedconcentrations of total tG, total tGEF and the available mGAP is the lower bound for thetG concentrations.Fig.2 illustrates the four Families (gray-colored “1+2+3” in the × table) promoted bythe feedforward connection mG ∗ → tGEF and the feedback loops tGEF → mGAP and tG ∗ → mGAP . Families 1 and 2 have a similar interpretation of the steady states 1and 2 obtained in Section 3.1 and Section 3.2. On the other hand, Families 3 and 4 wereobtained through contributions of the feedback loop tG ∗ → mGAP . These states share TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 21 F IGURE Possible steady states promoted by the feedforward con-nection and the feedback loops . For the three combinations of arrows(1, 1+2, and 1+2+3) chosen in our study, we calculate the steady statesolutions for the coupled GTPase circuit model. We then characterizefour steady state configurations that emerge from our analysis, depend-ing on low/high copy numbers of the m- and tGTPases. The feedforwardconnection mG ∗ → tGEF (represented by “1”) allows for the emer-gence of two steady states with low/high mG and mG* concentrations,while tG* steady state concentration remained high in both cases. On theother hand, low tG* steady state concentrations were obtained when thefeedback loop tGEF → mGAP was added to the system (representedby “1+2”). Finally, the feedback loop tG ∗ → mGAP allowed for theemergence of four parametrized families of steady state within the samelow/high configurations. (*) Families 2 and 4 can have higher tG* steadystates provided the inactive tGTPase concentrations assume strictly posi-tive values.the common feature of having tGTPase fully inactivated. As for steady states 3 and 4,Family 3 can be interpreted as a configuration where the copy numbers of mG and mG*are low, while in Family 4, those copy numbers are high.3.4. Numerical Simulations.
To complete our mathematical analysis, we numericallyinvestigate the range of initial conditions in which the trajectories of the system convergeto the different steady states. In particular, we illustrate the so-called basins of attraction [69] of the steady states, considering the same combination of connections between thetwo GTPase switches from §3.1 – §3.3. In Fig.3, we explore the case where the two GT-Pase switches are coupled through the feedforward connection mG ∗ → tGEF ( Fig.3A).We color the trajectories of the system according to the comparison between the initial conditions [ mG ](0) + [ mG ∗ ](0) and [ tGEF ∗ ](0) from the steady state analysis in Sec-tion 3.1. For fixed [ mG ∗ ](0) and [ tGEF ](0) values, we consider [ mG ](0) ranging from 0to 10 µM and therefore [ mG ](0) + [ mG ∗ ](0) can be less of higher than [ tGEF ](0) (blueor red-colored lines and dots). For all simulations, we plot the trajectories of the systemuntil equilibrium is reached. In the G ∗ × [ mG tot ] plane (Fig.3B), we observe a linearrelationship between these two quantities, where the black arrows indicate the time di-rection. If [ mG ](0) + [ mG ∗ ](0) < [ tGEF ](0) , the system converges to a state where noactive mGTPase exists (blue colored trajectories in Figs. 3B and 3C). On the other hand,if [ mG ](0) + [ mG ∗ ](0) > [ tGEF ](0) , the system converges a state where the concen-tration of the active and inactive mGTPase are positive at the final time (red-colored tra-jectories). To visualize these results in terms of dose-response curves, in Fig.3D we plotthe final-state values of [ mG tot ] and G ∗ (denoted by s.s) as a function of [ mG ](0) . Thetrajectories in the T ∗ × [ mG tot ] plane are shown in Fig.3E. We observe a detail showingthat T ∗ reaches a fixed final value around 0.97 when [ mG ](0) + [ mG ∗ ](0) > [ tGEF ](0) (see magnified view). We observe that the trajectories converge to steady states thatagree with the local stability results from §Section 3.1. This suggests that the condi-tions [ mG ](0) + [ mG ∗ ](0) < [ tGEF ](0) and [ mG ](0) + [ mG ∗ ](0) > [ tGEF ](0) are notonly valid in a neighborhood of the steady states, but also hold for other initial valuessatisfying those inequalities.Fig.4 illustrates the dynamics of the system when the feedback loop tGEF → mGAP (Arrow 2) is added to the feedforward connection (Fig.4A). In Fig.4B, we plot several [ tGEF ∗ ] trajectories starting at [ tGEF ∗ ] = 5 µM for different [ mG ](0) and [ mGAP ](0) values. The resulting rich variety of curves indicate the sensitivity of the system tothese initial conditions. In Fig.4C, different dose-response curves are generated to showthe steady state tGEF* values. If [ mGAP ](0) = 0 µM (blue and red dots), only thefeedforward connection affects the system, since mGAP cannot be activated by tGEF*.When [ mGAP ](0) = 1 (green squares), a similar steady state profile emerges, with [ tGEF ∗ ] s.s increasing for [ mG ](0) ≤ µM and remaining constant [ mG ](0) > .When [ mGAP ](0) = 8 µM , [ tGEF ∗ ] is zero for [ mG ](0) < µM and increases un-til [ mG ](0) < µM . For [ mG ](0) > , the steady state achieves its maximum valueslightly above [ tGEF ∗ ] > . Finally when [ mGAP ](0) = 11 µM , tGEF ∗ becomesfully recruited by mGAP and the [ tGEF ∗ ] s.s is zero for all [ mG ](0) values. In Fig.4D,we scan the space of initial amounts of mG and mGAP. When [ mGAP ](0) > µM , the [ tGEF ∗ ] s.s is zero, while for [ mGAP ](0) < µM is becomes nonzero and dependentof [ mG ](0) . In Figs. 4 E, F, and G, we analyze the [ tG ∗ ] concentration values and obtainsimilar results.Fig.5 illustrates the dynamics of the system when the feedback loops tGEF → mGAP and tG ∗ → mGAP are added to the feedforward connection (Fig.5A). In Fig.4B, weplot several [ tGEF ∗ ] trajectories starting at [ tGEF ∗ ] = 5 µM for different [ mG ](0) and [ mGAP ](0) values. In Fig.4C, different dose-response curves are generated to show TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 23 (A) (B) (C)(D) (E) F IGURE Trajectories of the system and steady states (arrow 1 (A)Schematics with the coupled GTPase switches and a feedforward connection mG ∗ → tGEF , represented by arrow 1. (B) [ mG ](0) was changed from 0 to10 µM and the trajectories of the system were calculated until equilibrium wasreached. In the G ∗ × [ mG tot ] plane, a linear relationship emerges. The blackarrows indicate the direction of time. If [ mG ](0) > µM , the system convergesto a final-state where the concentrations of the active and inactive mGTPase arenonzero. On the other hand, when [ mG ](0) < µM , the trajectories convergea final-state with no mGTPase exists (blue colored lines). (C) Trajectories ofthe active ( [ mG ] ) vs inactive mGTPase ([mG]) for [ mG ](0) . (D) Dose responsecurves show the steady states (denoted by s.s) for the total mGTPase concentra-tion and fraction of active tGEF ( G ) depending on [ mG ](0) in the two differ-ent scenarios. (E) The dynamics in the T × [ mG tot ] plane. Parameter values: [ mG ∗ ](0) = 0 µM, T ∗ (0) = 0 . , G ∗ (0) = 0 . , [ mGAP ] = 1 µM, [ mGEF ] =1 µM, k on = 3( s.µM ) − , k off = 1( s.µM ) − , [ tGAP ] = 1 µM, [ tGEF tot ] =10 µM, [ tG tot ] = 10 µM . Simulation time: s for panels B and E, s for panelsC, D, F, and G. Numerical simulations were performed using the solver ode15sin Matlab R2018a. All parameters were arbitrarily chosen only to illustrate thedynamic features of the model. the steady state tGEF* values. As in the previous case with only one feedback loop,if [ mGAP ](0) = 0 µM (blue and red dots), mGAP cannot be activated by tGEF*.When [ mGAP ](0) = 1 (green squares), a similar steady state profile emerges, with [ tGEF ∗ ] s.s increasing for [ mG ](0) ≤ µM and remaining constant [ mG ](0) > .When [ mGAP ](0) = 8 and 11 µM , [ tGEF ∗ ] increases until [ mG ](0) < µM . For [ mG ](0) > , the steady state achieves its maximum value. In Fig.5D, we scan the spaceof initial amounts of mG and mGAP and we observe a more graded response in compar-ison with Fig.4. In Figs. 5 E, F, and G, we analyze the [ tG ∗ ] concentration values andobtain similar results. (A)(B) (C)(E) (F) (D)(G) F IGURE Trajectories of the system and steady states (s.s) (arrows 1 and 2). (A)Schematics of the coupled GTPases with the feedforward connection mG ∗ → tGEF (arrow 1) and the feedback loop tGEF → mGAP (arrow 2). (B) [ tGEF ∗ ] trajecto-ries for [ mGAP ](0) =
0, 1, 8, and 11 µM . For each [ mGAP ](0) value, we plot twocurves for [ mG ](0) = µM (solid) (C) Dose response curves show [ tGEF ∗ ] s.s when [ mG ](0) ranges from 0 to 10 µM . If [ mGAP ](0) = 0 µM (blue andred dots), there will be no mGAP activation and therefore no effects of the feedback. For [ mGAP ](0) > µM , the feedback becomes effective and generate different [ tGEF ∗ ] responses. (D) Colormap for [ tGEF ∗ ] s.s concentrations for a range of [ mG ](0) and [ mGAP ](0) values. A sharp decrease on [ tGEF ∗ ] occurs when [ mGAP ](0) ≥ µM .When [ mGAP ](0) < µM , the [ tGEF ∗ ] s.s depend on [ mG ](0) . (E) [ tG ∗ ] trajecto-ries for [ mGAP ](0) =
0, 5, 9, and 11 µM and same [ mG ](0) . (F) Dose response curvesfor [ tG ∗ ] s.s depend on [ mGAP ](0) . (G) Colormap for [ tG ∗ ] s.s.; lower tG* concentra-tions result from higher [ mGAP ](0) values, since tGEF* is recruited for mGAP activa-tion. Parameter values: k on = 3( s.µM s ) − , k o f f = 1( s.µM ) − , [ mG ∗ ](0) = 0 µM , [ tGEF tot ](0) = 10 µM , [ tGEF ∗ ](0) = 5 µM , T ∗ (0) = 0 . , [ tG tot ] = 10 µM , [ mGAP ∗ ](0) = 1 µM , [ tGAP ∗ ](0) = 1 µM , [ mGEF ∗ ] = 1 µM . Simulation times: s (B and E) and s (C, D, F, and G). Numerical simulations were performed usingthe solver ode23s in Matlab R2018a. All parameters were arbitrarily chosen only toillustrate the dynamic features of the model. In Fig.6, we investigate the space of initial conditions for mG* and mGAP* in whichthe system converges to the different steady states. Fig.6A shows the simplest systemwhere the two GTPase switches are connected by the feedforward mG ∗ → tGEF . Twosteady states are obtained depending on the initial amount of mG*. For [ mG ∗ ](0) < [ tGEF ](0) − [ mGAP ∗ ](0) = 5 µM , the trajectories converge to steady state 1 with TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 25 (A)(B) (C)(E) (D)(F) (G) F IGURE Trajectories of the system and steady states (s.s) (arrows 1, 2, and3). (A) Schematics of the coupled GTPases with the feedforward connection mG ∗ → tGEF (arrow 1) and the feedback loops tGEF → mGAP (arrow 2) and tG ∗ → mGAP (arrow 3). (B) [ tGEF ∗ ] trajectories for [ mGAP ](0) =
0, 1, 8, and 11 µM .For each [ mGAP ](0) value, we plot two curves for [ mG ](0) = µM . (C)Dose response curves show [ tGEF ∗ ] s.s when [ mG ](0) ranges from 0 to 10 µM . If [ mGAP ](0) = 0 µM (blue and red dots), there will be no mGAP activation and there-fore no effects of the feedback loops. For [ mGAP ](0) > µM , the feedback becomeseffective and generate different [ tGEF ∗ ] responses. (D) Colormap for [ tGEF ∗ ] s.sconcentrations for a range of [ mG ](0) and [ mGAP ](0) values. A more graded decreaseon [ tGEF ∗ ] occurs when [ mGAP ](0) ≥ µM in comparison with Fig.4D. (E) [ tG ∗ ] trajectories for [ mGAP ](0) =
0, 5, 9, and 11 µM and same [ mG ](0) . (F) Dose re-sponse curves for [ tG ∗ ] s.s depend on [ mGAP ](0) and does not change significantly as [ mG ](0) increases. (G) Colormap for [ tG ∗ ] s.s.; lower tG* concentrations result fromhigher [ mGAP ](0) values, since tGEF* and tG* are recruited for mGAP activation.Parameter values: k on = 3( s.µM s ) − , k o f f = 1( s.µM ) − , [ mG ∗ ](0) = 0 µM , [ tGEF tot ](0) = 10 µM , [ tGEF ∗ ](0) = 5 µM , T ∗ (0) = 0 . , [ tG tot ] = 10 µM , [ mGAP ∗ ](0) = 1 µM , [ tGAP ∗ ](0) = 1 µM , [ mGEF ∗ ] = 1 µM . Simulation times: s (B and E) and s (C, D, F, and G) . Numerical simulations were performed usingthe solver ode23s in Matlab R2018a. All parameters were arbitrarily chosen only toillustrate the dynamic features of the model. no mG and mG* concentrations. On the other hand, for [ mG ∗ ](0) > [ tGEF ](0) − [ mGAP ∗ ](0) = 5 µM , then the system achieves the steady state 2 with non zero concen-trations of both m and t-GTPase. Fig.6B shows the results for the feedforward connec-tion, and feedback loops tGEF → mGAP (arrows 1+2). In this particular example, the four steady states can be achieved for [ mGAP ∗ ](0) and [ mG ∗ ](0) ranging from 0 to 12 µM and 0 and 10 µM , respectively. In the vertical direction, the initial amount of mG*governs the transitions from steady states 3 to 4 (lower [ mGAP ∗ ](0) ) and 1 to 2 (higher [ mGAP ∗ ](0) ). In both steady states 2 and 4 (Eqs. 3.34 and 3.36), the concentrations ofmGTPase are nonzero. Therefore, we predict that an increase of initial concentration ofmG* would favor the emergence of these two steady states. In the horizontal direction,when the initial amount of mGAP* increases, the available mGAP (inactive) decreasesas we set the total mGAP as 12 µM , which reduces the effects of the feedback loops andthus facilitates the emergence of the steady states 1 and 2 where the concentrations oftGTPase are nonzero.Fig.6C shows a similar colormap for the system with both feedback loops tGEF → mGAP and tG ∗ → mGAP . It is worth noticing the expansion of the basin of attractionof Families 1 and 2 compared to Fig.6A, while the basin of Families 3 and 4 shrinks.Remarkably, in both Figs. 6B and 6C, there is a critical point (represented by a blackcross) of intersection of the four basins of attraction. In this case, disturbances in theinitial conditions around that intersection point can drive the system to different steadystates. Thus, while coupling the two GTPase switches with a forward arrow only givestwo possible steady states, the negative feedback afforded by arrows 2 and 3 give riseto a larger range of possibilities. Additionally, the existence of a critical point emergesin the presence of the negative feedback suggesting a rich phase space for this coupledsystem. 4. D ISCUSSION
GTP-binding proteins (GTPases) regulate crucial aspects of numerous cellular events.Their ability to act as biochemical switches is essential to promote information process-ing within signaling networks. The two types of GTPases - monomeric (m) and het-erotrimeric (t) - have traditionally believed to function independently. More recently, aseries of experimental findings revealed that m- and tGTPases co-regulate each other inthe Golgi through a functionally coupled circuit [32]. In this work, we sought to un-derstand the dynamic properties of this dual GTPase circuit. To this end, we developeda system of differential equations that describes the evolution of two coupled GTPaseswitches. Our analysis provided insight into the emergence of steady state configura-tions that cannot be achieved in systems of isolated GTPase switches. To the best of ourknowledge, this is the first modeling effort that described coupled GTPase switches witha feedforward connection mG ∗ → tGEF and the two feedback loops tGEF → mGAP and tG ∗ → mGAP [32, 61, 30, 62, 63, 31].A major result from our analysis is a systematic characterization of the steady state con-centrations of both m- and tGTPases, as well as their GEFs and GAPs. Given the modelformulation and the fact that we do not know the various kinetic parameters, obtaining TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 27 (A) (B) (C) s.s 4 Family 4 s.s 3 Family 3 s.s 1 Family 1 s.s 2 Family 2s.s 1 s.s 2 F IGURE Basins of Attraction – dependency on [ mG ∗ (0)] and [ mGAP ∗ (0)] (A) The two steady states of the system with feedforward con-nection (Section 3.1) are only driven by changes in the initial amount of mG*(B) When the feedforward and both feedback loops tGEF → mGAP are con-sidered, we observe the emergence of four regions (green,yellow, dark blue andlight blue colored) corresponding to the four steady states from Section 3.2 (C)A similar result was found when we analyzed the system with the feedforwardand both feedback loops tGEF → mGAP and tG ∗ → mGAP . A black crossindicates a critical point at the intersection of the four basins of attraction. Pa-rameter values: k on = 3( s.µM ) − , k off = 1( s.µM ) − , [ mG ](0) = 0 µM , [ mG ∗ ](0) = 0 µM , [ tG ](0) = 5 µM , [ tG ∗ ](0) = 0 µM , [ tGEF ](0) = 5 µM , [ mGAP ](0) = 12 µM − [ mGAP ∗ ](0) , [ tGAP ∗ ](0) = 1 µM , [ mGEF ∗ ](0) =1 µM F IGURE Main results and conclusions from steady state analysis
We performed a steady state analysis of a GTPase coupled circuit thathas been observed experimentally. For three biologically relevant combi-nations among the feedforward and two feedback loops, we present thesteady states and their interpretation. Moreover, we established the re-quired initial conditions for the existence of the steady states. Each con-nection adds to the richness of the functioning of these coupled GTPaseswitches. these states is critical. We show the obtained steady states in all three arrow combina-tions that we chose to analyze (Fig.7). Remarkably, the different steady states show avariety of configurations in which both m- and tGTPase can be interpreted as havinglow or high concentration values. These results help us to understand how the feed-forward connection and feedback loops impact the coupled system. Once we know thesteady states, locally asymptotic stability ensures that the trajectories will converge backto the steady state, given sufficiently close initial conditions. In §3.1 – §3.2, we con-firmed that all steady states obtained with a feedforward connection and feedback loop tGEF → mGAP (arrow 1 or arrows 1+2 in Fig.1B) are locally asymptotically sta-ble. However, when the two feedback loops are considered along with the feedforward(arrows 1+2+3 in Fig.1B), the local stability analysis cannot be performed because thesteady states are not isolated. Instead, we obtain four one-parameter families that dependon the amount of inactive tGTPase. At this point, further investigation would be neededto determine the behavior of the system near those steady state families. Even as we aimto develop complex models that are refined with iterative experimental validations, wenote that our analysis gives insight to different steady states that emerge due to differentcouplings that may not exist in physiology. Such insights may become meaningful inthe context of disease pathogenesis where copy numbers of each player in the networkmotif may change relative to each other, and do so dynamically (e.g., when respondingto stress/stimuli), or disease-driving mutations alter their functions (e.g., activating andinactivating mutations in GTPases, GAPs, or GEFs).Even though the mathematical structure of the model may appear simple, we find a richphase space for the coupled GTPase switches by analyzing the combination of networkconnections that have more biological meaning. Future studies could also explore therole of an external stimulus in the emergence of ultrasensitive behavior [4] and spatialorganization of these switches. Moreover, our system is well suited for coupling withexperimental measurements [38], including dose-response curves, response times, andnoise fluctuations, as done recently in [70]. TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 29
5. A
CKNOWLEDGMENTS
This work was supported by Air Force Office of Scientific Research (AFOSR) Multidis-ciplinary University Research Initiative (MURI) grant FA9550-18-1-0051 (to P. Ranga-mani) and the National Institute of Health (CA100768, CA238042 and AI141630 to P.Ghosh). Lucas M. Stolerman acknowledges support from the National Institute of Health(CA209891). R EFERENCES [1] B. Alberts, A. Johnson, J. Lewis, D. Morgan, M. Raff, K. Roberts, and P. Walter. Molecular biologyof the cell. 2015.
Garland, New York , pages 139–194, 2013.[2] H. E. Hamm. The many faces of g protein signaling.
Journal of Biological Chemistry , 273(2):669–672, 1998.[3] H. R. Bourne, D. A. Sanders, and F. McCormick. The gtpase superfamily: a conserved switch fordiverse cell functions.
Nature , 348(6297):125–132, 1990.[4] A. Lipshtat, G. Jayaraman, J. C. He, and R. Iyengar. Design of versatile biochemical switches thatrespond to amplitude, duration, and spatial cues.
Proceedings of the National Academy of Sciences ,107(3):1247–1252, 2010.[5] K. L. Rossman, C. J. Der, and J. Sondek. Gef means go: turning on rho gtpases with guaninenucleotide-exchange factors.
Nature reviews Molecular cell biology , 6(2):167–180, 2005.[6] J. Wang, Y. Tu, S. Mukhopadhyay, P. Chidiac, G. H. Biddlecome, and E. M. Ross. Gtpase-activatingproteins (gaps) for heterotrimeric g proteins.
G Proteins: Techniques of Analysis , pages 123–151,1999.[7] J. L. Bos, H. Rehmann, and A. Wittinghofer. Gefs and gaps: critical elements in the control of smallg proteins.
Cell , 129(5):865–877, 2007.[8] J. Cherfils and M. Zeghouf. Regulation of small gtpases by gefs, gaps, and gdis.
Physiological re-views , 93(1):269–309, 2013.[9] D. P. Siderovski and F. S. Willard. The gaps, gefs, and gdis of heterotrimeric g-protein alpha subunits.
International journal of biological sciences , 1(2):51, 2005.[10] P. Ghosh, P. Rangamani, and I. Kufareva. The gaps, gefs, gdis and now, gems: New kids on theheterotrimeric g protein signaling block.
Cell cycle , 16(7):607–612, 2017.[11] I. Lopez-Sanchez, Y. Dunkel, Y.-S. Roh, Y. Mittal, S. De Minicis, A. Muranyi, S. Singh, K. Shan-mugam, N. Aroonsakool, F. Murray, et al. Giv/girdin is a central hub for profibrogenic signallingnetworks during liver fibrosis.
Nature communications , 5(1):1–18, 2014.[12] G. S. Ma, N. Aznar, N. Kalogriopoulos, K. K. Midde, I. Lopez-Sanchez, E. Sato, Y. Dunkel, R. L.Gallo, and P. Ghosh. Therapeutic effects of cell-permeant peptides that activate g proteins down-stream of growth factors.
Proceedings of the National Academy of Sciences , 112(20):E2602–E2610,2015.[13] H. Wang, T. Misaki, V. Taupin, A. Eguchi, P. Ghosh, and M. G. Farquhar. Giv/girdin links vascularendothelial growth factor signaling to akt survival signaling in podocytes independent of nephrin.
Journal of the American Society of Nephrology , 26(2):314–327, 2015.[14] A. Hartung, A.-M. Ordelheide, H. Staiger, M. Melzer, H.-U. H¨aring, and R. Lammers. The akt sub-strate girdin is a regulator of insulin signaling in myoblast cells.
Biochimica et Biophysica Acta(BBA)-Molecular Cell Research , 1833(12):2803–2811, 2013. [15] V. DiGiacomo, M. Maziarz, A. Luebbers, J. M. Norris, P. Laksono, and M. Garcia-Marcos. Probingthe mutational landscape of regulators of g protein signaling proteins in cancer.
Science signaling ,13(617), 2020.[16] D. Hanahan and R. A. Weinberg. The hallmarks of cancer. cell , 100(1):57–70, 2000.[17] G. A. Cardama, N. Gonz´alez, J. Maggio, P. L. Menna, and D. E. Gomez. Rho gtpases as therapeutictargets in cancer.
International journal of oncology , 51(4):1025–1034, 2017.[18] W. N. Liu, M. Yan, and A. M. Chan. A thirty-year quest for a role of r-ras in cancer: from an oncogeneto a multitasking gtpase.
Cancer letters , 403:59–65, 2017.[19] K. Sriram, K. Moyung, R. Corriden, H. Carter, and P. A. Insel. Gpcrs show widespread differentialmrna expression and frequent mutation and copy number variation in solid tumors.
PLoS biology ,17(11):e3000434, 2019.[20] V. Wu, H. Yeerna, N. Nohata, J. Chiou, O. Harismendy, F. Raimondi, A. Inoue, R. B. Russell,P. Tamayo, and J. S. Gutkind. Illuminating the onco-gpcrome: Novel g protein–coupled receptor-driven oncocrine networks and targets for cancer immunotherapy.
Journal of Biological Chemistry ,294(29):11062–11086, 2019.[21] M. O’hayre, J. V´azquez-Prado, I. Kufareva, E. W. Stawiski, T. M. Handel, S. Seshagiri, and J. S.Gutkind. The emerging mutational landscape of g proteins and g-protein-coupled receptors in cancer.
Nature Reviews Cancer , 13(6):412–424, 2013.[22] M. Garcia-Marcos, P. Ghosh, and M. G. Farquhar. Giv is a nonreceptor gef for g α i with a uniquemotif that regulates akt signaling. Proceedings of the National Academy of Sciences , 106(9):3178–3183, 2009.[23] P. Ghosh. Heterotrimeric g proteins as emerging targets for network based therapy in cancer: End ofa long futile campaign striking heads of a hydra.
Aging (Albany NY) , 7(7):469, 2015.[24] M. M. Papasergi, B. R. Patel, and G. G. Tall. The g protein α chaperone ric-8 as a potential therapeutictarget. Molecular pharmacology , 87(1):52–63, 2015.[25] E. Evers, G. Zondag, A. Malliri, L. Price, J.-P. Ten Klooster, R. Van Der Kammen, and J. Collard. Rhofamily proteins in cell adhesion and cell migration.
European journal of cancer , 36(10):1269–1274,2000.[26] S. Etienne-Manneville and A. Hall. Rho gtpases in cell biology.
Nature , 420(6916):629–635, 2002.[27] Y. Takai, T. Sasaki, and T. Matozaki. Small gtp-binding proteins.
Physiological reviews , 81(1):153–208, 2001.[28] A. G. Gilman. G proteins: transducers of receptor-generated signals.
Annual review of biochemistry ,56(1):615–649, 1987.[29] F. A. Barr, A. Leyte, and W. B. Huttner. Trimeric g proteins and vesicle formation.
Trends in cellbiology , 2(4):91–94, 1992.[30] J. L. Stow, J. B. De Almeida, N. Narula, E. J. Holtzman, L. Ercolani, and D. A. Ausiello. A het-erotrimeric g protein, g alpha i-3, on golgi membranes regulates the secretion of a heparan sulfateproteoglycan in llc-pk1 epithelial cells.
The Journal of Cell Biology , 114(6):1113–1124, 1991.[31] J. Cancino and A. Luini. Signaling circuits on the g olgi complex.
Traffic , 14(2):121–134, 2013.[32] I.-C. Lo, V. Gupta, K. K. Midde, V. Taupin, I. Lopez-Sanchez, I. Kufareva, R. Abagyan, P. A. Ran-dazzo, M. G. Farquhar, and P. Ghosh. Activation of g α i at the golgi by giv/girdin imposes finitenessin arf1 signaling. Developmental cell , 33(2):189–203, 2015.[33] U. Alon.
An introduction to systems biology: design principles of biological circuits . CRC press,2019.[34] J. M. Bower and H. Bolouri.
Computational modeling of genetic and biochemical networks . MITpress, 2001.
TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 31 [35] N. J. Eungdamrong and R. Iyengar. Modeling cell signaling networks.
Biology of the Cell , 96(5):355–362, 2004.[36] A. E. Cowan, I. I. Moraru, J. C. Schaff, B. M. Slepchenko, and L. M. Loew. Spatial modeling of cellsignaling networks. In
Methods in cell biology , volume 110, pages 195–221. Elsevier, 2012.[37] M. K. Morris, J. Saez-Rodriguez, P. K. Sorger, and D. A. Lauffenburger. Logic-based models for theanalysis of cell signaling networks.
Biochemistry , 49(15):3216–3224, 2010.[38] M. Getz, L. Swanson, D. Sahoo, P. Ghosh, and P. Rangamani. A predictive computational modelreveals that giv/girdin serves as a tunable valve for egfr-stimulated cyclic amp signals.
Molecularbiology of the cell , 30(13):1621–1633, 2019.[39] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon. Network motifs: simplebuilding blocks of complex networks.
Science , 298(5594):824–827, 2002.[40] A. Goldbeter and D. E. Koshland. An amplified sensitivity arising from covalent modification inbiological systems.
Proceedings of the National Academy of Sciences , 78(11):6840–6844, 1981.[41] J. E. Ferrell Jr. Self-perpetuating states in signal transduction: positive feedback, double-negativefeedback and bistability.
Current opinion in cell biology , 14(2):140–148, 2002.[42] N. T. Ingolia and A. W. Murray. Positive-feedback loops as a flexible biological module.
Currentbiology , 17(8):668–677, 2007.[43] J. E. Ferrell Jr and W. Xiong. Bistability in cell signaling: How to make continuous processes dis-continuous, and reversible processes irreversible.
Chaos: An Interdisciplinary Journal of NonlinearScience , 11(1):227–236, 2001.[44] U. Alon. Network motifs: theory and experimental approaches.
Nature Reviews Genetics , 8(6):450–461, 2007.[45] S. R. Neves, P. T. Ram, and R. Iyengar. G Protein Pathways. 296(5573):1636.[46] S. S. Shen-Orr, R. Milo, S. Mangan, and U. Alon. Network motifs in the transcriptional regulationnetwork of escherichia coli.
Nature genetics , 31(1):64–68, 2002.[47] U. S. Bhalla and R. Iyengar. Emergent Properties of Networks of Biological Signaling Pathways.283(5400):381–387.[48] U. S. Bhalla and R. Iyengar. Emergent properties of networks of biological signaling pathways.
Sci-ence , 283(5400):381–387, 1999.[49] U. S. Bhalla and R. Iyengar. Robustness of the bistable behavior of a biological signaling feedbackloop.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 11(1):221–226, 2001.[50] E. A. Logsdon, S. D. Finley, A. S. Popel, and F. M. Gabhann. A systems biology view of blood vesselgrowth and remodelling.
Journal of cellular and molecular medicine , 18(8):1491–1508, 2014.[51] S. D. Finley, L. J. Broadbelt, and V. Hatzimanikatis. Computational framework for predictivebiodegradation.
Biotechnology and bioengineering , 104(6):1086–1097, 2009.[52] S. D. Finley and A. S. Popel. Effect of tumor microenvironment on tumor vegf during anti-vegftreatment: systems biology predictions.
JNCI: Journal of the National Cancer Institute , 105(11):802–811, 2013.[53] P. Yen, S. D. Finley, M. O. Engel-Stefanini, and A. S. Popel. A two-compartment model of vegfdistribution in the mouse.
PloS one , 6(11):e27514, 2011.[54] G. Hornung and N. Barkai. Noise propagation and signaling sensitivity in biological networks: a rolefor positive feedback.
PLoS Comput Biol , 4(1):e8, 2008.[55] S. Hooshangi, S. Thiberge, and R. Weiss. Ultrasensitivity and noise propagation in a synthetic tran-scriptional cascade.
Proceedings of the National Academy of Sciences , 102(10):3581–3586, 2005.[56] T. Shibata and K. Fujimoto. Noisy signal amplification in ultrasensitive signal transduction.
Proceed-ings of the National Academy of Sciences , 102(2):331–336, 2005. [57] J. M. Pedraza and A. van Oudenaarden. Noise propagation in gene networks.
Science ,307(5717):1965–1969, 2005.[58] L. Qiao, W. Zhao, C. Tang, Q. Nie, and L. Zhang. Network topologies that can achieve dual functionof adaptation and noise attenuation.
Cell systems , 9(3):271–285, 2019.[59] S. H. Strogatz.
Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, ChemistryAnd Engineering . Westview Press, first edition edition edition, 1994.[60] L. Perko.
Differential equations and dynamical systems , volume 7. Springer Science & BusinessMedia, 2013.[61] C. Jamora, P. A. Takizawa, R. F. Zaarour, C. Denesvre, D. J. Faulkner, and V. Malhotra. Regulationof golgi structure through heterotrimeric g proteins.
Cell , 91(5):617–626, 1997.[62] J. L. Stow and K. Heimann. Vesicle budding on golgi membranes: regulation by g proteins andmyosin motors.
Biochimica et Biophysica Acta (BBA)-Molecular Cell Research , 1404(1-2):161–171,1998.[63] J. Stow. Regulation of vesicle trafficking by g proteins.
Curr. Opin. Nephrol. Hypertens , 4:421–425,1995.[64] D. T. Gillespie. Deterministic limit of stochastic chemical kinetics.
The Journal of Physical ChemistryB , 113(6):1640–1644, 2009.[65] S. K. Hahl and A. Kremling. A comparison of deterministic and stochastic modeling approaches forbiochemical reaction systems: On fixed points, means, and modes.
Frontiers in genetics , 7:157, 2016.[66] J.-P. Changeux, J. Thi´ery, Y. Tung, and C. Kittel. On the cooperativity of biological membranes.
Proceedings of the National Academy of Sciences of the United States of America , 57(2):335, 1967.[67] I. Famili and B. O. Palsson. The convex basis of the left null space of the stoichiometric matrix leadsto the definition of metabolically meaningful pools.
Biophysical journal , 85(1):16–26, 2003.[68] P. Rangamani and L. Sirovich. Survival and apoptotic pathways initiated by tnf- α : Modeling andpredictions. Biotechnology and bioengineering , 97(5):1216–1229, 2007.[69] H. E. Nusse and J. A. Yorke.
Dynamics: numerical explorations: accompanying computer programdynamics , volume 101. Springer, 2012.[70] K. R. Ghusinga, R. D. Jones, A. M. Jones, and T. C. Elston. Molecular switch architecture drivesresponse properties. bioRxiv , 2020.
TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 33 A PPENDIX
A. P
ROOF OF P ROPOSITION (cid:91) [ mG ] , (cid:92) [ mG ∗ ] , (cid:99) G ∗ , and (cid:99) T ∗ satisfying the following system: k on [ mGEF ∗ ] (cid:91) [ mG ] − k off [ mGAP ∗ ] (cid:92) [ mG ∗ ] − k on [ tGEF tot ](1 − (cid:99) G ∗ ) (cid:92) [ mG ∗ ] = 0 (A.1) k on [ tGEF tot ] (cid:99) G ∗ (1 − (cid:99) T ∗ ) − k off [ tGAP ∗ ] (cid:99) T ∗ = 0 (A.2) k on (1 − (cid:99) G ∗ ) (cid:92) [ mG ∗ ] = 0 (A.3) (cid:91) [ mG ] + (cid:92) [ mG ∗ ] + [ tGEF tot ] (cid:99) G ∗ = C. (A.4) From Eq. A.3, and since we assume k on > , we must have (cid:92) [ mG ∗ ] = 0 or (cid:99) G ∗ = 1 . Thuswe divide the steady state analysis in two cases.Case 1: (cid:92) [ mG ∗ ] = 0 . From Eq. A.1 we must have (cid:91) [ mG ] = 0 and from Eq. A.4, we obtain (cid:99) G ∗ = C [ tGEF tot ] .Since (cid:99) G ∗ ≤ by definition, we conclude that C ≤ [ tGEF tot ] . (A.5)Eq. A.5 is also sufficient for (cid:92) [ mG ∗ ] = 0 . Otherwise, if C ≤ [ tGEF tot ] and (cid:92) [ mG ∗ ] > ,then (cid:99) G ∗ = 1 (Eq. A.3) and from Eq. A.4, we would conclude that (cid:91) [ mG ] + (cid:92) [ mG ∗ ] ≤ ,which is imposible.Finally, by substituting (cid:99) G ∗ in Eq. A.2, we obtain (cid:99) T ∗ = koff [ tGAP ∗ ] konC and therefore thesteady state is given by (cid:16) (cid:91) [ mG ] , (cid:92) [ mG ∗ ] , (cid:99) T ∗ , (cid:99) G ∗ (cid:17) = (cid:32) , ,
11 + k off [ tGAP ∗ ] k on C , C [ tGEF tot ] (cid:33) Case 2: (cid:99) G ∗ = 1 In this case, (cid:92) [ mG ∗ ] ≥ and from Eqs. A.1 and A.4,we obtain (cid:92) [ mG ∗ ] = k on [ mGEF ∗ ] k on [ mGEF ∗ ] + k off [ mGAP ∗ ] ( C − [ tGEF tot ]) and (cid:91) [ mG ] = k off [ mGAP ∗ ] k on [ mGEF ∗ ] + k off [ mGAP ∗ ] ( C − [ tGEF tot ]) . In this case, since the steady state has to be nonnegative, we must have C ≥ [ tGEF tot ] . (A.6)which is also sufficient for (cid:99) G ∗ = 1 . Otherwise if C ≥ [ tGEF tot ] and (cid:99) G ∗ < , then (cid:92) [ mG ∗ ] = (cid:91) [ mG ] = 0 (Eq. A.1) and, from Eq. A.3, we would have C = (cid:91) [ mG ] + (cid:92) [ mG ∗ ] + [ tGEF tot ] (cid:99) G ∗ < [ tGEF tot ] , which is impossible. Finally, by substituting (cid:99) G ∗ = 1 in Eq. A.2, we obtain k on [ tGEF tot ](1 − (cid:99) T ∗ ) − k off [ tGAP ∗ ] (cid:99) T ∗ = 0 which gives us (cid:99) T ∗ = k on [ tGEF tot ] k on [ tGEF tot ]+ k off [ tGAP ∗ ] and therefore (cid:16) (cid:91) [ mG ] , (cid:92) [ mG ∗ ] , (cid:99) T ∗ , (cid:99) G ∗ (cid:17) = (cid:18) k off [ mGAP ∗ ] k on [ mGEF ∗ ] + k off [ mGAP ∗ ] ( C − [ tGEF tot ]) ,k on [ mGEF ∗ ] k on [ mGEF ∗ ] + k off [ mGAP ∗ ] ( C − [ tGEF tot ]) ,k on [ tGEF tot ] k on [ tGEF tot ] + k off [ tGAP ∗ ] , (cid:19) . (A.7) TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 35 A PPENDIX
B. P
ROOF OF T HEOREM C , C , and [ mGAP tot ] for the existence of eachsteady state. We then compute the Jacobian matrix of the system and determine the localstability of the steady state based on the classical linearization procedure [59]. Steady states.
We divide our analysis into four different cases that emerge from thepreliminary inspection of the system given by Eqs. 3.27–3.32.Case 1: (cid:92) [ mG ∗ ] = 0 and (cid:92) [ mGAP ∗ ] = [ mGAP tot ] .From Eq. 3.27, we have (cid:91) [ mG ] = 0 and from Eq. 3.31, (cid:92) [ tGEF ∗ ] = C − [ mGAP tot ] .Thus C ≥ [ mGAP tot ] since the steady state must be nonnegative. Now Eq. 3.32 gives (cid:92) [ tGEF ] = C − C and that implies C ≥ C .Finally, Eq. 3.28 yields ( C − [ mGAP tot ])(1 − (cid:99) T ∗ ) − k [ tGAP ∗ ] (cid:99) T ∗ = 0 and hence (cid:99) T ∗ = C − [ mGAP tot ]( C − [ mGAP tot ]) + k [ tGAP ∗ ] The steady state is therefore given by (cid:98) x = (cid:18) , , C − [ mGAP tot ]( C − [ mGAP tot ]) + k [ tGAP ∗ ] , C − C , C − [ mGAP tot ] , [ mGAP tot ] (cid:19) . We now observe that the two parameter relations C ≥ [ mGAP tot ] and C ≥ C (B.1)are sufficient for (cid:92) [ mG ∗ ] = 0 and (cid:92) [ mGAP ∗ ] = [ mGAP tot ] . In fact, if C ≥ C then (cid:92) [ mG ∗ ] = 0 from the same argument as in Case 3. Now, Eq. 3.31 gives (cid:92) [ tGEF ∗ ] = C − (cid:92) [ mGAP ∗ ] and from Eq. 3.30, we must have (cid:92) [ mGAP ∗ ] = [ mGAP tot ] or (cid:92) [ tGEF ∗ ] = 0 .If (cid:92) [ tGEF ∗ ] = 0 then (cid:92) [ mGAP ∗ ] = C ≥ [ mGAP tot ] and hence (cid:92) [ mGAP ∗ ] = [ mGAP tot ] .Therefore, we have shown that Eq. B.1 imply (cid:92) [ mG ∗ ] = 0 and (cid:92) [ mGAP ∗ ] = [ mGAP tot ] .Consequently, the steady state in this case must be given by Eq. 3.33.Case 2: (cid:92) [ tGEF ] = 0 and (cid:92) [ mGAP ∗ ] = [ mGAP tot ] From Eq. 3.32, (cid:92) [ tGEF ∗ ] = C − [ mGAP tot ] and hence [ mGAP tot ] ≤ C . From Eq.3.31, we must have (cid:91) [ mG ] + (cid:92) [ mG ∗ ] = C − C and that implies C ≥ C . Now, Eq. 3.27gives [ mGEF ∗ ] (cid:16) C − C − (cid:92) [ mG ∗ ] (cid:17) = k [ mGAP tot ] (cid:92) [ mG ∗ ] and therefore (cid:92) [ mG ∗ ] = [ mGEF ∗ ] ( C − C )[ mGEF ∗ ] + k [ mGAP tot ] and (cid:91) [ mG ] = k [ mGAP tot ] ( C − C )[ mGEF ∗ ] + k [ mGAP tot ] . From Eq. 3.28, we must have ( C − [ mGAP tot ]) (1 − (cid:99) T ∗ ) − k [ tGAP ∗ ] (cid:99) T ∗ = 0 from which we get (cid:99) T ∗ = C − [ mGAP tot ]( C − [ mGAP tot ]) + k [ tGAP ∗ ] and therefore the steady state is given by (cid:98) x = (cid:18) k [ mGAP tot ] ( C − C )[ mGEF ∗ ] + k [ mGAP tot ] , [ mGEF ∗ ] ( C − C )[ mGEF ∗ ] + k [ mGAP tot ] , (B.2) C − [ mGAP tot ]( C − [ mGAP tot ]) + k [ tGAP ∗ ] , , C − [ mGAP tot ] , [ mGAP tot ] (cid:19) . We now observe that the two parameter relations C ≥ [ mGAP tot ] and C ≥ C (B.3)are sufficient for (cid:92) [ tGEF ] = 0 and (cid:92) [ mGAP ∗ ] = [ mGAP tot ] .In fact, if C ≥ C then (cid:92) [ tGEF ] = 0 from the same argument as in Case 1. Now, Eq.3.32 gives (cid:92) [ tGEF ∗ ] = C − (cid:92) [ mGAP ∗ ] and from Eq. 3.30, we must have (cid:92) [ mGAP ∗ ] =[ mGAP tot ] or (cid:92) [ tGEF ∗ ] = 0 . If (cid:92) [ tGEF ∗ ] = 0 then (cid:92) [ mGAP ∗ ] = C ≥ [ mGAP tot ] (fromEq. B.3) and thus (cid:92) [ mGAP ∗ ] = [ mGAP tot ] . Therefore, we have shown that Eq. B.3imply (cid:92) [ tGEF ] = 0 and (cid:92) [ mGAP ∗ ] = [ mGAP tot ] . Consequently, the steady state in thiscase must be given by Eq. 3.34.Case 3: (cid:92) [ mG ∗ ] = 0 and (cid:92) [ tGEF ∗ ] = 0 .From Eq. 3.27, we have (cid:91) [ mG ] = 0 and from Eq. 3.28, we also get (cid:99) T ∗ = 0 since k and [ tGAP ∗ ] are strictly positive numbers. Now, Eq. 3.31 gives (cid:92) [ mGAP ∗ ] = C and thuswe must have C ≤ [ mGAP tot ] . Moreover, Eq. 3.32 results in (cid:92) [ tGEF ] = C − C andsince all steady states must be nonnegative, we obtain C ≥ C . In this case, the steadystate is given by (cid:98) x = (0 , , , C − C , , C ) (B.4) We now observe that the two parameter relations C ≤ [ mGAP tot ] and C ≥ C (B.5)are sufficient for (cid:92) [ mG ∗ ] = 0 and (cid:92) [ tGEF ∗ ] = 0 . In fact, by subtracting 3.31 from Eq.3.32, we obtain (cid:92) [ tGEF ] − (cid:91) [ mG ] + (cid:92) [ mG ∗ ] = C − C ≥ TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 37 and hence (cid:92) [ tGEF ] ≥ (cid:91) [ mG ] + (cid:92) [ mG ∗ ] . On the other hand, from Eq. 3.29, we must have (cid:92) [ tGEF ] = 0 or (cid:92) [ mG ∗ ] = 0 . Thus if (cid:92) [ tGEF ] = 0 then (cid:91) [ mG ] + (cid:92) [ mG ∗ ] ≤ and hence thenonnegativeness of the steady state implies (cid:91) [ mG ] = (cid:92) [ mG ∗ ] = 0 . Hence we conclude that C ≥ C implies (cid:92) [ mG ∗ ] = 0 .Now, Eq. 3.31 gives (cid:92) [ tGEF ∗ ] = C − (cid:92) [ mGAP ∗ ] and from Eq. 3.30, we must have (cid:92) [ tGEF ∗ ] = 0 or (cid:92) [ mGAP ∗ ] = [ mGAP tot ] . If (cid:92) [ mGAP ∗ ] = [ mGAP tot ] , then (cid:92) [ tGEF ∗ ] = C − [ mGAP tot ] ≤ (from Eq. B.5) and thus (cid:92) [ tGEF ∗ ] = 0 . Therefore, we have shownthat Eq. B.5 imply (cid:92) [ mG ∗ ] = 0 and (cid:92) [ tGEF ∗ ] = 0 . Consequently, the steady state in thiscase must be given by Eq. 3.35.Case 4: (cid:92) [ tGEF ] = 0 and (cid:92) [ tGEF ∗ ] = 0 From Eq. 3.32, we obtain (cid:92) [ mGAP ∗ ] = C and hence C ≤ [ mGAP tot ] . From Eq. 3.31,we have (cid:91) [ mG ] + (cid:92) [ mG ∗ ] = C − C and that implies C ≥ C since the concentrations atsteady state must be nonnegative. Eq. 3.27 then gives − [ mGEF ∗ ] (cid:16) C − C − (cid:92) [ mG ∗ ] (cid:17) + kC (cid:92) [ mG ∗ ] = 0 from which we obtain (cid:92) [ mG ∗ ] = [ mGEF ∗ ] ( C − C )[ mGEF ∗ ] + kC and (cid:91) [ mG ] = kC ( C − C )[ mGEF ∗ ] + kC . From Eq. 3.28, we have (cid:99) T ∗ = 0 and therefore the steady state is given by (cid:98) x = (cid:18) kC [ mGEF ∗ ] + kC ( C − C ) , [ mGEF ∗ ][ mGEF ∗ ] + kC ( C − C ) , , , , C (cid:19) . We now observe that the two parameter relations C ≤ [ mGAP tot ] and C ≥ C (B.6)are sufficient for (cid:92) [ tGEF ] = 0 and (cid:92) [ tGEF ∗ ] = 0 . In fact, if C ≥ C then by subtractingEq. 3.32 from Eq. 3.31, we have (cid:91) [ mG ] + (cid:92) [ mG ∗ ] − (cid:92) [ tGEF ] = C − C ≥ and hence (cid:91) [ mG ] + (cid:92) [ mG ∗ ] ≥ (cid:92) [ tGEF ] . On the other hand, from Eq. 3.29, we must have (cid:92) [ tGEF ] = 0 or (cid:92) [ mG ∗ ] = 0 . Thus if (cid:92) [ mG ∗ ] = 0 then (cid:91) [ mG ] = 0 (from Eq. 3.27)and hence the nonnegativeness implies (cid:92) [ tGEF ] = 0 . Hence we conclude that Eq. B.6guarantee (cid:92) [ tGEF ] = 0 .Now, Eq. 3.32 gives (cid:92) [ tGEF ∗ ] = C − (cid:92) [ mGAP ∗ ] and from Eq. 3.30, we must have (cid:16) [ mGAP tot ] − (cid:92) [ mGAP ∗ ] (cid:17) = 0 or (cid:92) [ tGEF ∗ ] = 0 . If (cid:92) [ mGAP ∗ ] = [ mGAP tot ] then (cid:92) [ tGEF ∗ ] = C − [ mGAP tot ] ≤ (from Eq. B.6) and thus (cid:92) [ tGEF ∗ ] = 0 . Therefore,we have shown that Eq. B.6 implies (cid:92) [ tGEF ] = 0 and (cid:92) [ tGEF ∗ ] = 0 . Consequently, thesteady state in this case must be given by Eq. 3.36. Local Stability Analysis.
We begin reducing the ODE system with the conservationlaws given by Eqs. 3.19 and 3.20. In fact, if we write [ mG ] = C − [ mG ∗ ] − [ tGEF ∗ ] − [ mGAP ∗ ] and [ tGEF ] = C − [ tGEF ∗ ] − [ mGAP ∗ ] then Eqs. 3.21 – 3.26 can be written in the form d [ mG ∗ ] dt = f ([ mG ∗ ] , T ∗ , [ tGEF ∗ ] , [ mGAP ∗ ]) (B.7) d T ∗ dt = f ([ mG ∗ ] , T ∗ , [ tGEF ∗ ] , [ mGAP ∗ ]) (B.8) d [ tGEF ∗ ] dt = f ([ mG ∗ ] , T ∗ , [ tGEF ∗ ] , [ mGAP ∗ ]) (B.9) d [ tGEF ∗ ] dt = f ([ mG ∗ ] , T ∗ , [ tGEF ∗ ] , [ mGAP ∗ ]) (B.10) where f ([ mG ∗ ] , T ∗ , [ tGEF ∗ ] , [ mGAP ∗ ]) = k on [ mGEF ∗ ] ( C − [ mG ∗ ] − [ tGEF ∗ ] − [ mGAP ∗ ]) − k off [ mGAP ∗ ][ mG ∗ ] − k on ( C − [ tGEF ∗ ] − [ mGAP ∗ ]) [ mG ∗ ] ,f ([ mG ∗ ] , T ∗ , [ tGEF ∗ ] , [ mGAP ∗ ]) = k on [ tGEF ∗ ](1 − T ∗ ) − k off [ tGAP ∗ ] T ∗ ,f ([ mG ∗ ] , T ∗ , [ tGEF ∗ ] , [ mGAP ∗ ]) = k on ( C − [ tGEF ∗ ] − [ mGAP ∗ ]) [ mG ∗ ] − k on [ tGEF ∗ ] ([ mGAP tot ] − [ mGAP ∗ ]) , and f ([ mG ∗ ] , T ∗ , [ tGEF ∗ ] , [ mGAP ∗ ]) = k on ([ mGAP tot ] − [ mGAP ∗ ]) [ tGEF ∗ ] . The eigenvalues of the Jacobian Matrix can be thus calculated for each one of the foursteady states given by Eqs. 3.33 – 3.36. We prove that all steady states are LAS byshowing that the eigenvalues of the Jacobian Matrix are all negative real numbers. Weperform the calculations with MATLAB’s R2019b symbolic toolbox and analyze eachcase separately (see supplementary file with MATLAB codes). In what follows, let k mGEF = k on [ mGEF ∗ ] and k tGAP = k off [ tGAP ∗ ] .(1) If C > [ mGAP tot ] and C > C , the Jacobian matrix evaluated at the steadystate given by Eq. 3.33 gives the eigenvalues λ = − k on ( C − [ mGAP tot ]) and λ = − k on ( C − [ mGAP tot ]) − k tGAP which are negative. Moreover, the other eigenvalues λ and λ are such that λ + λ = k on ( C − C ) − k mGEF − k off [ mGAP tot ] < TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 39 and λ λ = − k mGEF k on ( C − C ) > and thus λ and λ are negative and hence the steady state is LAS.(2) If C > [ mGAP tot ] and C > C , the Jacobian matrix evaluated at the steadystate given by Eq. 3.34 gives the eigenvalues λ = − k mGEF − k off [ mGAP tot ] , λ = − k tGAP − k on ( C − [ mGAP tot ]) ,λ = − k on ( C − [ mGAP tot ]) and λ = − k on k mGEF ( C − C ) k off [ mGAP tot ] + k on [ mGEF ∗ ] which are all negative and hence the steady state is LAS.(3) If C < [ mGAP tot ] and C > C ,the Jacobian matrix evaluated at the steadystate given by Eq. 3.35 gives the eigenvalues λ = k on ( C − [ mGAP tot ]) and λ = − k tGAP which are negative. Moreover, the other eigenvalues λ and λ are such that λ + λ = k on ( C − C ) − C k off − k mGEF < and λ λ = − k mGEF k on ( C − C ) > and thus λ and λ are negative and hence the steady state is LAS.(4) If C < [ mGAP tot ] and C > C , the Jacobian matrix evaluated at the steadystate given by Eq. 3.36 gives the eigenvalues λ = − k mGEF − C k off , λ = k on ( C − [ mGAP tot ]) , λ = − k tGAP and λ = − k on k mGEF ( C − C ) C k off + k on [ mGEF ∗ ] which are all negative and hence the steady state is LAS. A PPENDIX
C. P
ROOF OF T HEOREM ξ -dependent families of steady states,where ξ ≥ represent the tG concentration. We also obtain necessary relationships forthe conserved quantities ˜ C , ˜ C , and [ mGAP tot ] , as well as admissible intervals for ξ thatguarantee the existence of nonnegative steady states.Case 1: (cid:92) [ mG ∗ ] = 0 and (cid:92) [ mGAP ∗ ] = [ mGAP tot ] .From Eq. 3.47, we have (cid:91) [ mG ] = 0 and subtracting Eq. 3.46 from Eq. 3.45, we get (cid:92) [ tGEF ] = ˜ C − ˜ C ≥ only if ˜ C ≥ ˜ C . Substituting (cid:92) [ tGEF ] on the conservation lawgiven by Eq. 3.46 and using Eq. 3.48 to write (cid:91) [ tG ∗ ] = ξ [ (cid:92) tGEF ∗ ] k [ tGAP ∗ ] , we obtain ξ + ξ (cid:92) [ tGEF ∗ ] k [ tGAP ∗ ] + ( ˜ C − ˜ C ) + [ mGAP tot ] = ˜ C and hence (cid:92) [ tGEF ∗ ] = (cid:16) ˜ C − [ mGAP tot ] − ξ ] (cid:17) k [ tGAP ∗ ] k [ tGAP ∗ ] + ξ only if ˜ C − [ mGAP tot ] ≥ ξ. Therefore, in this case the ξ -dependent family of steadystates is given by (cid:98) x ξ = (cid:32) , , ξ, (cid:16) ˜ C − [ mGAP tot ] − ξ ] (cid:17) ξk [ tGAP ∗ ] + ξ , ˜ C − ˜ C , (cid:16) ˜ C − [ mGAP tot ] − ξ ] (cid:17) k [ tGAP ∗ ] k [ tGAP ∗ ] + ξ , [ mGAP tot ] (cid:33) Case 2: (cid:92) [ tGEF ] = 0 and (cid:92) [ mGAP ∗ ] = [ mGAP tot ] Using Eq. 3.47 to write (cid:91) [ mG ] = k ( ˜ C − ξ )[ mGEF ∗ ] (cid:92) [ mG ∗ ] and subtracting Eq. 3.46 from Eq. 3.45,we obtain the expressions for [ mG ∗ ] and [ mG ] (cid:92) [ mG ∗ ] = ( ˜ C − ˜ C )[ mGEF ∗ ] k [ mGAP tot ] + [ mGEF ∗ ] and (cid:91) [ mG ] = ( ˜ C − ˜ C ) k ( ˜ C − ξ ) k [ mGAP tot ] + [ mGEF ∗ ] and clearly we must have ˜ C ≥ ˜ C . Now looking at Eq. 3.46 and substituting (cid:91) [ tG ∗ ] = (cid:92) [ tGEF ∗ ] ξk [ tGAP ∗ ] , we obtain (cid:92) [ tGEF ∗ ] = (cid:16) ˜ C − [ mGAP tot ] − ξ ] (cid:17) k [ tGAP ∗ ] k [ tGAP ∗ ] + ξ TABILITY ANALYSIS OF A CIRCUIT WITH DUAL GTPASE SWITCHES 41 only if ˜ C − [ mGAP tot ] ≥ ξ . Therefore, in this case the ξ -dependent family of steadystates is given by (cid:98) x ξ = (cid:32) ( ˜ C − ˜ C ) k [ mGAP tot ][ mGEF ∗ ] + k [ mGAP tot ] , ( ˜ C − ˜ C )[ mGEF ∗ ][ mGEF ∗ ] + k [ mGAP tot ] ,ξ, (cid:16) ˜ C − [ mGAP tot ] − ξ ] (cid:17) ξk [ tGAP ∗ ] + ξ , , (cid:16) ˜ C − [ mGAP tot ] − ξ ] (cid:17) k [ tGAP ∗ ] k [ tGAP ∗ ] + ξ , [ mGAP tot ] (cid:33) . Case 3: (cid:92) [ mG ∗ ] = 0 and (cid:92) [ tGEF ∗ ] = 0 From Eqs. 3.47 and 3.48, we have (cid:91) [ mG ] = 0 and (cid:91) [ tG ∗ ] = 0 , respectively. SubtractingEq. 3.46 from Eq. 3.45, in this case we get (cid:92) [ tGEF ] = ˜ C − ˜ C ≥ only if ˜ C ≥ ˜ C .Now, from the conservation law given by Eq. 3.45, we obtain (cid:92) [ mGAP ∗ ] = ˜ C − ξ and (cid:92) [ mGAP ∗ ] ∈ [0 , [ mGAP tot ]] only if max(0 , ˜ C − [ mGAP tot ]) ≤ ξ ≤ ˜ C . In this case,the ξ -dependent family of steady states is given by (cid:98) x ξ = (cid:16) , , ξ, , ˜ C − ˜ C , , ˜ C − ξ (cid:17) . Case 4: (cid:92) [ tGEF ] = 0 and (cid:92) [ tGEF ∗ ] = 0 Eq. 3.48 gives (cid:91) [ tG ∗ ] = 0 and the conservation law given by Eq. 3.46 yields [ mGAP ∗ ] =˜ C − ξ . Now using Eq. 3.47 to write (cid:91) [ mG ] = k ( ˜ C − ξ )[ mGEF ∗ ] (cid:92) [ mG ∗ ] , the conservation law givenby Eq. 3.46 gives (cid:92) [ mG ∗ ] = ( ˜ C − ˜ C )[ mGEF ∗ ] k ( ˜ C − ξ ) + [ mGEF ∗ ] and (cid:91) [ mG ] = ( ˜ C − ˜ C ) k ( ˜ C − ξ ) k ( ˜ C − ξ ) + [ mGEF ∗ ] and since [ mGAP ∗ ] ∈ [0 , [ mGAP tot ]] and the steady states must be nonnegative, wemust have max(0 , ˜ C − [ mGAP tot ]) ≤ ξ ≤ ˜ C ≤ ˜ C . The ξ -dependent familiy of steady states is therefore given by (cid:98) x ξ = (cid:32) ( ˜ C − ˜ C ) k ( ˜ C − ξ )[ mGEF ∗ ] + k ( ˜ C − ξ ) , ( ˜ C − ˜ C )[ mGEF ∗ ][ mGEF ∗ ] + k ( ˜ C − ξ ) , ξ, , , , ˜ C − ξ (cid:33)(cid:33)