Pseudo-nullclines enable the analysis and prediction of signaling model dynamics

A powerful method to qualitatively analyze a 2D system is the use of nullclines, curves which separate regions of the plane where the sign of the time derivatives is constant, with their intersections corresponding to steady states. As a quick way to sketch the phase portrait of the system, they can be sufficient to understand the qualitative dynamics at play without integrating the differential equations. While it cannot be extended straightforwardly for dimensions higher than 2, sometimes the phase portrait can still be projected onto a 2-dimensional subspace, with some curves becoming pseudo-nullclines. In this work, we study cell signaling models of dimension higher than 2 with behaviors such as oscillations and bistability. Pseudo-nullclines are defined and used to qualitatively analyze the dynamics involved. Our method applies when a system can be decomposed into 2 modules, mutually coupled through 2 scalar variables. At the same time, it helps track bifurcations in a quick and efficient manner, key for understanding the different behaviors. Our results are both consistent with the expected dynamics, and also lead to new responses like excitability. Further work could test the method for other regions of parameter space and determine how to extend it to three-module systems.

Then, with the time derivative for Cdk1a equal to zero and replacing Cdk1i with the above, an expression that depends only on Cdk1a and Apca is obtained.By numerically scanning Cdk1a, an implicit equation for Apca is left.We solve this with MATLAB, obtaining different pairs (Apca,Cdk1a).This is the pseudo-nullcline for the first module, a result of the two differential equations for the Cdk complexes taken at zero.Now, we apply the same procedure to the Plxa and Apca equations.Taking the time derivative for Plxa equal to zero, we arrive at: Replacing Plxa with the above in the equation for Apca, we arrive at an explicit equation for Apca that depends only on Cdk1a.Once again, we numerically scan Cdk1a, obtaining different pairs (Apca, Cdk1a).This is the pseudo-nullcline for the second module, a result of the two differential equations for Plxa and Apca taken at zero.Apca and Cdk1a are the coupling variables of this system.
We take one final step before proceeding.Adding the equations for the Cdk complexes, we obtain: The sum of the two complexes defines the total cyclin (Cyctot), with above being its differential equation.When equal to zero: Now, we can perform a change of variables from Apca to Cyctot after arriving at the different pairs (Apca, Cdk1a).The variable Cyctot is the input of the system, and Cdk1a is the output responsible for mitotic phosphorylation.
Parameter values for the Tsai et.al model, taken from their work (Tsai, Theriot, and Ferrell 2014): Stable steady states in blue, unstable in red, maxima and minima of limit cycles in black circles.Both bifurcations are supercritical Hopf.Near the left-hand Hopf, the amplitude increases rapidly after a short parameter range with small values, as seen in the inset.
For the addition of the extra parameter: And with the differential equation equal to zero: The 5 conserved quantities are: We define X as:

