The three-way switch operation of Rac1/RhoA GTPase-based circuit controlling amoeboid-hybrid-mesenchymal transition

Metastatic carcinoma cells exhibit at least two different phenotypes of motility and invasion - amoeboid and mesenchymal. This plasticity poses a major clinical challenge for treating metastasis, while its underlying mechanisms remain enigmatic. Transitions between these phenotypes are mediated by the Rac1/RhoA circuit that responds to external signals such as HGF/SF via c-MET pathway. Using detailed modeling of GTPase-based regulation to study the Rac1/RhoA circuit's dynamics, we found that it can operate as a three-way switch. We propose to associate the circuit's three possible states to the amoeboid, mesenchymal and amoeboid/mesenchymal hybrid phenotype. In particular, we investigated the range of existence of, and the transition between, the three states (phenotypes) in response to Grb2 and Gab1 - two downstream adaptors of c-MET. The results help to explain the regulation of metastatic cells by c-MET pathway and hence can contribute to the assessment of possible clinical interventions.


Rac1/RhoA Regulatory Circuits Construction
In the core regulatory circuit (Fig. 1b), Rac1-GTP and RhoA-GTP can mutually inactivate each other and also promote their own activation, thus introducing more nonlinearity in the circuit. Each of these regulations has been indicated experimentally, as shown in Supplementary Table S1. For most of them, they are discovered in breast cancer cells, but a few of them are found in the other cell lines. Here, we assumed that the circuit we built is conserved in most of human cancer cells. Therefore, the study of this circuit can provide some useful insights into amoeboid to mesenchymal transition for cancer cells in general.

Theoretical framework for small GTPase-based Regulatory Circuits
To study the dynamics of Rac1/RhoA regulatory circuit, we developed a computational model for small GTPase-based Regulatory Circuits (GBC). Here, we show in details how we derived the model from binding and unbinding reactions at molecular level. We first focused on the process of cycling for a typical Rho GTPase 16 (denoted as Rho) (Fig. 2a).
The reactions for this process are where Rho-DNA represents the gene encoding this Rho GTPase, and Rho-GTP, Rho-GDP and Rho-GDI stand for Rho GTPase bound to GTP, GDP and GDI respectively.
Reaction SR1 is for the production of Rho GTPase. Reaction SR2 represents our assumption that Rho protein first binds to GDP instead of GTP when it is translated.
where R f and R g represent the Rho GTPase without binding to any molecules and the gene encoding Rho GTPase respectively. In order to simplify the equations, we proposed three assumptions here: 2) The levels of each GEF and GAP are constant except the ones activated by the signals (I 1 , I 2 ) we discussed later. Therefore, we can also integrate them into the reaction rate constants (such as k 8

3)
We assume that GDP is present in abundant amount. Rho can bind to GDP to produce inactive GDP-bound form (Rho-GDP) as soon as Rho protein is produced. Thus With these assumptions, the model can be converted to: ' 3 and k 4 are binding and releasing rate constant for GDI respectively. k 7 and k 5 ' are intrinsic GTP hydrolysis and loading rate constants respectively, while k ' 8 and k ' 6 are activated GTP hydrolysis and loading rate constants due to GAPs and GEFs respectively. GAPs; or internal such as the auto-regulations for Rac1 and RhoA. Therefore, the total GTP hydrolysis and loading rate constants (described by B and J functions respectively as given below) can be divided into three parts: the intrinsic ones, the activated ones by the GAPs and GEFs that are independent of the regulations of signals, and the activated ones by the GAPs and GEFs that depend on the signals. The first two parts, which can be grouped as basal GTP hydrolysis and loading rate constants, are both constant. Yet, the latter activated rate constants ( k A ) are related with the level of signals (I) and can be described by Hill function as: where k A max is the maximum activated rate constant. K d is the dissociation constant which represents the threshold at which the activated rate ( k A ) is at half maximum value.
n is the Hill coefficient determining the steepness of the function. Thus, the activated GTP hydrolysis and loading rate constants by signal I 1 and I 2 can be expressed as , respectively. Considering these signals into the GBC model above (Supplementary Equation (2)), the signals-driven GBC model can be given as: where we use k 8 '' and k 6 '' to represent the signal-independent activated rate constants for GTP hydrolysis and loading, respectively. To simplify Supplementary Equation (4), we defined two functions: one is B function standing for total GTP loading rate constant and the other is J function standing for total GTP hydrolysis rate constant. Both of them depend on the signals.
Integrating Supplementary Equation (4) and (5), we can get a generic deterministic model (Supplementary Equation (6)) to describe the transition among GDI state, GDP state and GTP state of a typical Rho GTPase (Fig. 2a).
Utilizing above methods, we can deduce the deterministic model for the detailed circuit ( Fig. 2), shown below: where Here Then we can get By substituting these expressions back into Supplementary Equation (7), we simplified the original model to four equations as below: where the first part such as

Parameters estimation
The values of most of the parameters considered in our model are not known exactly, but we hereby explain how we estimated the parameter values.

The roles of auto-regulations on Rac1 and RhoA
As showed in Fig.1c, Rac1/RhoA regulatory circuit usually can be a three-way switch (Tristability). In a narrow range of parameters, this circuit can even be quadra-stable ( Supplementary Fig. S2a). This multistability are mainly due to the auto-regulations in each side of the circuit, which can increase the nonlinearity in the system and open up additional stable steady state(s) 29 . Without these auto-regulations, only bistability can be found at most (Supplementary Fig. S2b). Yet, having the auto-regulation on at least one of the GTPases possibly makes the circuit tristable ( Supplementary Fig. S2c,b). Thus, the auto-regulations in this circuit are essential to the multistable behavior.

Bifurcation diagrams for the circuit driven by Grb2 and Gab1 signals.
Here, we calculated the one-parameter bifurcation diagrams for this regulatory circuit driven only by Grb2 signal or Gab1 signal. In each diagram, the increasing signals finally drive the cells to mono-stable phases, namely M or A phenotype. Since Gab1 activate both Rac1 and RhoA, when Gab1 increases further, the cells finally are induced to A/M phenotype ( Supplementary Fig. S4e,f).