Cusp Bifurcation in Metastatic Breast Cancer Cells

Ordinary differential equations (ODEs) can model the transition of cell states 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. We have found that concepts in bifurcation theory may be applied to model metastatic cell behavior. Our results show how a specific phenomenon called a cusp bifurcation describes metastatic cell state transitions, separating two qualitatively different transition modalities. Moreover, we show how the cusp bifurcation models other genetic networks, and we relate the dynamics after the bifurcation to observed phenomena in commitment to enter the cell cycle.


Introduction
Cell types have been modeled as stable steady states, or equilibria, of ordinary differential equationsẋ = V (x, α) 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. Here, to identify steady states, we seek 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 α varies with time because of genetic, epigenetic, spacial effects or other factors, then we may have stable state x α(t) depending on α(t) describing the behavior of cells transitioning from one type to another. We will give examples coming from the study of metastatic breast cancer. However it should be noted that such analyses may be generalized to other biological networks such as cell division [15] or synthetic gene circuits [29].
Regarding metastasis, much attention has been devoted to two cellular states: epithelial (E) and mesenchymal (M), each recognizable by the levels of specific proteins [22,21,7,32,11]. Normal epithelial (E) cells are not motile and can grow (divide) in response to growth signals, as opposed to mesenchymal (M) cells that can detach from epithelial layers and are motile, but less likely to divide [24]. This choice between cell division and movement is the "go-or-grow hypothesis" [10]. 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 [1]. Thus, EMT alone is not always sufficient for metastasis [9] if a binary classification into E and M 3 cell types is assumed since MET is required for growth in the new location. Moreover, EMT might not be necessary either, according to recent work [9] [16] that showed that most metastatic cells do not undergo 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 [6].
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 [5,13,12]. One or more intermediate, "hybrid" or "partial" EMT cell states with mixed E/M properties have been described [14,2,4]. Recent computational work varied the number of hidden intermediate states, aiming to improve fits to experimental data [3]. The results indicated that intermediate states can accelerate EMT. It seems possible that these intermediate EMT states are sufficient for metastasis [8]. We are interested in studying the transitions which involve pro-metastatic and anti-metastatic mono-stables states of cells, analogous to EMT and MET, using mathematical models derived from bifurcation theory. In the modelẋ = V (x, α), the variables x(t) and parameters α(t) are time dependent and represent proteins involved in metastatic cell transitions. We decide to view the stable equilibria of the system as "cell types" as in Smale and Rajapakse [29] where they identify conditions in biological networks for which a pitchfork bifurcation in 2 and 3 variable systems with one parameter exist. The pitchfork bifurcation 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. We use known methods [28] to verify that a pitchfork exists in the metastasis model by Lee et al.[6], which was strongly suggested by Figure 1. Moreover we have verified that all the examples of pitchfork bifurcations proven by Smale and Rajapakse in [29] concerning "Repressillator" and "Toggle" synthetic gene circuits are actually one dimensional sub-bifurcations of cusp bifurcations. Furthermore, we plot the cusp curve which is the projection of the fold onto the parameter space and see that it divides regions in the parameter space. The regions separated by the curve and cusp point may correspond to biological transitions of two types. The first is what we call a continuous transition, as it relates to the biological phenomenon where there may be "hybrid" or "partial" cell types and occurs when there is only one stable equilibrium in the dynamical system. As we mentioned previously, this transition type may be sufficient for metastasis, however the state types and whether or not they are a continuum are unknown. The second transition we have called a discontinuous transition between binary cell types, or anti-metastatic and pro-metastatic cells, and this happens when a stable equilibrium bifurcates and two stable equilibria and one unstable equilibrium (saddle) appear.
Further experimentation and testing must be done to verify that this analysis of the mathematical dynamics of cusp bifurcations corresponds to the biological phenomenon. Additionally, we believe that these findings can be generalized and that cusp bifurcations may exist in models which describe other biological networks such as cell division [15]. In what follows we show a novel application of bifurcation theory in biology, address the question of non-binary transitions in metastasis, similar to the behavior of transition between E and M cells, and discuss further applications of these concepts. We have also provided an Appendix as a tutorial in bifurcation theory and any relevant code at the end of the paper.