𝒅𝒕
. The last two terms are zero, following the quasisteady state approximation.Taking the differential equation for K, one can re-write it as it follows: The last two terms are auxiliary.From quasi-steady state approximation, the differential equations for the first-level intermediate complexes are equal to zero, which leads to: Where   = .
Similar steps are applied to the second level, arriving at two differential equations for A and App.
To recap: The resulting four-dimensional system is: Setting the first two differential equations to zero, we have the pseudo-nullclines for X and K0.From each, we arrive at two expressions for Z, the coupling variable.Taking both equal to one another, the following explicit equation for K0 is obtained (we use the same notation for the variables and their concentrations): With MATLAB, we perform a scan of K1, obtaining K0.At this point, two unknown variables remain: X and Z.Using the conservation equation for Ktot, replacing each intermediate complex with equations resulting from the quasi-steady state approximation (see above), and replacing Z, an implicit equation for X is reached.After solving it with MATLAB, all that remains is the calculation of Z from either of the previously mentioned expressions.
For the second-level pseudo-nullcline, we rewrite Z as the combination of A and Ap, use the differential equations for A and App equal to zero, plus the conservation equation of Atot.This last equation can be written entirely as a function of Ap and App.Doing a scan of App, we find Ap.Then, we obtain A, which leads to Z and finally to X.
Parameter values for the 2+2 model: In the first level, every rate constant is equal to 1 with the following exceptions: l1=8, l6=50.
The input E1tot is changed and informed in each result.
Supplementary Figure 2. One-dimensional bifurcation diagram for the 2+2 model with the total amount of input kinase E1tot as the parameter, and the double-phosphorylated substrate App as the output.Stable steady states in blue, unstable in red, maxima and minima of limit cycles in black circles.Very small input ranges with oscillations are found starting from the Hopf bifurcations (H), which we did not consider in the results.
Basic scheme for the method, representing the two modules and their interconnections.(2), one has also that: Thus X(a * ) and Y (b * ) constitute a steady state of the coupled modules since the following equations are satisfied: In conclusion the intersections of the pseudo-nullclines C 1 and C 2 are in one-toone correspondence with the steady states of the coupled modules.Another advantage of this geometrical method is that it is able to reveal a limit point bifurcation, like a saddle-node bifurcation.As it shown in the S.I., this occurs when a steady state corresponds to a tangential intersection of the pseudonullclines.In particular this feature enables to distinguish between a SNIC bifurcation or a Hopf bifurcation because in the first case oscillations appear through a tangent bifurcation, whereas in the second case the pseudo-nullclines intersect transversely.Both cases are illustrated by applying our method to different signaling motifs studied in the result section. [

Tangency and bifurcation
In this section of the Supplementary Material, we prove the following result: if a stationary state of the system corresponds to a limit point bifurcation, then the pseudo-nullclines C 1 and C 2 become tangent at their intersection.
Referring to the main text, the starting hypothesis is that the stationary states of the system are solutions of the equations: where, to fix the ideas, x and f were supposed to be n-dimensional, and y and g to be m-dimensional.Moreover, we assume that a given solution of this system of n + m equations can be written as : X(y 1 ) − x = 0 and Y (x 1 ) − y = 0.One can consider these two equations as the zero of a unique function H(x, y) = 0, where H is a n + m-dimensional function.Now, given a stationary state (x * , y * ), it corresponds to a limit point bifurcation if the derivative DH(x * , y * ) is not invertible.(Otherwise, by the implicit function theorem, the solution (x * , y * ) can be uniquely continued, and there is no limit point bifurcation).DH(x * , y * ) is not invertible if and only if the determinant of the jacobian matrix DH(x * , y * ) equals zero.Therefore, let us compute the determinant of: The result can be decomposed in 2 terms A + B: the first, A, equals to the element DH 11 (equal to −1) times the corresponding minor.But this minor is the determinant of an upper triangular matrix having only elements −1 on its diagonal.Thus, the result is A = (−1) n+m .The second term, B, can be computed as DH 1,n+1 (= X ′ 1 ) times the cofactor (−1) n times the corresponding minor, that is the determinant of the following matrix: This result can be obtained by suppressing the line and the column of the element Y ′ 1 (with a cofactor (−1) n−1 ), leaving a (n + m − 4) diagonal matrix of −1.Therefore, canceling the determinant of matrix ( 4) is equivalent to writing: On the other hand, by construction the pseudo-nullclines C 1 and C 2 are parametrized by (X 1 (y 1 ), y 1 ) and by (x 1 ,Y 1 (x 1 )).So, at a given intersection (x * 1 , y * 1 ) the components of their tangent vectors can be computed respectively as (X ′ 1 , 1) and (1,Y ′ 1 ).Therefore canceling the determinant: gives a condition for the alignement of these tangent vectors, and thus for the tangency of the pseudo-nullclines at their intersection.
One-dimensional bifurcation diagram for the Tsai.et al model with ksynth as the parameter, and cdk1a as the output, using the authors parameter values.
Therefore in the plane of coordinates (a, b), steady states (a * , b * ) lie in the intersection set of two curves C 1 and C 2 whose graphs are respectively given by the parametrizations (A(b), b) and (a, B(a)).This remark motivates the appellation of "pseudo-nullclines" for these curves.Conversely, if (a * , b * ) belongs to the intersection of the two curves C 1 and C 2 , then a * = A(b * ) and b * = B(a * ).Using functions X and Y defined by eqs.