Now with the two nullclines, we can find all steady states by the intersections of the two curves. We will discuss the numerical methods to identify all steady states and the way to determine their stability.

Find intersections

We consider two curves \(A\) and \(B\), each of which is specified by a series of 2D vectors, i.e., \((\mathbf{A_1}, \mathbf{A_2}, ...,\mathbf{A_n}, \mathbf{A_{n+1}}, ...)\) for \(A\) and \((\mathbf{B_1}, \mathbf{B_2}, ...,\mathbf{B_n}, \mathbf{B_{n+1}}, ...)\) for \(B\). Here, \(\mathbf{A_n}\) represents the coordinates of the nth point along the curve \(A\), i.e., \((X(A_n), Y(A_n))\). The data frames null1_cont and null2_cont obtained in Part 3A are examples of the two curves. Now, we consider the four points \(\mathbf{A_n}\), \(\mathbf{A_{n+1}}\), \(\mathbf{B_n}\) and \(\mathbf{B_{n+1}}\), as shown in the figure below. The curve between two consecutive points can be approximated as a straight line. We will numerically evaluate (1) whether an intersection is formed by the two segments and if so (2) the coordinate of the intersection.

Figure 1
Figure 1


The figure shows two situations in which the two lines intersect (left) or not (right). Any point \(\mathbf{C}\) along the line (\(\mathbf{A_n}\), \(\mathbf{A_{n+1}}\)) can be written as

\[ \mathbf{C} = \mathbf{A_n} + \alpha (\mathbf{A_{n+1}} - \mathbf{A_n}) \tag{1} \] When \(0 \leq\alpha\leq 1\), \(\mathbf{C}\) is between \(\mathbf{A_n}\) and \(\mathbf{A_{n+1}}\). Otherwise, \(\mathbf{C}\) is outside of the two points.

If the point \(\mathbf{C}\) is also along the line (\(\mathbf{B_n}\), \(\mathbf{B_{n+1}}\)). We get

\[ \mathbf{C} = \mathbf{B_n} + \beta (\mathbf{B_{n+1}} - \mathbf{B_n}) \tag{2} \] From Equations (1) and (2), we have

\[ \alpha (\mathbf{A_{n+1}} - \mathbf{A_n}) - \beta (\mathbf{B_{n+1}} - \mathbf{B_n}) = \mathbf{B_n} - \mathbf{A_n} \tag{3} \] Equation (3) has two variables \(\alpha\) and \(\beta\) and two equations. Thus, \(\alpha\) and \(\beta\) can be determined as below.

For a 2-D vector \(\mathbf{V} = (V_X, V_Y)\), a perpendicular vector would be \(\mathbf{V_{\perp}} = (\frac{1}{V_X}, -\frac{1}{V_Y})\). That is because

\[\mathbf{V} \cdot \mathbf{V_{\perp}} = 0\] An easy way to solve \(\alpha\) in Equation (3) is to perform the dot product operation to the vector perpendicular to \(\mathbf{B_{n+1}} - \mathbf{B_n}\) to both sides of the equation. The vector, denoted as \(\mathbf{V_{B\perp}}\) here, can be expressed as

\[\mathbf{V_{B\perp}} = (\frac{1}{X_{B_{n+1}} - X_{B_n}}, -\frac{1}{Y_{B_{n+1}} - Y_{B_n}})\] With the dot product operation, we have

\[ \alpha (\frac{X_{A_{n+1}} - X_{A_n}}{X_{B_{n+1}} - X_{B_n}} - \frac{Y_{A_{n+1}} - Y_{A_n}}{Y_{B_{n+1}} - Y_{B_n}}) = (\frac{X_{B_n} - X_{A_n}}{X_{B_{n+1}} - X_{B_n}} - \frac{Y_{B_n} - Y_{A_n}}{Y_{B_{n+1}} - Y_{B_n}})\]

We then multiple both hand sides with \((X_{B_{n+1}} - X_{B_n})(Y_{B_{n+1}} - Y_{B_n})\) and get

\[\alpha = \frac{(X_{B_n} - X_{A_n})(Y_{B_{n+1}} - Y_{B_n}) - (Y_{B_n} - Y_{A_n})(X_{B_{n+1}} - X_{B_n}) }{(X_{A_{n+1}} - X_{A_n})(Y_{B_{n+1}} - Y_{B_n}) - (Y_{A_{n+1}} - Y_{A_n})(X_{B_{n+1}} - X_{B_n}) } \tag {4}\] Similarly, we get

\[\beta = \frac{(X_{B_n} - X_{A_n})(Y_{A_{n+1}} - Y_{A_n}) - (Y_{B_n} - Y_{A_n})(X_{A_{n+1}} - X_{A_n}) }{(X_{A_{n+1}} - X_{A_n})(Y_{B_{n+1}} - Y_{B_n}) - (Y_{A_{n+1}} - Y_{A_n})(X_{B_{n+1}} - X_{B_n}) } \tag {5}\] The intersection \(C\) can be then calculated from either Equation (1) or (2).