Cusp bifurcation and metastatic cell transitions 2.1 Metastatic cell transition as an ODE model
In the model from the paper "Network of mutually repressive metastasis regulators can promote cell heterogeneity and metastatic transition" [6] Lee et al. considered three differential equations with two parameters V (R, L, B, ρ, k) where R, L, B are real positive variables and ρ, k are real positive parameters. R, L, B represent the proteins RKIP, let-7 and BACH1 which interact in the cell and are highly relevant for determining breast cancer metastasis. The parameters ρ and k describe the instability of RKIP and insensitivity of BACH1 to self-regulation, respectively. The equations The constants are set to s = .02, S = 20, c = 200, m = 2, b = 3, r = 5, p = 10, a = 1000. Let x = (R, L, B) so for convenience we may write V ( x, ρ, k) and take derivatives with respect to x.
Observable cell states are equilibria of the system V ( x, ρ, k). For ρ, k fixed, equilibria are points x such that V ( x, ρ, k) = 0. The equilibrium is stable if the real parts of the eigenvalues of D x V ( x, ρ, k) are negative. All solutions of the ODE which start near a stable equilibrium tend to the equilibrium as time increases. An equilibrium may lose its stability and a bifurcation may occur as we vary (ρ, k) if one of the eigenvalues tends to have zero real part or more specifically, if the eigenvalue becomes zero. In other words, to find bifurcations, we are looking for the solutions of the Via an intricate analysis of the ODE, Lee et al. divide a region in the (ρ, k) parameter plane into three sub-regions ( Figure 1). One region with a single stable equilibrium corresponding to an anti-metastatic state of the cell, one with a single stable equilibrium corresponding to a prometastatic state of the cell and one with three equilibria, two of which are stable. This suggests that as ρ, k vary the state of the cell may start in one mono stable region and pass through a bi-stable region to the other mono stable region. Thus the boundary separating the bi-stable and mono-stable regions is a curve of interest. It is defined by the solution of 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 under determined system [27] to (locally) locate the set of solutions. In this case this is a curve which has a tangent direction. The curve projected onto the (ρ, k)-plane is smooth and locally separates the regions if the tangent vector to the curve is not vertical. This type of point is called a saddle node point. If the tangent vector to the curve is vertical then the projection of the point to the (ρ, k)-plane is a cusp point. The shape of the bi-stable 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 D x V ( x, ρ, k). We assume ( v) is not zero. The rank of D x V ( x, ρ, k) is two where the equations are satisfied because D x V ( x, ρ, k) is a 3 by 3 matrix with an eigenvalue equal to zero. Since the map has maximal rank it is stable even if the parameters vary slightly.
The solution to the equations (1) define the values of proteins R, L, B and the parameter values ρ, k required at the cusp point. In order to observe metastasis, cells must traverse from the bottom to the top and back in Figure 1. To the left of the cusp, which is the grey area of Figure 1, we have cells that can transition continuously from one state to another ( Figure 2a). To the right of the cusp, cells must cross the bistable region in yellow in the upward direction and then again in the downward direction, where the cell plots are discontinuous ( Figure 2b). We explain this further in the Discussion section below. Mathematically, these processes can be computed by finding solutions 8 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 be 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.

Methods and model verification
We built a Newton's Method algorithm in MATLAB for an under determined system of equations to validate the ρ, k solution values identified by Lee et al. We used MATCONT [31], a continuation toolbox for ODEs, to isolate a 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 Appendix.     In Figure 7, 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 2a. 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. as demonstrated in Figure 5. Furthermore, if we project the fold lines of the surface onto the ρ, k plane we end up with a curve and cusp point as seen in Figure 4. Figure 9 illustrates these three perspectives and summarizes how the solution and cusp curves relate to the surface.

Gene circuits
The cusp point is relevant to many areas of biology besides metastasis. It should be found for gene networks described earlier as bistable or multistable. This applies to the toggle switch, one of the foundational gene circuits in synthetic biology. Thus, we considered the 2-gene and a 3-gene networks from Appendix 1A and 1C of Smale and Rajapakse [29], as other examples of a cusp bifurcation in biological systems. Smale and Rajapakse showed that pitchfork bifurcations exist for these networks. Below we give an analysis of the 2-gene network.
We used MATCONT to show that a cusp point occurs for the systeṁ In Figure 10(a) 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 10(b) we set α 1 = 2, m = k > 0 and varied α 1 and m. 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 [29]. 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 earlier in Figure 5. We plot the pitchfork using MATCONT in Figure 11.

15
In Figure 12 we used the same fimiplicit3 plot in MATLAB as we did for Figure 8. 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.

Cell Division
Another important network for which our analysis applies drives the cell cycle. Rajapakse and Smale [30] 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 [15]).
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 [29] 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 steady state corresponds to differentiation, which is known to be antagonistic with cell cycling [23] [20].
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 13 we produce a vector plot using the 2-gene Toggle network from Smale and Rajapakse [29]. 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. Figure 13: Vector plot of 2 gene circuit system to show 3 equilibria [29] 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 the 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 [29] follows:ẋ To find the behavior in the bistable region we set α = 2 and m = 3. In Figure 14 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 [15]). 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   14(b), 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 14(b), 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 [15]. In the examples the line in 2-dimension and the disc or eigenspace in 3-dimension approximate the true manifold 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 [29,15]. Furthermore, one could demonstrate that a pitchfork or cusp bifurcation exists in cell division as well.

