# Cusp Bifurcation in Metastatic Breast Cancer Cells

Brenda Delamonica  Applied Mathematics and Statistics Department, Stony Brook
University, Stony Brook, NY 11794, USA  Gábor Balázsi∗ The Louis and Beatrice
Laufer Center for Physical and Quantitative Biology, Stony Brook University,
Stony Brook, NY 11794, USA
Department of Biomedical Engineering Department, Stony Brook University, Stony
Brook, NY 11794, USA  Michael Shub∗∗ Department of Mathematics,City College
and the Graduate Center of CUNY

###### Abstract

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.

$\ast,\ast\ast$ Corresponding authors<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Contents

  1. 1 Introduction
  2. 2 Cusp bifurcation and metastatic cell transitions
    1. 2.1 Metastatic cell transition as an ODE model
    2. 2.2 Methods and model verification
    3. 2.3 Results and cusp point analysis
    4. 2.4 Cusp Bifurcation in other networks
      1. 2.4.1 Gene circuits
      2. 2.4.2 Cell Division
  3. 3 Discussion
  4. 4 Acknowledgments
  5. 5 Appendix A: Cusp Bifurcation Methods
  6. 6 Appendix B: Matlab Code

##  1 Introduction

Cell types have been modeled as stable steady states, or equilibria, of
ordinary differential equations

$\dot{x}=V(x,\alpha)$

