Cusp bifurcation in a metastatic regulatory network

Understanding the potential for cancers to metastasize is still relatively unknown. While many predictive methods may use deep learning or stochastic processes, we highlight a long standing mathematical concept that may be useful for modeling metastatic breast cancer systems. Ordinary differential equations (ODEs) can model cell state transitions by considering the pertinent environmental variables as well as the paths systems take over time. Bifurcation theory is a branch of dynamical systems which studies changes in the behavior of an ODE system while one or more parameters are varied. Many studies have applied concepts in one-parameter bifurcation theory to model biological network dynamics, and cell division. However, studies of two-parameter bifurcations are much more rare. Two-parameter bifurcations have not been studied in metastatic systems. Here we show how a specific two-parameter bifurcation phenomenon called a cusp bifurcation separates two qualitatively different metastatic cell state transitions modalities and propose a new perspective on defining such transitions based on mathematical theory. We hope the observations and verification methods detailed here may help in the understanding of metastatic potential from a basic biological perspective and in clinical settings.


Introduction
Different cell states can emerge during disease progression, such as cancer metastasis. Regarding metastatic cancer, much attention has been devoted to two cellular states: epithelial (E) and mesenchymal (M), each recognizable by the levels of specific proteins, which correspond to the steady states of multistable gene regulatory networks [31,29,14,42,18]. Normal epithelial (E) cells are not motile and can grow (divide) in response to growth signals, as opposed to mesenchymal (M) cells that do not form epithelial layers and are motile, but less likely to divide [34]. The assumption of a binary choice between cell division and movement is the "go-or-grow hypothesis" [17]. Since metastasis requires departure from a primary site and growth at a different site, both the epithelialmesenchymal transition (EMT), and mesenchymal-epithelial transition (MET) seem to be required [4]. Thus, EMT alone is not always sufficient for metastasis [16] if a binary classification into E and M cell types is assumed since MET is also required for growth in the new location. Moreover, EMT might not be necessary either, according to recent work [16] [23] showing that most metastatic cells do not undergo full EMT, and EMT inhibition does not reduce metastasis. Furthermore, metastasis relies on detached cells invading their tissue neighborhood and accessing the bloodstream, which occurs on top of EMT, under the control of different genes called pro-and antimetastatic regulators, such as BACH1 and RKIP [13]. We believe that we can pinpoint the threshold when these regulators cause a qualitative shift in the transition between the two states, and describe a new view that may be crucial to understanding whether and how metastasis will occur.
Many recent studies indicate that EMT and MET are more complex than binary processes, i.e., they are transitions between more than two distinct, well-defined cellular states [12,20,19]. One or more intermediate, "hybrid" or "partial" EMT cell states with mixed E/M properties have been described [21,9,11]. Recent computational work that varied the number of hidden intermediate states, aiming to improve fits to experimental data [10] found that intermediate states can accelerate 4 EMT. It seems possible that only these intermediate EMT states, instead of full EMT, are necessary and sufficient for metastasis [15]. In  Finally, how can these theoretical questions help us to understand the biology of metastasis?
The use of equilibrium states of differential equations as a model for biological or chemical systems including metastasis-regulatory networks, has a long history which we do not try to survey here, except for a few instances. [46,7,43,44,32,5,45,8] The study of bistable systems has played an important role, especially noteworthy are the works on the toggle switch [5,45] and on direct and indirect fully positive feedback [3] [30]. (A deeper example of these models and how the relate to cusp bifurcations can be found in the Supplementary Materials Appendix A) Mathematical studies of genetic regulatory networks have usually relied on solving the corresponding differential equations numerically or with the use of topological or fixed point techniques and theorems from one-parameter bifurcation theory to prove the existence of solutions with particular properties. One example of the use of topological methods describing critical points and associated Boolean networks is by Glass [24]. An example of the extensive use of numerics is Lee et al [13]. By contrast, the use of two-parameter bifurcation theory rarely if ever enters into the biological theory, even in papers with the word bifurcation in the title as in Spencer et al. [22] or in the text as in Rajapakse et al. [39] The cusp bifurcation is a concept from two-parameter bifurcation theory, found by solving a system of equations in the state and parameter variables. Here we present methods for finding cusp points and give a self-contained elementary derivation of a set of equations and their solution for finding cusp points. We provide below a simple computational framework to find accurate solutions to such systems. We show how this provides a broader view on the emergence of bistability in biological systems, by dividing the two-dimensional parameter space into regions with distinct transition types: continuous and discontinuous transitions. This should improve the conceptual understanding of cell state transitions in metastatic gene networks and other biological systems.