find_intersection <- function(lineA, lineB, nA, nB) {
  dAX = lineA[nA+1,1] - lineA[nA,1]
  dAY = lineA[nA+1,2] - lineA[nA,2]
  
  dBX = lineB[nB+1,1] - lineB[nB,1]
  dBY = lineB[nB+1,2] - lineB[nB,2]
  
  dABX = lineB[nB,1] - lineA[nA, 1]
  dABY = lineB[nB,2] - lineA[nA, 2]
  
  d = dAX * dBY - dAY * dBX
  
  alpha = (dABX * dBY - dABY*dBX)/d
  beta = (dABX * dAY - dABY*dAX)/d
  
  if((alpha*(1-alpha) >= 0) &
     (beta*(1-beta) >= 0)) {   # check whether there is an intersection in between
    intersection = c((1-alpha)*lineA[nA,1] + alpha *lineA[nA+1,1], 
                    (1-alpha)*lineA[nA,2] + alpha *lineA[nA+1,2],1)
    # 1st and 2nd elements: X & Y
    # 3rd element: 1 means found, 0 means not found
  }
  else{
    intersection = c(alpha,beta,0)
  }
  return(intersection)
}

# test example with three points in each line. They intersect by the first segments.
lineA = array(c(0,1,2,0,0,0), c(3,2)) 
lineB = array(c(0.5,0.5,0.5,1,-1,-3), c(3,2))

# Plotting lines in R
plot(lineA, type='l', col='red', lwd=2, pch=16, ylim=c(min(lineA[,2], lineB[,2]), max(lineA[,2], lineB[,2])),
     xlab='X', ylab='Y')
lines(lineB, type='l', col='blue', lwd=2, pch=16)
points(lineA, col='red', pch=16)
points(lineB, col='blue', pch=16)

# Adding labels and legend
legend('topright', legend=c('Line A', 'Line B'), col=c('red', 'blue'), lty=1, lwd=2, pch=16)

find_intersection(lineA, lineB, 1, 1) # Intersect at (0.5, 0)
## [1] 0.5 0.0 1.0
find_intersection(lineA, lineB, 1, 2) # No intersection
## [1]  0.5 -0.5  0.0
find_intersection(lineA, lineB, 2, 2) # No intersection
## [1] -0.5 -0.5  0.0

We can then find all intersections by running the find_intersection function through all possible line segments, as shown in the code below.

find_intersection_all <- function(lineA, lineB) {
  lineA_ind = seq_len(nrow(lineA)-1) 
  lineB_ind = seq_len(nrow(lineB)-1) 
  lines_all = expand.grid(lineA_ind, lineB_ind)   # all combinations
  results = apply(lines_all, MARGIN = 1, function(inds) {
    return(find_intersection(lineA, lineB, inds[1], inds[2]))})
  return(t(results[1:2,which(results[3,] == 1)])) # subsetting the matrix, keep intersections only 
}

load("./extra/data/03A/null_03A.Rdata")

ss_all = find_intersection_all(null1_cont, null2_cont)
ss_all
##              X         X
## [1,]  52.90118 361.57784
## [2,] 189.96888 126.64398
## [3,] 542.35734  35.27042

This approach works fine. However, it is very slow, as it has to go through all pairs of line segments. To address this issue, we can search along of the first nullcline \(\frac{dX}{dt}=0\), and check for sign flipping of \(\frac{dY}{dt}\) between two consecutive points. Similarly, we can search along the second nullcline \(\frac{dY}{dt}=0\), and check for sign flipping of \(\frac{dX}{dt}\) between two consecutive points. Once we find those cases, we then identify the intersection value for any of the combinations.

# the same derivative function as the one in Part 3A
hill_inh <- function(X,X0,n) {
  a = (X/X0)**n
  return(1/(1+a))
}

derivs_ts <- 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))
}

#dX/dt along the nullcline dY/dt = 0
dxdt_y_null = apply(null2_cont, MARGIN = 1, function(Xs) {
  d = derivs_ts(0, Xs)
  return(d[1])})
#dY/dt along the nullcline dX/dt = 0
dydt_x_null = apply(null1_cont, MARGIN = 1, function(Xs) {
  d = derivs_ts(0, Xs)
  return(d[2])})

find_intersection_all_fast <- function(lineA, lineB, dxdt, dydt) {
  small = 10^-3    # small is used to allow some tolerance of numerical errors in detecting sign flipping
  lineA_ind = which(dydt * c(dydt[-1],NA)<=small)  # find sign flipping between two consecutive points
  lineB_ind = which(dxdt * c(dxdt[-1],NA)<=small) 
  
  lines_all = expand.grid(lineA_ind, lineB_ind)   # all combinations of the above
  results = apply(lines_all, MARGIN = 1, function(inds) {
    return(find_intersection(lineA, lineB, inds[1], inds[2]))})
  return(t(results[1:2,which(results[3,] == 1)]))
}