such that $x\in\mathbb{R}^{n}$ is a concentration vector representing the
cell’s molecular composition and $\alpha\in\mathbb{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_{\alpha}$ which satisfy $V(x_{\alpha},\alpha)=0$. If the eigenvalues of the
derivative of $V(x_{\alpha},\alpha)$, denoted $D_{x}V(x_{\alpha},\alpha)$, all
have negative real part then nearby solutions all tend to $x_{\alpha}$ as time
increases. It is possible that $V(x,\alpha)$ may have a unique stable state,
multiple stable states or even more complicated dynamical behavior. Each
stable equilibrium $x_{\alpha,i}$ corresponds to a cell type. If the parameter
$\alpha$ varies with time because of genetic, epigenetic, spacial effects or
other factors, then we may have stable state $x_{\alpha(t)}$ depending on
$\alpha(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 epithelial-mesenchymal
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 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]. In general, the
number of such intermediate states is unknown, raising the question: Can the
number of intermediate states go to infinity, allowing continuous transitions?
And what defines the boundary between such continuous transitions versus the
widely-studied discrete transitions, with a finite number of distinguishable
steady states?

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
$\dot{x}=V(x,\alpha)$, the variables $x(t)$ and parameters $\alpha(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.

##  2 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,\rho,k)$ where $R,L,B$ are real positive variables and $\rho,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 $\rho$ and $k$ describe the instability of
RKIP and insensitivity of BACH1 to self-regulation, respectively. The
equations are

$\displaystyle\frac{dR}{dt}$ $\displaystyle=\frac{1}{1+B}-\rho R$
$\displaystyle\frac{dL}{dt}$ $\displaystyle=\frac{aR^{r}}{m^{r}+R^{r}}-L-cLB$
$\displaystyle\equiv V(R,L,B,\rho,k)$ $\displaystyle\frac{dB}{dt}$
$\displaystyle=s+\frac{(S-s)k^{b}}{k^{b}+B^{b}}-B-cLB$

The constants are set to $s=.02,S=20,c=200,m=2,b=3,r=5,p=10,a=1000$. Let
$\vec{x}=(R,L,B)$ so for convenience we may write $V(\vec{x},\rho,k)$ and take
derivatives with respect to $\vec{x}$.

Observable cell states are equilibria of the system $V(\vec{x},\rho,k)$. For
$\rho,k$ fixed, equilibria are points $\vec{x}$ such that
$V(\vec{x},\rho,k)=0$. The equilibrium is stable if the real parts of the
eigenvalues of $D_{x}V(\vec{x},\rho,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 $(\rho,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 determinant
$Det[D_{x}V(\vec{x},\rho,K)]=0$.

Figure 1: Map of dynamics from Lee et al. [6]. Grayscale shading indicates
monostability. Darker shading corresponds to higher BACH1 levels. Yellow color
indicates bistabiltiy. The red star indicates the expected cusp point where
($\rho$,K) = (1.0024,0.6595).

Via an intricate analysis of the ODE, Lee et al. divide a region in the
$(\rho,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 pro-metastatic
state of the cell and one with three equilibria, two of which are stable. This
suggests that as $\rho,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

$\displaystyle V(\vec{x},\rho,k)=0$ $\displaystyle
Det[D_{x}V(\vec{x},\rho,k)]=0$

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 $(\rho,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 $(\rho,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 $(\rho,k)$ parameter coordinates of a non-degenerate
solution of the following five by five system of equations.

$\begin{split}V(\vec{x},\rho,k)&=0\\\ Det[D_{x}V(\vec{x},\rho,k)]&=0\\\
\nabla_{x}Det[D_{x}V(\vec{x},\rho,k)]\bullet(\vec{v})&=0\end{split}$ (1)

Here $\nabla_{x}$ denotes the gradient in terms of $\vec{x}$ which we take of
the determinant of the derivative. We then calculate the dot product with
$(\vec{v})$ which is the first column of the adjugate matrix of
$D_{x}V(\vec{x},\rho,k)$. We assume $(\vec{v})$ is not zero. The rank of
$D_{x}V(\vec{x},\rho,k)$ is two where the equations are satisfied because
$D_{x}V(\vec{x},\rho,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 $\rho,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 of the system where either $k$ or $\rho$ 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.

(a) Path of $\alpha$ outside of bistable region

(b) Path of $\alpha$ through bistable region

Figure 2: Two different modes of cell state transitions between monostable
regions

###  2.2 Methods and model verification

We built a Newton’s Method algorithm in MATLAB for an under determined system
of equations to validate the $\rho,k$ solution values identified by Lee et al.
We used MATCONT[31], a continuation toolbox for ODEs, to isolate a cusp point
on the $(\rho,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.

###  2.3 Results and cusp point analysis

To validate this is in fact a cusp point, we used Newton’s Method and verified
$(\vec{x},\rho,k)=(0.9321,2.2184,0.0435,1.0281,0.1343)$ is a solution of the
equations (1). In Figure 3 the results of the Newton’s Method algorithm show
various iterations where all the values are very close to 0. When plotting
$\rho,k$ we found a similar separation of the regions defined by Lee et al.
(Figure 4).

Figure 3: Results using Newton’s Method algorithm

(a) MATCONT

(b) Newtons Method

Figure 4: Cusp Point plots

Next, we numerically verified that the five by five system has an invertible
derivative at the solution. At a cusp point the tangent to the cusp exhibits a
pitchfork bifurcation. [28]

Figure 5: Pitchfork Bifurcation

Figure 5 is a plot of the pitchfork bifurcation around the cusp point
$(\rho,k)=(1.0281,0.1343)$ where the tangent vector is (-0.1542,1). For values
of $T$ moving in the positive direction of the tangent, $B$ consistently has 3
solutions, two stable equilibria which are the upper and lower branches of the
pitchfork and one unstable equilibria in between. In the negative direction of
$T$ we find only 1 solution.

We explore how varying values of $\rho$ affects $B$ in Figures 6 and 7. First
we note that near the cusp point for different values of $k>.13$ we see a
similar curve (Figure 6). If we follow a solution from the lower monostable
region it disappears as $\rho$ 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$ crosses from the bistable to the monostable region. As $\rho$
increase $B$ stays in the upper branch of the curve. As $\rho$ is approaching
zero $B$ tends towards the upper limit point. The limit points are where
$\rho$ crosses the cusp line for fixed $k$, thus demonstrating hysteresis. We
may interpret this as having a concentration of antimetastatic cells for
values of $\rho=\\{0,1.024\\}$ and a concentration of prometastatic cells from
$\rho=\\{.761,\infty\\}$ which “mix” in the bistable region.

(a) k=10

(b) k=20

Figure 6: Vary $\rho$ for fixed k = 10 and k = 20

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 $\rho$
increases. As we continuously increase $\rho$ 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.

Figure 7: Vary $\rho$ for fixed k = .13 Figure 8: Cusp fold plot of $\rho$ and
k

Finally, in Figure 8 we plot the surface $B$. If we intersect the solution
plane at a fixed value of $k$ the solution curves are similar to Figures 6 and
7 depending on the choice of $k$. If we intersect the plane at the cusp point
where the surface folds we have a solution curve which looks like a pitchfork,
as demonstrated in Figure 5. Furthermore, if we project the fold lines of the
surface onto the $\rho,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.

Figure 9: Cusp fold plot analysis

###  2.4 Cusp Bifurcation in other networks

####  2.4.1 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.

(a) Plot varying $\alpha_{1},m$ where $\alpha_{2}=2$

(b) Plot varying $\alpha_{2},m$ where $\alpha_{1}=2$

Figure 10: MATCONT 2-gene cusp results

We used MATCONT to show that a cusp point occurs for the system

$\displaystyle\dot{x}$ $\displaystyle=\frac{\alpha_{1}}{1+y^{m}}-x$
$\displaystyle\dot{y}$ $\displaystyle=\frac{\alpha_{2}}{1+x^{k}}-y$

In Figure 10(a) we set $\alpha_{2}=2$ and $m=k>0$. As we varied $\alpha_{1}$
and $m$ we found a cusp exists at $(x,y,\alpha_{1},m)=(1.00,1.00,2.00,1.99)$.
In Figure 10(b) we set $\alpha_{1}=2$, $m=k>0$ and varied $\alpha_{1}$ and
$m$. A cusp point exists at $(x,y,\alpha_{2},m)=(1.00,1.00,2.00,1.99)$. This
aligns with Smale and Rajapakse’s results [29]. They set
$\alpha_{1}=\alpha_{2}=2$, $m=k>0$ and proved that a pitchfork bifurcation
occurs at the point $m=2$. They note that for all $0\leq 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.

(a) Plot of $x,m$

(b) Plot of $y,m$

Figure 11: MATCONT pitchfork bifurcation

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 $(\alpha_{1},m)$ and $(\alpha_{2},m)$ respectively. The
surfaces both are folds which we would expect to see for a cusp bifurcation.
Figures 10(a) and (b) are the projection of the fold on the $(\alpha_{1},m)$
and $(\alpha_{2},m)$ plane. As for the metastasis-regulatory network, the cusp
point separates the parameter ranges for continuous versus discontinuous
transitions.

(a) 2-gene x, $\alpha_{1},m$ plot

(b) 2-gene y, $\alpha_{2},m$ plot

Figure 12: Matlab fimplicit3 plots

####  2.4.2 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. 1.

There are 3 equilibria, two sinks and a saddle.

  2. 2.

The forward orbit of any point tends to one of the equilibria.

  3. 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:

$\dot{x}=\frac{\alpha}{1+z^{m}}-x$ $\dot{y}=\frac{\alpha x^{m}}{1+x^{m}}-y$
$\dot{z}=\frac{\alpha}{1+y^{m}}-z$

To find the behavior in the bistable region we set $\alpha=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.

(a)

(b)

Figure 14: Matlab plot of eigenspace and solution curves in blue. Lines in
magenta are the eigenvectors of the saddle point which determine the
eigenspace.

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,\alpha_{1}=\alpha_{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.

##  3 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 $\alpha(t)$ in the space
of parameters $\rho$ and $k$ from Figure 1.

(a)

(b)

Figure 15: (a) Time-evolution exclusively through monostable regions. Every
point on the path is monostable.(b) Image of monostable cells. Pro-metastatic
cells in black, Anti-metastatic cells in white and intermediate monostable
cells in continuous transition region in grey.

The time-evolution of the parameters $\alpha(t)=(\rho,k)$ shows possible
transitions between stable cell types. Here, the green line indicates how
values of $\alpha(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 $\alpha(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.

(a)

(b)

Figure 16: (a) Transition between monostable regions traversing bistable mixed
region. (b) Image of Discontinuous transition. Monostable Pro-metastatic cells
in black, Monostable Anti-metastatic cells in white and cells in mixed
Bistable region in both white and black.

In Figure 11 the path of $\alpha(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^{\prime}(t_{0})$. Over time the point $x^{\prime}(t)$
splits into a new stable point $x^{\prime}_{1}(t)$ and a saddle
$x^{\prime}_{2}(t)$. In the bistable region there are three equilibria $x(t)$,
$x^{\prime}_{1}(t)$ and $x^{\prime}_{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^{\prime}_{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^{\prime}_{1}(t)$ as $x(t)$ tends to $x^{\prime}_{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.

Figure 17: 1-dimensional time phases of stable states

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 $\alpha(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.

##  4 Acknowledgments

This research was supported by the Rich Summer Internship program and City
College. Michael Shub’s research was partially supported by the Smale
Institute. GB acknowledges support by the National Institutes of Health, NIGMS
MIRA Program (R35GM122561) and by the Laufer Center for Physical and
Quantitative Biology.

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.

##  5 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 equations
$\dot{x}=V(x,\alpha)$ as above where $x\in\mathbb{R}^{n}$ and $\alpha$ is a
parameter in $\mathbb{R}^{j}$ and
$V:\mathbb{R}^{n}\times\mathbb{R}^{j}\to\mathbb{R}^{n}$. In our case $\alpha$
is a real parameter so $j=1$.

Bifurcation theory studies how the behavior of the solution of $V$ changes as
$\alpha$ 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,\alpha)=0$ defines the equilibria. Recall that in
general

$V(x,\alpha):\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}$

and we assume that when $V(x,\alpha)=0$

equation 2) $DV(x,\alpha)$ goes from $\mathbb{R}^{n+1}\to\mathbb{R}^{n}$ and
the rank of the derivative is $n$.

So the set of zeros is locally a curve near $(x,\alpha)$ by the implicit
function theorem.

We assume one of the eigenvalues of the derivative $D_{x}V(x,\alpha)$ is $0$,
so we add

equation 3) $Det[D_{x}V(x,\alpha)]=0$

and now we assume that

$(V(x,\alpha),Det[D_{x}V(x,\alpha)]):\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}\times\mathbb{R}$

has rank $n+1$.

So $x,\alpha$ is an isolated solution of equations 2) and 3), and the curve of
zeros is tangent to the $\alpha=$constant plane. A saddle of index $k$ and one
of index $k-1$ collide and annihilate or appear if the parameter runs in the
opposite direction.

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 $\rho$ and $k$

$V(x,\rho,k):\mathbb{R}^{n}\times\mathbb{R}^{2}\to\mathbb{R}^{n}$

Then

1′)$V(x,\rho,k)=0$

defines the set of equilibria and we assume

2′)the rank of $D(V(x,\rho,k))=n$

so the set of equilibria near $(x,\rho,k)$ is a 2-dimensional surface. On the
surface we have those equilibria with a zero eigenvalue where the number of
equilibria might be changing, these are described by adding the equation

3′) $Det[D_{x}V(x,\rho,k)]=0$

and we assume that the map $(x,\rho,k)\to(V(x,\rho,k),Det[D_{x}V(x,\rho,k)])$
has rank $n+1$ where the two equations are satisfied. Hence the solution is a
curve near $(x,\rho,k)$.

When we project the curve on the parameter plane, $(\rho,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.

Figure 20: Fold projected on parameter plane

In case A) there are two possibilities:

Aa) 0 is a simple eigenvalue of $D_{x}(V(x,\rho,k)$ along the curve.In this
case the projected curve divides the plane near $(\rho,k)$.On one side there
is no equilibrium on the other side, there are two of index $j$ and $j-1$. So
in a region of interest to us, a saddle and a sink. If a curve in the
parameter plane crosses the projected curve we see a saddle-node bifurcation.

This is called a fold bifurcation.

Ab)There is a point along the curve where the rank of $D_{x}V(x,\rho,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 curve is the null space of the
derivative of $(V(x,\rho,k),Det[D_{x}V(x,\rho,k)])$. Now for a vertical vector
$(v,0,0)$ to be in the null space it is necessary and sufficient that

4′) $D_{x}(V(x,\rho,k)(v)=0$ and

5′) $\nabla_{x}Det[D_{x}V(x,\rho,k)]\bullet(v)=0$

The references above give methods for determining that these equations have a
solution. When equations 1′) and 3′) are satisfied $D_{x}V(x,\rho,k)$ has a
kernel as $Det[D_{x}V(x,\rho,k)]=0$ so we would like to express the vector $v$
in terms of $(x,\rho,k)$. Assuming that $D_{x}V(x,\rho,k)$ has rank $n-1$ when
$V(x,\rho,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,\rho,k)$.

So a cusp point is defined as a non-degenerate solution of the $n+2$ by $n+2$
system of equations

$\begin{split}V(x,\rho,k)&=0\\\ Det[D_{x}V(x,\rho,k)]&=0\\\
\nabla_{x}Det[D_{x}V(x,\rho,k)]\bullet(v)&=0\end{split}$ (2)

Where $v$ is a non-zero column of the adjugate matrix of $D_{x}V(x,\rho,k)$.

Here are some pictures of the cusp bifurcation when $n=1$.

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 saddle-node 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,\alpha)=0$ and $Det[D_{x}V(x,\alpha)]=0$ for
saddle node bifurcations and of the system of equations (2) where
$D_{x}V(x,\alpha)(v)=0$, $Det[D_{x}V(x,\alpha)]\bullet(v)=0$ and $v$ is a
column of the adjugate matrix of $D_{x}V(X,\alpha)$ for cusp bifurcations. The
region in Figure 1 looks like there might be a cusp bifurcation which puts a
sharper definition of the regions and supports the cell transitions in Figure
2.

##  6 Appendix B: Matlab Code

Figure 1

    %function zB = funB(t,B)
    function zB = funB(B)
    %no shBACH1 here

    global s S c K m b r rho a

    zB = s+(S-s)*K^b/(K^b+B^b) - B*( 1+ (a/(m^r*rho^r*(1+B)^r+1)) * (c/(1+c*B)) );

    function [Ld,Lt] = LOSSderiv(B)

    global s S c K m b r rho a

    % loss term
    Lt = B.*(1+ (a./(m^r*rho^r*(1+B).^r+1)).*(c./(1+c*B)));

    fderiv_num1=(m^r*rho^r*(1+B).^r+1).^2.*(1+c*B).^2;
    fderiv_num2=a*c*(m^r*rho^r*(1+B).^r+1).*(1+c*B);
    fderiv_num3=a*c^2*B.*(m^r*rho^r*(1+B).^r+1);
    fderiv_num4=a*c*m^r*rho^r*r*B.*(1+B).^(r-1).*(1+c*B);
    fderiv_den=(m^r*rho^r*(1+B).^r+1).^2.*(1+c*B).^2;

    % fd1=1+a*c./((m^r*rho^r*(1+B).^r+1).*(1+c*B));
    % fd2=B.*(a*c^2*(m^r*rho^r*(1+B).^r+1))./((m^r*rho^r*(1+B).^r+1).^2.*(1+c*B).^2);
    % fd3=B.*(a*c*m^r*rho^r*r*(1+B).^(r-1).*(1+c*B))./((m^r*rho^r*(1+B).^r+1).^2.*(1+c*B).^2);

    % derivative of loss term
    Ld=(fderiv_num1+fderiv_num2-fderiv_num3-fderiv_num4)./fderiv_den;
    %Ld=fd1-fd2-fd3;

    clear all; close all;

    global s S c K m b r rho a alpha

    L=3;
    paramrange=logspace(-L,L,100);

    B=logspace(-5,3,100);

    colorsetR=[0.5 0 0; 1 0 0; 1 .5 .8];
    colorsetB=[0 0 0.5; 0 0 1; 0 1 1];

    clear isbist;clear lowBACH1;
    s=0.02;
    S=20;
    c=200;
    m=2;
    b=3;
    r=5;
    a=1000;
    % scan rho: RKIP degradation rate
    % scan alpha: shBACH1 effect
    Krange=paramrange;
    rhorange=paramrange;
    ctr1=1;
    %figure;
    for rho=rhorange
        ctr2=1;
        for K=Krange
            lhs=B.*(1+ (a./(m^r*rho^r*(1+B).^r+1)).*(c./(1+c*B)));
            rhs=s+(S-s)*K^b./(K^b+B.^b);
            %subplot(7,7,(ctr1-1)*7+ctr2);plot(B,lhs,’LineWidth’,2);hold on;plot(B,rhs,’r’,’LineWidth’,2);set(gca,’XScale’,’Log’,’YScale’,’Log’);

            [Ld,Lt]=LOSSderiv(B);   % get the derivative of the loss term, and the loss term

            seekdownslope=1-sign(Ld);   %TRUE is loss term’s derivative is negative or 0
            midP=mean(B(find(seekdownslope)));
            %[a K any(seekdownslope) 3]

            %chance for bistable region next
            if(any(seekdownslope) && Ld(floor(mean(find(seekdownslope))))*Ld(end)<0 && Ld(floor(mean(find(seekdownslope))))*Ld(1)<0)

                zD1=fzero(’LOSSderiv’,[B(1) midP]);
                zD2=fzero(’LOSSderiv’,[midP B(end)]);

                [Ld1,Lt1]=LOSSderiv(zD1);
                [Ld2,Lt2]=LOSSderiv(zD2);
                Gt1=(s+(S-s)*K^b./(K^b+zD1.^b));    % gain term 1
                Gt2=(s+(S-s)*K^b./(K^b+zD2.^b));    % gain term 2
                isbist(ctr1,ctr2)=(Lt1-Gt1>0)*(Gt2-Lt2>0)*(zD2-zD1>0);
                lowBACH1(ctr1,ctr2)=(Gt2-Lt2<0);
                %[a K (Lt1-Gt1) (Gt2-Lt2)]
                %sBACH1a(ctr1,ctr2)=Gt1-Lt1;  % gain term - loss term = dB/dt
                %sBACH1b(ctr1,ctr2)=Gt2-Lt2;  % gain term - loss term = dB/dt
            %monostability is certain
            else
                isbist(ctr1,ctr2)=0;
                zD=fzero(’funB’,[0 B(end)*100]);
                lowBACH1(ctr1,ctr2)=(zD<0.2);
                %sBACH1(ctr1,ctr2)=zD<0.2;
            end;
            ctr2=ctr2+1;
        end;
        ctr1=ctr1+1;
    end;

    xvec=log10(rhorange);yvec=log10(Krange(1:47));
    [xx,yy]=ndgrid(xvec,yvec);
    % figure;mesh(xvec,yvec,atan2(yy,xx));
    atval=atan2(yy,xx);
    atval(48,length(Krange))=0;

    % figure;imagesc(log10(rhorange),log10(Krange),atval/min(min(atval)));colorbar;
    % figure;imagesc(log10(rhorange),log10(Krange),isbist);colorbar;

    B = bwboundaries(isbist);
    %mask=isbist>0.5;

    bistsurf=isbist*0.5+lowBACH1;
    %bistsurf=lowBACH1;
    %bistsurf(:,1:47)=(atval(:,1:47)/min(min(atval))).^0.75;
    bistsurf(:,1:47)=(atval(:,1:47)/min(min(atval)));

    figure;
    I=imagesc(log10(rhorange),log10(Krange),bistsurf);
    patch(log10(Krange(B{:}(:,2))),log10(rhorange(B{:}(:,1))),[0.8 0.7 0.2])
    %imagesc(log10(rhorange),log10(Krange),sBACH1);colorbar;
    hold on;cb=colorbar;
    %patch(log10(rhorange(isbist)),log10(Krange(isbist)),’g’);

    % plot(log10(10),log10(.05),’o’,’Color’,colorsetB(1,:),’LineWidth’,2,’MarkerSize’,20);
    % plot(log10(10),log10(.5),’o’,’Color’,colorsetB(2,:),’LineWidth’,2,’MarkerSize’,20);
    % plot(log10(10),log10(5),’o’,’Color’,colorsetB(3,:),’LineWidth’,2,’MarkerSize’,20);
    %
    % plot(log10(10),log10(.5),’^’,’Color’,colorsetR(3,:),’LineWidth’,2,’MarkerSize’,16);
    % plot(log10(.1),log10(.5),’^’,’Color’,colorsetR(1,:),’LineWidth’,2,’MarkerSize’,16);
    % plot(log10(.1),log10(5),’^’,’Color’,colorsetB(3,:),’LineWidth’,2,’MarkerSize’,16);

    plot(log10((Krange(47)+Krange(48))/2),log10((rhorange(50)+rhorange(51))/2),’*’,’Color’,’r’,’LineWidth’,2,’MarkerSize’,20);
    [(Krange(47)+Krange(48))/2 (rhorange(50)+rhorange(51))/2]

    set(gca,’FontSize’,16,’YDir’,’normal’);
    xlabel(’log_{10}(K): BACH1 insensitivity to self-repression’,’FontSize’,16);
    ylabel(’log_{10}(\rho): RKIP instability’,’FontSize’,16);
    colormap(gray);grid on;caxis([0 1]);
    text(0,2,{’Monostable’,’Pro-metastatic’},’Color’,[1 1 1],’FontSize’,18,’FontWeight’,’Bold’);
    text(0.1,-2.2,{’Monostable’,’Anti-metastatic’},’Color’,[0 0 0],’FontSize’,18,’FontWeight’,’Bold’);
    text(1.5,-.5,{’Bistable’,’Mixed’},’Color’,[0 0 0],’FontSize’,18,’FontWeight’,’Bold’);
    set(cb,’TickLabels’,flipud(get(cb,’TickLabels’)))

Figure 3

    function [lastDW,lastN,N,Wx,DWd,DWw]=NewtonsMethodfcn(eqns,x0,vars,count)

    NTable=[];
    WxTable=[];
    vars0=x0;
    W=double(subs(eqns,vars,vars0));
    DW=jacobian(eqns,vars);
    DWsub=double(subs(DW,vars,vars0));

    if rank(DW)==size(DW,2)
        DWdagger=inv(DWsub);
    else
        DWdagger=DWsub’*inv(DWsub*DWsub’);
    end
    DWW=DWdagger*W;
    Nprime=vars0’-DWW;
    WNprime=subs(eqns,vars,Nprime’);
    NTable=[vars0’,Nprime];
    WxTable=[W,WNprime];
    for j=1:count
        W=double(WNprime);
        DWsub=double(subs(DW,vars,Nprime’));
        if rank(DW)==size(DW)
            DWdagger=inv(DWsub);
        else
            DWdagger=DWsub’*inv(DWsub*DWsub’);
        end
        DWW=DWdagger*W;
        Nprime=Nprime-DWW;
        NTable=[NTable,Nprime];
        WNprime=double(subs(eqns,vars,Nprime’));
        WxTable=[WxTable,WNprime];
    end
    lastDW=DWsub;
    lastN=Nprime;
    N = double(NTable);
    Wx = double(WxTable);
    DWd=DWdagger;
    DWw=DWW;

    %script to generate table using Newtons Method function
    s=.02;
    S=20;
    c=200;
    m=2;
    b=3;
    r=5;
    a=1000;

    syms R L B rho K

    assume([R L B], ’Real’)
    assume([R L B], ’Positive’)

    eqn1=(1/(1+B))-rho*R;
    eqn2=((a*R^r)/(m^r+R^r))-L-c*L*B;
    eqn3=s+(((S-s)*K^b)/(K^b+B^b))-B-c*L*B;

    %a11*(a22*a33 - a23*a32) + a13*a21*a32 from Jacobian matrix
    eqn4=(-rho)*((-B*c-1)*(-L*c-(B^(b-1)*K^b*b*(S-s))/(B^b+K^b)^2-1)
        -(-L*c)*(-B*c))+(-1/(B+1)^2)*((R^(r-1)*a*r)/(R^r+m^r)-(R^r*R^(r-1)
        *a*r)/(R^r+m^r)^2)*(-B*c);
    eqns=[eqn1;eqn2;eqn3;eqn4];
    vars=[R,L,B,rho,K];

    %table format
    SizeW=[4 6];
    varTypesW={’double’,’double’,’double’,’double’,’double’,’double’};
    varsW={’W(X0)’,’W(X1)’,’W(X2)’,’W(X3)’,’W(X4)’,’W(X5)’};
    RowsW={’eqR’,’eqL’,’eqB’,’eqDet’};

    vars1=[0.932146,2.218409,0.043501,1.028071,0.134353]; %cusppt from matcont
    [newDW,newN,N,Wx,DWd,DWw]=NewtonsMethodfcn(eqns,vars1,vars,5);
    insertT=Wx(:,1:6);
    WX=table(’Size’,SizeW,’VariableTypes’,varTypesW,’VariableNames’,
        varsW,’RowNames’,RowsW);
    WX(:,:)=array2table(insertT)

Figure 4

    s=.02;
    S=20;
    c=200;
    m=2;
    b=3;
    r=5;
    a=1000;
    syms R L B rho K
    assume([R L B], ’Real’)
    assume([R L B], ’Positive’)
    eqn1=(1/(1+B))-rho*R;
    eqn2=((a*R^r)/(m^r+R^r))-L-c*L*B;
    eqn3=s+(((S-s)*K^b)/(K^b+B^b))-B-c*L*B;
    eqn4=(-rho)*((-B*c-1)*(-L*c-(B^(b-1)*K^b*b*(S-s))/(B^b+K^b)^2-1)
        -(-L*c)*(-B*c))+(-1/(B+1)^2)*((R^(r-1)*a*r)/(R^r+m^r)-(R^r*R^(r-
        1)*a*r)/(R^r + m^r)^2)*(-B*c);
    eqns=[eqn1;eqn2;eqn3;eqn4];
    vars=[R,L,B,rho,K];
    rhoKFW=[]; %list to save curve coordinates in forward direction
    rhoKBW=[]; %list to save curve coordinates in backward direction
    vars1=[0.932146,2.218409,0.043501,1.028071,0.134353]; %cusppt from matcont
    [newDW,newN,N,Wx,DWd,DWw]=NewtonsMethodfcn(eqns,vars1,vars,5);
    rhoKFW=[N(4:5,end)];
    rhoKBW=[N(4:5,end)];

    %adjust index range to plot more points
    for k=1:5
            deg=10^(-k);
            ker=null(newDW)*(deg);
            null(newDW);
            x0p=newN+ker;
            x0m=newN-ker;
            x1p=x0p;
            x1m=x0m;

            for j=1:20
                [newDW2,newN2,N2,Wx2,DWd2,DWw2]=
                    NewtonsMethodfcn(eqns,x1p’,vars,1);
                [newDW3,newN3,N3,Wx3,DWd3,DWw3]=
                    NewtonsMethodfcn(eqns,x1m’,vars,1);

                ker2=null(newDW2)*deg;
                x1p=newN2+ker2;
                rhoKFW=[rhoKFW,N2(4:5,end)];

                ker3=null(newDW3)*deg;
                x1m=newN3-ker3;
                rhoKBW=[rhoKBW,N3(4:5,end)];

            end
     end

    plot(rhoKFW(2,:),rhoKFW(1,:),’o’,rhoKBW(2,:),rhoKBW(1,:),’o’,
        0.134353,1.028071,’d’)
    xlabel(’K’)
    ylabel(’\rho’)

Figure 5

    s=.02;
    S=20;
    c=200;
    m=2;
    b=3;
    r=5;
    p=10;
    a=1000;
    syms R L B rho K
    assume([R L B], ’Real’)
    assume([R L B], ’Positive’)
    eqn1=(1/(1+B))-rho*R;
    eqn2=((a*R^r)/(m^r+R^r))-L-c*L*B;
    eqn3=s+(((S-s)*K^b)/(K^b+B^b))-B-c*L*B;
    %rho,K values from MATCONT CONTINUATION
    rho0=1.028071;
    K0=0.134353;
    D=jacobian([eqn1;eqn2;eqn3],[R,L,B,rho,K]);
    vars=[R,L,B,rho,K];

    %cusp pt from MATCONT
    varscusp=[0.932146,2.218409,0.043501,1.028071,0.134353];

    subs(D,vars,varscusp);
    null(subs(D,vars,varscusp));
    nullsp=null(double(subs(D,vars,varscusp)));
    nullsp(4,1);
    nullsp(5,1);

    vector=[nullsp(4,1)/nullsp(5,1);1];

    tvec=linspace(-.01,.01,100);
    Bvec=[];

    for i=1:5
        rhoK=vector*tvec(i)+[rho0;K0];
        X=vpasolve(subs([eqn1;eqn2;eqn3],[rho,K],[rhoK(1),rhoK(2)]));
        Bvec=[Bvec,X.B];
    end

    Bvec3=[];
    for i=51:100
        rhoK=vector*tvec(i)+[rho0;K0];
        X=vpasolve(subs([eqn1;eqn2;eqn3],[rho,K],[rhoK(1),rhoK(2)]));
        Bvec3=[Bvec3,X.B];
    end

    T=[tvec(1:50),tvec(51:100),tvec(51:100),tvec(51:100)];
    Bt=[Bvec(1:50),Bvec3(1,:),Bvec3(2,:),Bvec3(3,:)];

    figure(2)
    plot(T,Bt,’o’)
    xlabel(’T’)
    ylabel(’B’)

Figure 8

    s=.02;
    S=20;
    c=200;
    m=2;
    b=3;
    r=5;
    a=1000;

    syms R L B rho K
    assume([R L B rho K], ’Real’)
    assume([R L B rho K], ’Positive’)

    eqnB = @(K, rho, B) s + (S-s)*K^b/(K^b+B^b) - B*(1+a*c/((m^r*rho^r*(1+B)^r+1)*(1+c*B)));

    fimplicit3(eqnB)
    xlabel(’K’);
    ylabel(’rho’);
    zlabel(’B’);

    set(gca, ’Xscale’,’log’)
    set(gca, ’Yscale’,’log’)
    zlim([0,10])

Figure 12

    syms x y a m
    assume([x y a m],’Real’)
    assume([x y a m],’Positive’)

    figure
    fimplicit3(@(a,m,x) (a/(1+(2/(1+x^m))^m)) - x, [0 3 0 3 0 3])
    xlabel(’\alpha_1’); ylabel(’m’); zlabel(’x’);
    title(’2-gene network, plot x vary \alpha_1,m’)

    figure
    fimplicit3(@(a,m,y) (2/(1+(a/(1+y^m))^m)) - y, [0 3 0 3 0 3])
    xlabel(’\alpha_2’); ylabel(’m’); zlabel(’y’);
    title(’2-gene network, plot y vary \alpha_2,m’)

Figure 13

    alpha=2;
    x = -1:.2:3;
    y = -1:.2:3;
    [X,Y] = meshgrid(x,y);
    m = 10; %vary m
    V1 = alpha./(1+Y.^m)-X
    V2 = alpha./(1+X.^m)-Y
    uV1=V1./sqrt(V1.^2+V2.^2);
    uV2=V2./sqrt(V1.^2+V2.^2);
    hold on
    quiver(X,Y,uV1,uV2,.5,’c’)
    plot(1,1,"*r",2, 0,"ro",0, 2,"ro",x,x,’r’)
    str={’saddle’,{’sink’,’high CDK2’,’low Rb’},{’sink’,’low CDK2’,’high Rb’}}
    text([1.1 2.1 0.1],[1 0 2],str)
    text(-.5,-.5,’x=y’,’Color’,’red’)
    xlabel(’x’)
    ylabel(’y’)
    title(’2-gene vector field, m=10, \alpha=2’)
    hold off

Figures 14 - 16

    alpha=2;
    m=3;

    %% Solution curves and eigenspace

    figure(1)
    vec = linspace(0,4,3);
    plot3(1,1,1,’k*’)
    hold on

    syms x y z

    %eigenvectors of the saddle
    x0=[1,1,1] % saddle pt
    DxV=jacobian([
        (alpha/(1+z^m))-x;
        ((alpha*x^m)/(1+x^m))-y;
        (alpha/(1+y^m))-z],[x,y,z]);
    A=double(subs(DxV,[x,y,z,alpha,m],[x0,alpha,m]));
    [evec,eval]=eig(A);
    real1=(evec(:,1)+evec(:,2))/2;
    real2=A*real1;
    %eigenspace
    xs=linspace(0,2,5);
    ys=linspace(0,2,5);
    [Xs,Ys]=meshgrid(xs,ys);

    a=((real1/norm(real1))’+x0)-x0;
    b=((real2/norm(real2))’+x0)-x0;
    n=cross(a,b)
    Z=(-n(1).*(Xs-1)-n(2).*(Ys-1))./n(3) + 1;

    surf(Xs,Ys,Z,’FaceColor’,[1 0 0],’EdgeColor’,’none’,’FaceAlpha’,.4)
    eigdir=[((real1/norm(real1))’+x0);x0;((real2/norm(real2))’+x0)];
    plot3(eigdir(:,1),eigdir(:,2),eigdir(:,3),’m-’,’LineWidth’,2)

    clear x y z
    t=[0,20];
    for i=1:length(vec)
        for j=1:length(vec)
            for k=1:length(vec)
                X0=[vec(i),vec(j),vec(k)];
                [t,x_out]=ode45(@xdot, t, X0);
                plot3(x_out(:,1),x_out(:,2),x_out(:,3),’b:’)
                plot3(x_out(end,1),x_out(end,2),x_out(end,3),’k*’)
            end
        end
    end

    [t,x_out]=ode45(@xdot, [0 100], [0.5 0.5 2.5]);
    sink1=x_out(end,:);
    [t,x_out]=ode45(@xdot, [0 100], [2.5,2.5,.5]);
    sink2=x_out(end,:)-[0,0,.05];

    grid on
    xlim([0 2])
    ylim([0 2])
    zlim([0 2])
    xlabel(’x’)
    ylabel(’y’)
    zlabel(’z’)
    title(’Plot of solution curves divided by the eigenspace’)
    legend(’Equilibria’,’Eigenspace’)
    hold off

    %% Surfaces

    [x,y] = meshgrid(linspace(-1,10));

    yy =  (2*x.^m)./(x.^m+1);
    zz =  2./(x.^m+1);
    xx =  2./(x.^m+1);

    figure(2)
    plot3(1,1,1,’r*’,’MarkerSize’,12)
    hold on
    surf(xx,y,x,’FaceAlpha’,0.2,’FaceColor’,’b’,’LineStyle’,’none’)
    surf(x,yy,y,’FaceAlpha’,0.2,’FaceColor’,’r’,’LineStyle’,’none’)
    surf(x,y,zz,’FaceAlpha’,0.2,’FaceColor’,’g’,’LineStyle’,’none’)

    plot3(sink1(1),sink1(2),sink1(3),’ko’,’MarkerSize’,12)
    plot3(sink2(1),sink2(2),sink2(3),’ko’,’MarkerSize’,12)
    %surf(x,yy,xx)
    xlim([0 2])
    ylim([0 2])
    zlim([0 2])
    xlabel(’x’)
    ylabel(’y’)
    zlabel(’z’)
    grid on
    legend(’saddle’,’x’,’y’,’z’,’sink’)
    title(’3-gene solution surfaces and nullclines’)
    hold off

    %% Vector plot
    clear x y z
    figure(3)
    x = 0:.5:2.5;
    y = 0:.5:2.5;
    z = 0:.5:3.5;
    [X,Y,Z] = meshgrid(x,y,z);
    V1 = alpha./(1+Z.^m)-X;
    V2 = alpha*X.^m./(1+X.^m)-Y;
    V3 = alpha./(1+Y.^m)-Z;
    uV1=V1./sqrt(V1.^2+V2.^2+V3.^2);
    uV2=V2./sqrt(V1.^2+V2.^2+V3.^2);
    uV3=V3./sqrt(V1.^2+V2.^2+V3.^2);
    quiver3(X,Y,Z,uV1,uV2,uV3,.3)
    hold on
    plot3(1,1,1,’r*’)
    plot3(sink1(1),sink1(2),sink1(3),’ko’)
    plot3(sink2(1),sink2(2),sink2(3),’ko’)
    xlim([0 2])
    ylim([0 2])
    zlim([0 2])
    xlabel(’x’)
    ylabel(’y’)
    zlabel(’z’)
    title(’3-gene vector field, m=3, \alpha=2’)
    hold off

    %% Projections
    clear x y z X Y Z
    x = linspace(0,4,20);
    y = linspace(0,4,20);
    z = linspace(0,4,20);
    Y =  (2*x.^m)./(x.^m+1);
    Z =  2./(y.^m+1);
    X =  2./(z.^m+1);

    figure(4)
    % XY
    nexttile
    plot(x,Y)
    title(’Projection on x-y plane’)
    xlabel(’x’)
    ylabel(’y’)
    % ZX
    nexttile
    plot(X,z)
    title(’Projection on z-x plane’)
    xlabel(’x’)
    ylabel(’z’)
    % YZ
    nexttile
    plot(y,Z)
    title(’Projection on y-z plane’)
    xlabel(’y’)
    ylabel(’z’)

    %% ODE function

    function dx=xdot(t,x)
    alpha = 2;
    m = 3;
    dx=zeros(3,1);
    dx(1)=(alpha/(1+x(3)^m))-x(1);
    dx(2)=((alpha*x(1)^m)/(1+x(1)^m))-x(2);
    dx(3)=(alpha/(1+x(2)^m))-x(3);
    end

## References

  * [1] Thomas Brabletz.  Emt and met in metastasis: where are the cancer stem cells?  Cancer cell, 22(6):699–701, 2012.
  * [2] Gianluca Selvaggio et al.  Hybrid epithelial-mesenchymal phenotypes are controlled by microenvironmental factors.  Cancer research, 80(11):2407–2420, 2020.
  * [3] Hanah Goetz et al.  A plausible accelerating function of intermediate states in cancer metastasis.  PLoS computational biology, 16(3):e1007682, 10 Mar. 2020.
  * [4] Ievgenia Pastushenko et al.  Identification of the tumour transition states occurring during emt.  Nature, 556(7702):463–468, 2018.
  * [5] Jingyu Zhang et al.  Tgf-$\beta$-induced epithelial-to-mesenchymal transition proceeds through stepwise activation of multiple feedback loops.  Science signaling, 7(345):ra91, 30 Sep. 2014.
  * [6] Jiyoung Lee et al.  Network of mutually repressive metastasis regulators can promote cell heterogeneity and metastatic transition.  PNAS, vol. 111,3 (2014): E364-73. doi:10.1073/pnas.1304840111.
  * [7] Joseph B Addison et al.  Functional hierarchy and cooperation of emt master transcription factors in breast cancer metastasis.  Molecular cancer research : MCR, 19(5):784–798, 2021.
  * [8] Kamen P Simeonov et al.  Single-cell lineage tracing of metastatic cancer reveals selection of hybrid emt states.  Cancer cell, 39,8:1150–1162, (2021).
  * [9] Kari R Fischer et al.  Epithelial-to-mesenchymal transition is not required for lung metastasis but contributes to chemoresistance.  Nature, 527(7579):472–6, 2015.
  * [10] Keith S. Hoek et al.  In vivo switching of human melanoma cells between proliferative and invasive states.  Cancer Res, 68(3):650–6, 2008 Feb 1.
  * [11] Megan Mladinich et al.  Tackling cancer stem cells via inhibition of emt transcription factors.  Stem cells international, 2016:5285892, 2016.
  * [12] Mingyang Lu et al.  Microrna-based regulation of epithelial-hybrid-mesenchymal fate determination.  Proceedings of the National Academy of Sciences of the United States of America, 110(45):18144–9, 2013.
  * [13] Mohit Kumar Jolly et al.  Emt and met: necessary or permissive for metastasis?  Molecular oncology, 11(7):755–769, 2017.
  * [14] Nicole M Aiello et al.  Emt subtype influences epithelial plasticity and mode of cell migration.  Developmental cell, 45(6):681–695.e4, 2018.
  * [15] Sabrina L. Spencer et al.  The proliferation-quiescence decision is controlled by a bifurcation in cdk2 activity at mitotic exit.  Cell, 155(2):369–383, 2013.
  * [16] Xiaofeng Zheng et al.  Epithelial-to-mesenchymal transition is dispensable for metastasis but induces chemoresistance in pancreatic cancer.  Nature, 527,7579:525–530, (2015).
  * [17] Willy J.F Govaerts.  Numerical Methods for bifurcations of dynamical equilibrium.  SIAM, 2000.
  * [18] John Guckenheimer and Philip Holmes.  Non-linear Oscillations, Dynamical Systems and Bifurcation of Vector Fields.  Springer, 1983.
  * [19] John Guckenheimer and Yuri A. Kuznetsov.  Cusp bifurcation.  http://www.scholarpedia.org/article/Cusp_bifurcation, 2013.
  * [20] Pierre Gönczy and Yemima Budirahardja.  Coupling the cell cycle to development.  Development, vol. 136,17: 2861-72., 2009.
  * [21] Douglas Hanahan and Robert A Weinberg.  The hallmarks of cancer.  Cell, 100(1):57–70, 2000.
  * [22] Raghu Kalluri and Robert A. Weinberg.  The basics of epithelial-mesenchymal transition.  The Journal of clinical investigation, 119(6):1420–8, 2009.
  * [23] Marc W Kirschner and Li Victor C.  Molecular ties between the cell cycle and differentiation in embryonic stem cells.  Proceedings of the National Academy of Sciences of the United States of America, 111,26: 9503-8, 2014.
  * [24] Abraham Q Kohrman and David Q Matus.  Divide or conquer: Cell cycle regulation of invasive behavior.  Trends in Cell Biology., 27(1):12–25, 2017 Jan.
  * [25] Yuri Kuznetsov.  Elements of applied Bifurcation Theory, 3rd edition.  Springer, 2004.
  * [26] Yuri A. Kuznetsov.  Saddle-node bifurcation.  http://www.scholarpedia.org/article/Saddle-node_bifurcation, 2015\.
  * [27] Charles Pugh.  Mathematical Analysis.  Springer, 2015.
  * [28] Enrique R. Pujals, Michael Shub, and Yun-Jung Yang.  Stable and non-symmetric pitchfork bifurcations.  arXiv: Dynamical Systems, 2018.
  * [29] Indika Rajapakse and Steven Smale.  Mathematics of the genome.  arXiv, 2016.
  * [30] Indika Rajapakse and Steven Smale.  The pitchfork bifurcation.  International Journal of Bifurcation and Chaos, 27.09:1750132, 2017\.
  * [31] A. Dhooge W. Govaerts Yu.A. Kuznetsov H.G.E. Meijer B. Sautois.  New features of the software matcont for bifurcation analysis of dynamical systems.  MCMDS, Vol. 14, No. 2:pp 147–175, 2008.
  * [32] Jean Paul Thiery and Jonathan P Sleeman.  Complex networks orchestrate epithelial-mesenchymal transitions.  Nature reviews. Molecular cell biology, 7(2):131–42, 2006.