The cusp point separates two kinds of cell state transitions
Cell states or "cell types" [39] have been modeled as stable steady states, or equilibria, of such that x ∈ R n is a concentration vector representing the cell's molecular composition and α ∈ R j is a vector of j parameters, representing internal and external characteristics, including reaction rates. Steady states correspond to the values of x α which satisfy V (x α , α) = 0. If the eigenvalues of the derivative of V (x α , α), denoted D x V (x α , α), all have negative real part then nearby solutions all tend to x α as time increases. It is possible that V (x, α) may have a unique stable state, multiple stable states or even more complicated dynamical behavior. Each stable equilibrium x α,i corresponds to a cell type. If the parameter α depends on some external factor f ext because of sensory, genetic, epigenetic, spatial or other effects, then the number and value of stable states x α(fext) may depend on α(f ext ), which describes the dynamical behavior of cells transitioning from one type to another dependent on f ext , which could be the time variable. We are interested in studying the transitions which involve pro-metastatic and anti-metastatic mono-stable states of cells, analogous to EMT and MET, using mathematical models derived from bifurcation theory. In the modelẋ = V (x, α), the variables x(t) that represent proteins involved in metastatic cell transitions are time-dependent. The parameters α(f ext ) depend on an external factor, which could be a chemical concentration, cell size, or time. As long as α(f e xt) is a continuous function of f e xt, the theory we present is valid. If the 6 factor f e xt = t is time, as we assume in the following for simplicity, the adjustment time scale of the steady state would be faster than the time scale of α(t), as usually assumed in bifurcation theory. We will mainly focus on cell types in metastatic breast cancer, but such analyses may also be generalized to other metastatic cancers, biological networks orchestrating events such as cell division, [22] or synthetic gene circuits [39].
Smale and Rajapakse refer to cell states as "cell types" in [39] where they identify conditions in biological networks for which a pitchfork bifurcation in 2 and 3 variable systems with one parameter exist. Yet, the pitchfork is a one-parameter bifurcation that may not be stable. There is a stable two parameter bifurcation called the cusp bifurcation, which includes a pitchfork as a one dimensional sub-bifurcation. The applicability of the cusp bifurcation or two-parameter bifurcations to metastatic transitions has not been widely investigated. Here we use known methods [38] [39] concerning "Repressillator" and "Toggle" synthetic gene circuits are actually one dimensional sub-bifurcations of cusp bifurcations (See Supplementary Materials). We plot the cusp curve in the metastatic breast cancer model, which is the projection of the fold onto the parameter space and observe that it divides the parameter space into two regions separated by the curve and cusp point, which correspond to biological transitions of two types. These two transition types correspond to the two types of paths taken by the curve α(t). The first goes around the the cusp point, and has no bistable points, whereas the second crosses the bistable region. The first we call a continuous transition, as it may capture the biological phenomenon of not just one, but any number of "hybrid" or "partial" cell types, differently than in previous studies. Specifically, instead of multiple equilibria, we find that there is always only one stable equilibrium in the dynamical system, corresponding to a continuum of partial cell types, which may be sufficient for metastasis. The second transition is a 7 discontinuous transition between binary cell types, or anti-metastatic and pro-metastatic cells, which happens when a stable equilibrium bifurcates and two stable equilibria and one unstable equilibrium (saddle) appear.
Overall, we show a novel application of bifurcation theory in biology, propose a shift from continuous (non-binary) to discrete (binary) transitions at the cusp point, and discuss further applications of these concepts in metastasis and other biological phenomena. We have also provided a tutorial in bifurcation theory and any relevant code in the Supplementary Materials, which others may find useful for their own research.
There are now 4 equations in 5 unknowns to solve. We assume that the rank of the derivative of this system is 4 when the equations are satisfied. By the implicit function theorem we may graph a curve for an underdetermined system [37] to (locally) locate the set of solutions. The curve projected onto the (ρ, k)-plane is smooth and locally separates the regions of the plane. Under certain conditions this curve may meet at a cusp point (See Appendix B). The shape of the bistable region is such that we suspected that there is a cusp point, which we do in fact find to be close to the red star in Figure 1.
A cusp point is the (ρ, k) parameter coordinates of a non-degenerate solution of the following five by five system of equations.
Here ∇ x denotes the gradient in terms of x which we take of the determinant of the derivative.
We then calculate the dot product with ( v) which is the first column of the adjugate matrix of Mathematically, these paths can be observed by finding solutions of the system where either k or ρ are fixed. In the first case, k is fixed and must be less than its cusp point value, and in the second k must become greater than the cusp point value. In Figure 2 we draw sketches of what we expect these solutions to look like and how they relate to metastasis.

