Dynamic modelling of the mTOR signalling network reveals complex emergent behaviours conferred by DEPTOR

The mechanistic Target of Rapamycin (mTOR) signalling network is an evolutionarily conserved network that controls key cellular processes, including cell growth and metabolism. Consisting of the major kinase complexes mTOR Complex 1 and 2 (mTORC1/2), the mTOR network harbours complex interactions and feedback loops. The DEP domain-containing mTOR-interacting protein (DEPTOR) was recently identified as an endogenous inhibitor of both mTORC1 and 2 through direct interactions, and is in turn degraded by mTORC1/2, adding an extra layer of complexity to the mTOR network. Yet, the dynamic properties of the DEPTOR-mTOR network and the roles of DEPTOR in coordinating mTORC1/2 activation dynamics have not been characterised. Using computational modelling, systems analysis and dynamic simulations we show that DEPTOR confers remarkably rich and complex dynamic behaviours to mTOR signalling, including abrupt, bistable switches, oscillations and co-existing bistable/oscillatory responses. Transitions between these distinct modes of behaviour are enabled by modulating DEPTOR expression alone. We characterise the governing conditions for the observed dynamics by elucidating the network in its vast multi-dimensional parameter space, and develop strategies to identify core network design motifs underlying these dynamics. Our findings provide new systems-level insights into the complexity of mTOR signalling contributed by DEPTOR.


S1.1. Model reaction rates and equations
. Reactions and rates of the DEPTOR-mTOR models (see Fig. 1b, main text and Fig. S4a for the corresponding reaction schemes of the closed and open models, respectively).

Reactions
Reaction rates Parameter Values  S4a). Table S2. Ordinary differential equations (ODEs) of the DEPTOR-mTOR models. The reaction rates are given in Table S1.  S4a).

S1.2. Guidance for selection of parameter values
As precise values of the kinetic parameters in the DEPTOR-mTOR system are largely unknown at present for many different cell types, the parameter values used as the reference set for model analysis and simulations were guided by typical ranges of physiological values and constrained by data where possible (Table S1). For examples, protein dissociation constants in binding events typically lie in the low nanomolar range for strong bonds and in the low micromolar range for weak bonds. The kon rates are limited by the rate of collisions, which is limited by the rate of diffusion approximately ranging from 0.1 to 10 nM -1 s -1 1 . Michaelis-Menten (MM) constants (Km) can typically vary over a broad range and to scan a wide space we set if from 1 to 1000 nM. Catalytic constants range in MM (kc) was set from 0.0001 to 1 s -1 and maximal velocities (Vm) from 0.001 to 10 nMs -1 . Moreover, it is important to note that a major aim of modelling is to provide a basis for guiding experimental analysis and testing explicit hypotheses; a model by itself is not an objective "truth," but it can be used to falsify or confirm a specific hypothesis. Therefore, comprehensive systematic parameter exploration compatible with experimentally observed behaviour constitutes a suitable approach to mitigate the lack of measured parameters, and comprehensively characterise systems dynamic of the networks of interest. This is the approach we adopted in this study.
Plausible ranges of protein abundances of DEPTOR, mTORC1/2 and other proteins were based on those reported by quantitative proteomic studies previously reported for a range of cell lines, e.g. HeLa 2 and U2OS 3 cell lines. The cell volumes used for calculation in HeLa and U2OS cell lines were 2.6x10 3 μm 3 4 and 4x10 3 μm 3 3 respectively. For each protein, we divide the molecular weight (given in Da or g mol -1 ) by Avagadro's constant (6.02×10 23 mol -1 ) to obtain the mass per molecule in grams. This is then multiplied by the number of copies per cell to obtain the mass of protein per cell, m. Therefore, the molarity is obtained from the equation where p is the number of copies per cell, V is the cell volume (L) and A is Avagadro's constant. As a result, we typically varied the protein concentrations from 0 to 1000 nM. It is however very important to note that systems dynamics are not typically determined by the absolute abundance of the network nodes of but mainly by the relative abundances between them (e.g. if DEPTOR is 3 folds more than mTORC1). Thus, the fact that we varied protein concentrations within the range [0, 1000] (nM) for all nodes certainly captured all possible abundance ratios between the nodes.

S2.1. Bifurcation analyses in low-dimensional parameter spaces
To show that bistability exists for a wide range of parameters and is not specific to the reference parameter sets chosen for simulation, we constructed 1D and 2D bifurcation diagrams using the continuation software XPPaut, where XPP (X-windows phase space) provides a graphical interface to AUTO. XPPaut finds fixed points of a system and tracks them as a parameter is varied, giving lines along which equilibrium points exist for the set of parameters being examined 5,6 . Other temporal simulations are implemented using Wolfram Mathematica 10.1.0.0 (Wolfram Research, Inc).