Discussion
We have identified a cusp bifurcation in the model provided by Lee et. al. To observe such a bifurcation the critical points of a system must be stable. Here we will explain the mathematical process to find such a bifurcation as it applies to our model in some more detail.
In Figure 10 we have superimposed potential paths of α(t) in the space of parameters ρ and k from Figure 1.
The time-evolution of the parameters α(t) = (ρ, k) shows possible transitions between stable cell types. Here, the green line indicates how values of α(t) begin in the monostable anti-metastatic region (A), traverses an ambiguous monostable region and arrives into the monostable pro-metastatic region (P). These transitions we will call APT, or PAT depending on the direction. The unique stable state of the differential equations x(t) evolves together with α(t) but the possible cellular type ultimately transitions between pro-metastatic and anti-metastatic states. Therefore this image represents a continuous transition depending on which way the path is traversed. It is important to note that both transitions traverse the ambiguous region to the left of the bistable region. Cells in this ambiguous region are neither pro-nor anti-metastatic. They are in an intermediary biologically ambiguous state where the pro-and antimetastatic states become indistinguishable. Such ambiguous parameter regions must also appear in all other analyses describing bistable or multistable systems, but they have not been mathematically characterized from the perspective of metastasis. Analogies may exist with the physics of phase transitions, e.g., liquid and gas phases becoming indistinguishable beyond the critical point of water.
In Figure 11 the path of α(t) crosses the boundary between the monostable and bistable regions. Similar transitions have been extensively described in the EMT/MET literature. At first we consider the path as it goes up transitioning from anti-metastatic to pro-metastatic (APT). As the path crosses the boundary at t 0 , the edge of the bistable mixed region, the stable state x(t) continues on while a new equilibrium point is created which we denote x (t 0 ). Over time the point x (t) splits into a new stable point x 1 (t) and a saddle x 2 (t). In the bistable region there are three equilibria x(t), x 1 (t) and x 2 (t). If we cross the higher boundary of the bistable region a stable point and a saddle collide and annihilate each other that is, x(t) and x 2 (t) collide.
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.
Here is a picture in 1-dimensions of what is happening. Now if we imagine that metastasis will require an APT transition followed by a PAT transition we see that the PAT transition may require that the parameter cross the whole bistable region in reverse. The region near the star in Figure 1 would seem to provide the most fertile region for such a transition in both directions. This applies to both continuous and discrete transitions, to the left and right of the red star. We will call the star point the cusp point, the reason will become clear below. Two questions are immediate: 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. 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 binary transition may not sufficiently describe the biological phenomenon and we think that it is possible to show mathematically that these systems bifurcate.
Further interesting work would be to solve the equations in Lee et al. symbolically to make sure we have all possible solutions and it would be valuable to rigorously prove that the rank of the derivative is always 4. A forthcoming analysis of the initial results in Figures 13-16 could be made to generalize the observation of cusp bifurcations in other biological networks and extend the theory in a variety of biological phenomenon. Finally, in order to have confidence that the dynamics we described in this paper relates to metastatic transitions more experimentation is required to validate the theory. We thank Indika Rajapakse for helpful discussion as well as his software recommendation, and we thank him for pointing out some of the literature concerning cell division as well as open problems to us. We also thank Chris York's assistance in producing a final cusp fold image in MATLAB.

Appendix A: Cusp Bifurcation Methods
So far our analysis has been based on the rather fine numerical work of Lee et. al. [6] represented in Figure 1. 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 equationṡ x = 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 [18,25,17]. 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. [26,19] 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. 25 We assume one of the eigenvalues of the derivative D x V (x, α) is 0, so we add and now we assume that Here is an illustration of the saddle node bifurcation from Scholarpedia in one dimension and one parameter. Figure 18: Saddle-node bifurcation in the one-dimensional system [26] 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 19: Saddle-node bifurcation on the plane in the system [26] 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. 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 28 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.
29 Figure 21: Cusp bifurcation in the one-dimensional system [19] The surface below is the surface of equilibria. Outside of the cusp region there is one equilibrium and inside 3. Figure 22: Equilibrium manifold near a cusp bifurcation [19] As a path crosses the cusp curve away from the cusp point itself the equation has a saddlenode bifurcation. The picture of cubics illustrates again what is happening in 1-dimension. The two curves of the cusp are tangent at the cusp point. So the cusp point is a place where the hysteresis effect may in some sense be briefest. See questions 1 and 2 above.
A strategy for the analysis of equilibria might be to look for non -degenerate solutions of V (x, α) = 0 and Det[D x V (x, α)] = 0 for saddle node bifurcations and of the system of equations