Numerical approach for finding the cusp point and model verification
We built a Newton's Method algorithm in MATLAB for an underdetermined system of equations to validate the ρ, k solution values we identified from Lee et al. We used MATCONT [41], a continuation toolbox for ODEs, to isolate a more accurate set of values for the cusp point on the (ρ, k)-plane and plot the projection of the cusp curve. We also built codes in MATLAB to visualize the various behaviors of the system which can be found in the supplementary materials.
To validate this is in fact a cusp point, we solved the system of equations in MATCONT to find that the cusp is located at ( x, ρ, k) = (0.9321, 2.2184, 0.0435, 1.0281, 0.1343). Note that the values from MATCONT return ρ, k = (1.0281, 0.1343) whereas the previous values in Figure 1 were   We explore how varying values of ρ affects B in Figure 5. First we note that near the cusp point for different values of k > .13 we see a similar curve. If we follow a solution from the lower monostable region it disappears as ρ increases, at 1.024. This is the first limit point (LP) or equilibrium in the graph. Then B increases towards the value corresponding to the other stable equilibrium which is the upper limit point. The curve between the two limit points represents the bistable region which we showed in Figure 2b. The bistable region in Figure 1 corresponds to the unstable equilibrium and the space between the two limit points. At the second limit point B   was fixed and ρ varied. In both plots the red stars are limit points, desmonstrating discontinuous solutions to the system.
In Figure 6, if k = .13 or less than the value required at the cusp, the solution curve is S shaped and B is increasing continuously as ρ increases. As we continuously increase ρ we pass from the anti-metastatic to pro-metastatic cell states as in Figure 2b. In this scenario, if points on the solution curve correspond to various cell types, we could posit that we transition to many intermediate cell states as we go from one mono-stable state to another. This differs from the case when k > .13 where we see two distinct sets of cell types. Figure 6: A continuous solution curve plot from MATCONT's Equilibrium curve differential equations solver. This image is meant to demonstrate how varying ρ for a fixed value of k < .13 will yield a continuous uninterrupted curve without any cusp or limit points.  in metastasis more generally, where the path of the system will determine final metastatic states. It may even be possible to find this behavior in other cancers or caner-related phenomena governed by bistable or multistable regulatory networks. Below we describe one possible interpretation of how the cusp may determine certain distinct system pathways towards metastasis.
In Figure 8 we have superimposed potential paths of α(t) in the space of parameters ρ and k from Figure 1. This would be a discontinuous transition in the direction A to P, since at every point along the path where the A and P states co-exist, they correspond to separate sets of variables x A and x P , which are clearly distinguishable. There is no ambiguity at any point about an individual cell being in one state or another. But now as we run this path backwards to produce a PAT we see that the state of the cell corresponds to x 1 (t) as x(t) tends to x 1 (t) in the bistable region.
Thus traversing the path in one direction and then the other exhibits hysteresis. Hysteresis is another hallmark of discontinuity. Unlike continuous transitions, discontinuous EMT/MET has been extensively investigated for many regulatory networks (See also Figure 2).
Here is a picture in 1-dimensions of what is happening. 1. What determines the path α(t) of continuous or discontinuous transition that the cell will trace out in the parameter space?
2. Will the paths of the cells stay away from or pass close to the cusp point? What is the biological significance of this alternative? Will paths that pass close to the cusp be more likely to undergo an APT transition, followed by an PAT transition, and thus be more likely to establish distant tumors?
We addressed these questions in connection with metastatic cell state transitions above. We also proposed the possibility of seeing the same behavior in 2 or 3-gene circuits and in cell division dynamics that mimic a "Toggle" network (See supplementary materials). A general and in depth perspective on the mathematics behind finding a cusp bifurcation is introduced in the appendix which we hope could be applied to other biological networks. It is possible even more biological systems exhibit cusp and pitchfork bifurcations where a simple one-parameter binary transition may not sufficiently describe the biological phenomenon and we think that it is possible to show mathematically that these systems can undergo two-parameter bifurcations. If this is the case then we believe it is important to revisit existing studies of gene networks and apply these findings from two-parameter bifurcation theory to system models and predictive techniques.
We note that there is another approach to the use of bifurcation theory in biology. This is the Catastrophe theory of Rene Thom [44]. The developed theory concerns the zeros of gradient vector fields and their bifurcations. There is a cusp bifurcation which is very much the same in its geometric features. However, one has to be a little careful here since the bifurcation theory of differential equations and gradient differential equations have some subtle differences. For example the general cusp bifurcation has the possibility of exhibiting a periodic solution which the gradient system does not. The survey paper by Rand et al [6] has recent updates to this theory. Given a differential equation satisfying certain properties, there is a gradient system which shares the 22 asymptotic behavior of the original system. Now Rand adds the unstable manifold geometry behavior to the analysis. A drawback of this theory may be that the gradient system is not immediately at hand, whereas we work with the equations directly.
Further work as it relates to metastatic breast cancer would be to incorporate these models to existing predictive processes that assess metastatic potential. Furthermore, it would be valuable   We used MATCONT to show that a cusp point occurs for the systeṁ Figure 11a we set α 2 = 2 and m = k > 0. As we varied α 1 and m we found a cusp exists at (x, y, α 1 , m) = (1.00, 1.00, 2.00, 1.99). In Figure 11b we set α 1 = 2, m = k > 0 and varied α 1 and m.

