In the previous section, we show typical dynamical analysis based on phase plane for two-variable non-linear dynamical systems. In the case of toggle switch gene circuit, the vector field and time trajectories indicates that there are two basins. A term named separatrix is a curve in the phase plane that defines the boundary of the basins. Starting from any initial condition from one side of the separatrix will always lead to the corresponding steady state (attractor). Usually it is non-trivial to obtain separatrix in general. But for certain cases like saddle-node systems, there is an easy trick to obtain the separatrix.

We again consider the following ODEs,

\[\begin{equation} \begin{cases} \frac{dX}{dt} = f_X(X,Y) \\ \frac{dY}{dt} = f_Y(X,Y) \end{cases} \end{equation}\]

But now we modify both equations by changing the sign only:

\[\begin{equation} \begin{cases} \frac{dX}{dt} = -f_X(X,Y) \\ \frac{dY}{dt} = -f_Y(X,Y) \end{cases} \end{equation}\]

We first plot the vector field of the modified ODEs in the phase plane.

library(ggplot2)
hill_inh <- function(X,Xth,n) {
  a = (X/Xth)**n
  return(1/(1+a))
}

derivs_ts_revised <- function(t, Xs) {
  X = Xs[1]
  Y = Xs[2]
  dxdt = -(5 + 50 * hill_inh(Y, 100, 4) - 0.1 * X)
  dydt = -(4 + 40 * hill_inh(X, 150, 4) - 0.12 * Y)
  return(c(dxdt, dydt))
}

X_all = seq(0, 600, by=20)   # all X grids
Y_all = seq(0, 650, by=20)   # all Y grids
XY_all = expand.grid(X_all, Y_all)   # all combinations of X and Y

results_unit = t(apply(XY_all, MARGIN=1, function(Xs) {
  v = derivs_ts_revised(0, Xs)
  v_norm = v / sqrt(v[1]**2 + v[2]**2)
  return(c(Xs, v_norm))})) # generate all vector field data
colnames(results_unit) = c("X", "Y", "dX", "dY")

p1 = ggplot(data=as.data.frame(results_unit), aes(x=X, y=Y)) + 
  geom_segment(aes(xend=X+20*dX, yend=Y+20*dY), arrow = arrow(length = unit(0.05,"in")))
p1

Now, the middle stead state remains a saddle point. But instead of being “repulsive”, the vectors point towards the middle steady state along the direction of the three steady states. Once reaching the steady state, the vectors point outwards along the boundary, which is the separatrix. Thus, we can sample points around the saddle point as the initial conditions for the modified ODEs to numerical obtain the separatrix.

# 4th order Runge-Kutta (RK4) for a generic multi-variable system
RK4_generic <- function(derivs, X0, t.total, dt, ...){
  # derivs: the function of the derivatives 
  # X0: initial condition, a vector of multiple variables
  # t.total: total simulation time, assuming t starts from 0 at the beginning
  # dt: time step size 
  t_all = seq(0, t.total, by=dt)
  n_all = length(t_all)
  nx = length(X0)
  X_all = matrix(0, nrow = n_all, ncol = nx)
  X_all[1,] =  X0
  for (i in 1:(n_all-1)) {
    t_0= t_all[i]
    t_0.5 = t_0 + 0.5*dt
    t_1 = t_0 + dt
    k1 = dt * derivs(t_0,X_all[i,],...)
    k2 = dt * derivs(t_0.5,X_all[i,] + k1/2,...)
    k3 = dt * derivs(t_0.5,X_all[i,] + k2/2,...)
    k4 = dt * derivs(t_1,X_all[i,] + k3,...)
    X_all[i+1,] = X_all[i,] + (k1+2*k2+2*k3+k4)/6
  }
  return(cbind(t_all, X_all))   # the output is a matrix of t & X(t) for all time steps
}

derivs_ts_orig <- function(t, Xs) {
  X = Xs[1]
  Y = Xs[2]
  dxdt = 5 + 50 * hill_inh(Y, 100, 4) - 0.1 * X
  dydt = 4 + 40 * hill_inh(X, 150, 4) - 0.12 * Y
  return(c(dxdt, dydt))
}

set.seed(77)   # set the seed for the random number generator, remove this line to generate more random initial conditions
t.total = 100
dt  = 0.01
X_range = c(0,600)
Y_range = c(0,650)

# Simulations of the original system
X_init_1 = array(runif(20, 0, 600), dim = c(10,2)) # generate 10 random initial conditions
for(i in seq_len(nrow(X_init_1))){
  results = RK4_generic(derivs_ts_orig, X_init_1[i,], t.total, dt)
  colnames(results) = c("t", "X", "Y")
  p1 = p1 + geom_path(data=as.data.frame(results), aes(x=X, y=Y), color = "blue")
}

# Simulation to find separatrix
X_init_2 = array(runif(20, -1, 1), dim = c(10,2)) # generate 10 random initial conditions
for(i in seq_len(nrow(X_init_2))){
  X_init = X_init_2[i,] + c(189.96888, 126.64398)
  results = RK4_generic(derivs_ts_revised, X_init, t.total, dt)
  results = results[which(((results[,2] - X_range[1])*(results[,2] - X_range[2]) <0) &
          ((results[,3] - Y_range[1])*(results[,3] - Y_range[2]) < 0)),] # filter points outside of the plotting area
  colnames(results) = c("t", "X", "Y")
  p1 = p1 + geom_path(data=as.data.frame(results), aes(x=X, y=Y), color = "red")
}

p1 + xlim(0, 600) + ylim(0, 650)
## Warning: Removed 124 rows containing missing values (geom_segment).

In the above plot, the trajectories representing the separatrix (red) are shown together with the vector field of the modified ODEs (black) and the trajectories of the original ODEs (blue). We can observe that the blue lines never go across the red line.