S2.2. Multi-dimensional analysis of systems dynamics using DYVIPAC
Although the dynamical properties of a system are often examined based on conventional 2D bifurcation analysis (e.g. using XPPaut and AUTO 5,6 and presented on 2D bifurcation diagrams, the fact that usually only two parameters are varied at a time (while remaining parameters are set at fixed values) poses concrete limitation in our effort to obtain a global, multi-dimensional picture of the systems dynamics.
To overcome this limitation, we employ a software we recently developed called DYVIPAC 7 where multiple model parameters can be simultaneously sampled; and thus the dynamic behaviour of the studied system can be probed over a much wider region of the multidimensional parameter space. Importantly, DYVIPAC adapted the Parallel Coordinates plots 8 as an innovative way to represent the multi-dimensional data from the ensemble dynamical analysis 7 . The classification of the network's dynamic behaviour at each sampled parameter set by DYVIPAC is based on the concept of linear stability analysis, briefly explained below. For details, please consult 7 and other texts on this topic.

S2.3. Dynamical assessment based on linear stability analysis
Tools of nonlinear dynamics provide a useful framework to assess the dynamical properties of the DEPTOR-mTOR models. We are mainly interested in the asymptotic states and transitions between different dynamic regimes, such as bistable, oscillatory and fixed-point dynamics. The steady states of the system are obtained by equating the right hand side (RHS) of the system ODEs to zero and solving for the values of the model states' concentrations.
A solution different from the trivial one (all the species equal to zero) can exist and the temporal evolution of the system can be described by the independent selected species. If we assume that Sindep is the vector composed by these independent concentrations, the temporal evolution of the systems is completely determined by where f(Sindep) is the vector formed by the reaction rates of the independent selected species. The steady states denoted by S 0 indep are obtained by solving f(S 0 indep)=0. The asymptotic stability of these steady states can be determined by linear analysis upon perturbations. Thus, in the vicinity of any of these steady states the temporal evolution of a perturbation from this state, denoted by ΔS 0 indep =Sindep -S 0 indep, is given by where J(S 0 indep) is the Jacobian matrix of the reduced system evaluated at the considered steady state. Thus, the dynamical behaviour of the system is entirely specified by the Jacobian matrix and by its eigenvalues and eigenvectors. If we assume that the eigenvalues of J(S 0 indep) for a given steady state S 0 indep, denoted by λi (i=1,2,..n where n is the number of independent species) are ordered decreasingly by the values of their real part, Re(λi), then the dominant growth term of the perturbation is governed by where E1 is the eigenvector corresponding to the eigenvalue with the largest real part. The state S 0 indep is asymptotically stable if, and only if, Re(λ1) is negative and ΔS 0 indep tends to exponentially decrease in the course of time. If there is at least one positive eigenvalues' real part (i.e. at least Re(λi) is positive) then the perturbation grows exponentially and S 0 indep is unstable in response to the perturbation.
In order to study the asymptotic states and the transition between different dynamical behaviours, DYVIPAC numerically solves the system for any set of parameter values given. Once the solutions are obtained, they are substituted in the Jacobian Matrix and the eigenvalues are numerically calculated to classify the different dynamical behaviours. When a unique steady state exists, the sign of the largest real part of the eigenvalue associated with that solution allows us to classify that steady state as stable or unstable indicating the presence of sustained oscillations. The same analysis is possible when three steady states coexist. Evaluating the Jacobian matrix and calculating the eigenvalues for each solution enabled classification of the possible dynamical behaviours. For the considered ranges of parameters, different cases were reported: (1) two stable and one unstable solution indicating a bistable dynamical behaviour; and (2) one stable and two unstable solutions suggesting a co-existing bistable/oscillatory behaviour where the system is basically bistable with one fixed-point steady-state branch and one periodic branch. Figures   Fig. S1. Related to Fig. 3c, main text. Dependence of steady-state pmTORC2 and pmTORC1 (shown in Fig. 3c) on increasing DEPTOR abundance, V1 = 3. The remaining parameter values used are given in Table S1 and Table S2.  Fig. 6f but for bistable sets (b), co-existing BS/OS (c) and overlay of all three dynamics (d).

S3. Supplementary
Red ovals indicate strong restriction of parameter values that underlying a specific dynamics. Fig. 6 and Figs. 5b, e, main text. (a-b) PC plots showing oscillation and bistability parameter sets resulting from DYVIPAC-based analysis for combined changes of DEPTOR's affinities with mTORC1 and 2 (k13f, k14f) and strength of mTORC1-induced negative feedback loop (k15c) On short timescales (0-60 minutes) the system that includes protein synthesis and degradation behaves almost identically to the system where synthesis and degradation are neglected. In contrast, on long timescales (>> 1 hrs) when synthesis and degradation are included, a unique steady state (whose value depends on the synthesis and degradation rates) is reached. Fig. S5. Illustration of a therapeutically-relevant expression window for DEPTOR, within this range DEPTOR level is sufficiently abundant to effectively inhibit mTORC1 and mTORC2, but not too high to trigger hyperactivation of Akt, another pro-suvival kinase. In contrast DEPTOR levels below this range do not effectively suppress mTORC1 or mTORC2; and DEPTOR levels higher than this range unexpectedly activate Akt.  Here we chose to display representative sets with qualitatively different behaviours for illustrative purpose: (i) In the first parameter set (first row), both models display two separate steady states (one fixed-point, one periodic) at short timescale but in the open model, pmTORC1 converges to a single steady state in the long timescale; (ii) In the second parameter set, both models display two fixed-point steady states at short timescale, but the open model converges to a single fixed-point. However this only happens at a very long time (>40hrs). (iii) In the third parameter set, both models are al most indistinguisable over a relatively long time (>12hrs). Taken together, both models behave similarly over short timescale for a variety of parametric conditions. Rank is the number of elements in the largest linearly independent set of reaction vectors iii)

Fig. S3. Related to
Deficiency is a specialized measure of the linear independence of the reaction vectors defined as: (number of complexes) -(number of linkage classes) -(rank of network)