Acknowledgments
A cusp point exists at (x, y, α 2 , m) = (1.00, 1.00, 2.00, 1.99). This aligns with Smale and Rajapakse's results [39]. They set α 1 = α 2 = 2, m = k > 0 and proved that a pitchfork bifurcation occurs at the point m = 2. They note that for all 0 ≤ m < 2 there is only one solution to the one parameter system which is (x, y) = (1, 1). It can be shown that for values past the cusp, or when m > 2, there are three solutions and we have the same pitchfork scenario as we described in the paper. We plot the pitchfork using MATCONT in Figure 12.
In Figure 13 we used the same fimiplicit3 plot in MATLAB as we did in the paper. We plot the solution surface of x and y separately and vary the two relevant parameters (α 1 , m) and (α 2 , m) respectively. The surfaces both are folds which we would expect to see for a cusp bifurcation. Figure 11 shows the projection of the fold on the (α 1 , m) and (α 2 , m) plane. As for the 25 metastasis-regulatory network, the cusp point separates the parameter ranges for continuous versus discontinuous transitions.

Cell Division
Another important network for which our analysis applies drives the cell cycle. Rajapakse and Smale [40] suggest that the pitchfork bifurcation may explain some of the properties of cellular division. The decision of cells between cycling and senescence is driven by a bistable network, the core of which is a variant of the toggle switch involving CDK2 and Rb (Figure 2.A in [22]).
Considering the above analysis for the toggle switch, the cell cycle network might also be modelled as a Toggle, with a pitchfork or cusp bifurcation. Thus, we use the 2-gene Toggle network from Smale and Rajapakse [39] to demonstrate what could be happening during cell division and how it relates to a cusp bifurcation. We use the toggle switch to illustrate some biologically relevant features of the dynamics and note that the bifurcations of the cell cycle network will occur at different, currently unknown parameter values.
Since the toggle switch has two steady states, one of which corresponds to entering the cell cycle, we need to assign a biological meaning to the other steady state. We propose that the other 26 steady state corresponds to differentiation, which is known to be antagonistic with cell cycling [33] [28].
The dynamics of a differential equation with fixed parameters after the pitchfork or cusp bifurcations are: 1. There are 3 equilibria, two sinks and a saddle.
2. The forward orbit of any point tends to one of the equilibria.
3. The orbits which tend to the saddle, called the stable manifold of the saddle, form an (n-1) dimensional disc which separates the two basins of the stable equilibria.
In Figure 14 we produce a vector plot using the 2-gene Toggle network from Smale and Rajapakse [39]. We see all three behaviors listed in this image, where the upper sink could be low CDK2/high Rb, the lower sink is high CDK2/low Rb and the x = y line divides paths a cell may take as they differentiate if they initiate above the line or proliferate if they begin on a path below the line. 27 Almost all orbits tend to one of the stable equilibria. An orbit starting near the red line will have very slow dynamics. Paths close to the unstable equilibria and far from the stable equilibria will take a long time to arrive near a stable equilibrium. In the case of cell division one stable equilibrium corresponds to the cell being proliferative, and the other to cell differentiation. Points close to the stable manifold will have fast dynamics as they tend to the sink. However, points close to the unstable equilibria and far from the stable equilibria will define cells that may appear quiescent for a long time, although ultimately they will commit to cycling. Cells distributed along the green arrows could have drastically long waiting times before entering the cell cycle, as observed in the experiments.
We see the same behavior in 3 dimensions. The 3 dimensional system is defined by Smale and Rajapakse as [39] follows:ẋ To find the behavior in the bistable region we set α = 2 and m = 3. In Figure 15 we show the 3 equilibria as black stars, and the forward orbits are the blue dotted lines each tending towards one of the equilibria. The eigenspace of the saddle is a plane that divides two regions in three dimensions, separating the two basins of stable equilibria. This plane is an approximation of the stable manifold of the saddle. We use ode45 in MATLAB to solve the system starting from various points in the space. Most points stay on the same side of the manifold as where they started. We suggest that the same behavior is seen in the system involving CDK2, Rb and Ef2 ( Figure   2.A in [22]). In our model we may interpret x, y, z as EF2, CDK2, Rb respectively. In this space as z increases x, y decrease and solutions tend to the stable sink in the upper left hand corner of Figure   15b, which corresponds to differentiation. Conversely, as z decreases, x, y increase and solutions tend to the sink in the lower right hand corner of Figure 15b, which corresponds to cell cycle entry.
This reflects a switch behavior in cell division with 3 variables.
There is a symmetry in the 2-gene model when m > 2, α 1 = α 2 . In two dimensions the line x = y is the stable manifold that separates the two cell behaviors. Any initial point with x > y or x < y will converge to the same stable equilibirum. Similarly when considering the 3 dimensional model the stable manifold of the saddle separates the space such that all orbits tend to the sink of the same side. While our model is the linear approximation of the stable manifold it is clear that for values of x that are large enough, all orbits tend to the stable state which corresponds to high CDK2 concentrations. This is the observed situation with CDK2 and proliferation [22]. In the examples the line in 2-dimension and the disc or eigenspace in 3-dimension approximate the true manifold 29 and are based off simple gene networks, however, it could be interesting to develop a more accurate model of cell division to combine the findings in the two papers referenced [39,22]. Furthermore, one could demonstrate that a pitchfork or cusp bifurcation exists in cell division as well.