ss_all = find_intersection_all_fast(null1_cont, null2_cont, dxdt_y_null, dydt_x_null)
ss_all
##              X         X
## [1,]  52.90118 361.57784
## [2,] 189.96888 126.64398
## [3,] 542.35734  35.27042

The above two methods give the same results. But the second method is way faster.

Stability

We will now explore the stability of the identified steady states.

Consider a system following

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

It is at a steady state \((X_s, Y_s)\) with small perturbations, \(X = X_s + \Delta X\) and \(Y = Y_s + \Delta Y\). By Taylor expansion and \(f_X(X_s, Y_s) = 0\), \(f_Y(X_s, Y_s) = 0\), we get

\[\begin{equation} \frac{d}{dt} \begin{pmatrix} \Delta X \\ \Delta Y \end {pmatrix}= \begin{pmatrix} \frac{\partial f_X}{\partial X} & \frac{\partial f_X}{\partial Y} \\ \frac{\partial f_Y}{\partial X} & \frac{\partial f_Y}{\partial Y} \end {pmatrix} \begin{pmatrix} \Delta X \\ \Delta Y \end {pmatrix}\tag{6} \end{equation}\]

The matrix

\[\begin{equation} \mathbf{J} = \begin{pmatrix} \frac{\partial f_X}{\partial X} & \frac{\partial f_X}{\partial Y} \\ \frac{\partial f_Y}{\partial X} & \frac{\partial f_Y}{\partial Y} \end {pmatrix} \end{equation}\]

is the Jacobian matrix that we have discussed previously.

We now set \(\Delta X(t) = \delta x e^{\lambda t}\), where \(\delta x\) and \(\lambda\) are two constants. Similarly, \(\Delta Y(t) = \delta y e^{\lambda t}\). Together with Equation (6), we have

\[\begin{equation} \lambda \begin{pmatrix} \delta x \\ \delta y \end {pmatrix}= \mathbf{J} \begin{pmatrix} \delta x \\ \delta y \end {pmatrix}\tag{7} \end{equation}\]

This is a typical eigenvalue problem. The solutions of \(\lambda\) are the eigenvalues. Now, we implement this calculation. Instead of using the analytical forms of the Jacobian matrix, we will estimate it numerically from the rate equations.

# A generic function to check stability for a 2D ODE sytem
stability_2D <- function(derivs, ss, ...) { # ss is a vector of steady state values X_S, Y_S
  delta = 0.001
  f_current = derivs(0,ss, ...)   # f(x,y) this is computed, just in case it is not exactly 0
  f_plus_dx = derivs(0,ss + c(delta,0), ...) # f(x+dx, y)
  f_plus_dy = derivs(0,ss + c(0, delta), ...) # f(x, y+dx)
  
  # finite difference to approximate Jacobian
  dfxdx = (f_plus_dx[1] - f_current[1])/delta
  dfxdy = (f_plus_dy[1] - f_current[1])/delta
  dfydx = (f_plus_dx[2] - f_current[2])/delta
  dfydy = (f_plus_dy[2] - f_current[2])/delta
  
  jacobian = array(c(dfxdx, dfydx, dfxdy, dfydy), c(2,2))
  lambda = eigen(jacobian)$values
  if(class(lambda[1]) == "complex") {
    if(Re(lambda[1]) < 0){
      stability = 4   # stable spiral
    }else{
      stability = 5   # unstable spiral
    }
  }
  else{
    if((lambda[1] < 0) & (lambda[2] <0)){
      stability = 1   # stable 
    }else if((lambda[1] >= 0) & (lambda[2] >= 0)){
      stability = 2   # unstable 
    }else{
      stability = 3   # saddle
    }
  }
  return(stability)
}

ss_with_stability =  t(apply(ss_all, MARGIN = 1, function(ss) {
                          return(c(ss,stability_2D(derivs = derivs_ts, ss = ss)))}))
colnames(ss_with_stability) = c("X", "Y", "Stability")
ss_with_stability
##              X         Y Stability
## [1,]  52.90118 361.57784         1
## [2,] 189.96888 126.64398         3
## [3,] 542.35734  35.27042         1

The stability of a steady state can be determined by the eigenvalues according to following table. Note that the eigenvalues of a Jacobian matrix are not necessarily real numbers. See Part 3G for a few other examples of gene circuits with different types of steady-state stability.

Conditions Stability
\(\lambda_1 < 0\) & \(\lambda_2 < 0\) Stable
\(\lambda_1 \geq 0\) & \(\lambda_2 \geq 0\) Unstable
\(\lambda_1 \lambda_2 < 0\) Saddle point
\(Re(\lambda_1) < 0\) & \(Re(\lambda_2) < 0\) Stable spiral
\(Re(\lambda_1) \geq 0\) & \(Re(\lambda_2) \geq 0\) Unstable spiral

In the current example, two steady states are stable and the other one is a saddle point.