Appendix B: Cusp Bifurcation Methods
So far our analysis has been based on the rather fine numerical work of Lee et. al. [13] represented in Figure 1 of the paper. But we may give a more analytical treatment of the figures which in this case and much more generally establish the existence of the cusp point and give a sharper description of the geometry of the division into monostable and multistable regions near it. The discussion of how the stable states change along paths in the parameter space will be the same in this greater generality. The subject matter is called bifurcation theory. We are given smooth differential equationsẋ = V (x, α) as above where x ∈ R n and α is a parameter in R j and V : R n × R j → R n . In our case α is a real parameter so j = 1.
Bifurcation theory studies how the behavior of the solution of V changes as α changes. This is a very well developed subject [26,35,25]. The two stable bifurcations which will interest us are the saddle-node and cusp bifurcations. We will use the referenced numerical methods to find them. [36,27] The equation 1) V (x, α) = 0 defines the equilibria. Recall that in general V (x, α) : R n × R → R n and we assume that when V (x, α) = 0 equation 2) DV (x, α) goes from R n+1 → R n and the rank of the derivative is n.
So the set of zeros is locally a curve near (x, α) by the implicit function theorem.
We assume one of the eigenvalues of the derivative D x V (x, α) is 0, so we add Here is an illustration of the saddle node bifurcation from Scholarpedia in one dimension and one parameter. Figure 16: Saddle-node bifurcation in the one-dimensional system [36] If n = 2 or higher, here is the dynamical picture in the plane. In general the downward arrow represents an (n − 1) dimensional plane of the non-zero eigenvalues. Figure 17: Saddle-node bifurcation on the plane in the system [36] Now if we have two parameters ρ and k When we project the curve on the parameter plane, (ρ, k) there is either A) no singularity of the projection, i.e. no tangent vector to the curve is vertical or B) there is a singularity of the projection, that is a tangent vector to the curve which is vertical.
In case A) there are two possibilities: This is called a fold bifurcation.
Ab)There is a point along the curve where the rank of D x V (x, ρ, k) is n − 1 but a second eigenvalue is 0. Now when a curve crosses the projected curve we may create a saddle of index j and another equilibria of index j − 1 or j + 1 depending on where the crossing takes place and for some values of the parameter there are limit cycles. We don't see this behavior in the data of Lee et al. and we will not discuss it. This bifurcation is called a Bogdanov-Takens bifurcation.
There are generically five types of bifurcations of equilibria depending on two parameters (See Kuznetsov). All might be interesting for biology.
Case B) is the one that interests us. It is the cusp bifurcation. There is a singularity of the projection, that is a tangent vector to the curve which is vertical. Recall that the tangent to the 33 curve is the null space of the derivative of (V (x, ρ, k), Det[D x V (x, ρ, k)]). Now for a vertical vector (v, 0, 0) to be in the null space it is necessary and sufficient that The references above give methods for determining that these equations have a solution.
When equations 1 ) and 3 ) are satisfied D x V (x, ρ, k) has a kernel as Det[D x V (x, ρ, k)] = 0 so we would like to express the vector v in terms of (x, ρ, k). Assuming that D x V (x, ρ, k) has rank n − 1 when V (x, ρ, k) = 0 then at least one of the (n − 1) by (n − 1) minors has non-vanishing determinant and the adjugate matrix is not zero. By Cramer's Rule any non-zero column of the adjugate matrix is in the kernel of D x V (x, ρ, k).
So a cusp point is defined as a non-degenerate solution of the n + 2 by n + 2 system of equations V (x, ρ, k) = 0 Where v is a non-zero column of the adjugate matrix of D x V (x, ρ, k).
Here are some pictures of the cusp bifurcation when n = 1. Figure 19: Cusp bifurcation in the one-dimensional system [27] The surface below is the surface of equilibria. Outside of the cusp region there is one equilibrium and inside 3.