In Part 3, we will generalize what we learn from Part 2 to ODEs with two variables and apply them to model a gene circuit with two genes. We will focus on phase plane analysis for a visual characterization of nonlinear dynamical systems.
We first model the toggle switch circuit with two genes \(X\) and \(Y\), whose dynamics can be described by a set of two ODEs below (i.e., Equation (13) in Part 2A).
\[\begin{equation} \begin{cases} \frac{dX}{dt} = g_{X0} + g_{X1}\frac{1}{1+(Y/Y_{th})^{n_Y}} - k_XX \\ \frac{dY}{dt} = g_{Y0} + g_{Y1}\frac{1}{1+(X/X_{th})^{n_X}} - k_YY \end{cases} \end{equation}\]
Below shows functions to compute the derivatives of the toggle switch circuit. For illustration, the kinetic parameters of the circuit are hard coded in the script.
hill_inh <- function(X,X_th,n) { # inhibitory Hill function
# X_th: Hill threshold, n: Hill coefficient
a = (X/X_th)**n
return(1/(1+a))
}
derivs_ts <- function(t, Xs) { # Calculate derivative function for a toggle switch circuit
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))
}
We can simulate the dynamics of \(X\) and \(Y\) with 4th order Runge-Kutta (RK4). The function RK4_generic is a generic version of RK4 for a multi-variable system, very similar to the function RK4 in Part 2B. Now, the time trajectory is stored in a matrix called X_all.
# 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
}
set.seed(77) # set the seed for the random number generator, remove this line to generate more random initial conditions
X_init_all = array(runif(20, 0, 600), dim = c(10,2)) # generate 10 random initial conditions
t.total = 100
dt = 0.01
plot(NULL, xlab="t (Minute)", ylab="Levels (nM)", xlim=c(0,80), ylim=c(0,650))
for(i in seq_len(nrow(X_init_all))){
results = RK4_generic(derivs = derivs_ts, X0 = X_init_all[i,], t.total = t.total, dt = dt)
lines(results[,1], results[,2], col = 1)
lines(results[,1], results[,3], col = 2)
}
legend("top", inset=0.02, legend = c("X", "Y"),
col=1:2, lty=1, cex=0.8)
We randomly selected 10 initial conditions and simulate the ODEs starting from each. As shown in the above plot, the system can reach to two stable steady states, i.e., two sets of \(X\) and \(Y\).
A powerful way to characterize the dynamic behavior of such a two-variable system is the phase plane analysis. A phase plane describes the states of a system by two axes representing, e.g. in this case, the levels of the two variables. Here, we plot the time trajectories from the previous simulations to the phase plane of \(X\) and \(Y\).
plot(NULL, xlab="X (nM)", ylab="Y (nM)", xlim=c(0,600), ylim=c(0,650))
for(i in seq_len(nrow(X_init_all))){
results = RK4_generic(derivs = derivs_ts, X0 = X_init_all[i,], t.total = t.total, dt = dt)
lines(results[,2], results[,3], col = 1)
}
The trajectories shown in the phase plane illustrate again the existence of two stable steady states.
We can also plot arrows for points in the phase plane to represent the direction of changes.
\[\begin{equation} \begin{cases} V_X = \frac{dX}{dt} \\ V_Y = \frac{dY}{dt} \end{cases} \end{equation}\]
library(ggplot2)
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 = t(apply(XY_all, MARGIN=1, function(Xs) {return(c(Xs, derivs_ts(0, Xs)))})) # generate all vector field data
colnames(results) = c("X", "Y", "dX", "dY")
ggplot(data=as.data.frame(results), aes(x=X, y=Y)) +
geom_segment(aes(xend=X+dX, yend=Y+dY), arrow = arrow(length = unit(0.05,"in")))
In the vector filed, vectors are pointing towards two stable steady states. There is also another steady state in the middle of the plot, where the nearby vectors point toward it from one direction and point outwards from another direction. This steady state is a saddle point, a type of unstable steady state.
We can choose to plot the unit vectors of the vector field to highlight the directions.
results_unit = t(apply(XY_all, MARGIN=1, function(Xs) {
v = derivs_ts(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
Another convenient way to evaluate the behavior of a nonlinear dynamical system is the techniques of nullcline analysis. A nullcline corresponds to the steady state condition of one ODE. Thus, there are two nullclines for (1) \(\frac{dX}{dt} = 0\) and (2) \(\frac{dY}{dt} = 0\). In the current example, the first nullcline satisfies
\[f_X(X, Y) = g_{X0} + g_{X1}\frac{1}{1+(Y/Y_{th})^{n_Y}} - k_XX = 0 \tag{1} \] The second one satisfies
\[ f_Y(X, Y) = g_{Y0} + g_{Y1}\frac{1}{1+(X/X_{th})^{n_X}} - k_YY = 0 \tag{2} \]
The curve for an equation like \(f_X(X, Y) = 0\) can be obtain numerically, as we will discussed very soon. But for this relatively simple example, we can get the curve right away by separation of variables.
From Equation (1), we can move \(k_XX\) to the right hand side, and we have
\[X(Y) = \frac{g_{X0} + g_{X1}\frac{1}{1+(Y/Y_{th})^{n_Y}}}{k_X} \tag{3}\]
Similarly, from Equation (2), we can also get
\[ Y(X) = \frac{g_{Y0} + g_{Y1}\frac{1}{1+(X/X_{th})^{n_X}}}{k_Y} \tag{4} \]
The first nullcline (Equation (3)) can also be expression as \(Y(X)\). Although we don’t have to, as \(X(Y)\) is sufficient to characterize the nullcline curve numerically.
X_all = seq(0, 600, by=1) # all X grids
Y_all = seq(0, 650, by=1) # all Y grids
nullcline1 <- function(Y) {
X= (5 + 50 * hill_inh(Y, 100, 4))/0.1
return(cbind(X, Y))
}
nullcline2 <- function(X) {
Y = (4 + 40 * hill_inh(X, 150, 4))/0.12
return(cbind(X, Y))
}
null1 = nullcline1(Y_all)
null2 = nullcline2(X_all)
p1 +
geom_path(data = as.data.frame(null1), aes(x=X, y=Y, colour = "dX/dt=0"), size=1) +
geom_path(data = as.data.frame(null2), aes(x=X, y=Y, colour = "dY/dt=0"), size=1)
The intersections of the two nullclines are the steady states of the system. By checking the number of intersections, we can identify the multi-stability property easily. The separation of variables should be applied whenever possible. However, it may not be feasible in many situations. We then have to rely on other approaches, such as contour method and numerical continuation.
One of such approaches relies on the contour analysis. We consider \(Z = f_X(X,Y)\) and construct a landscape of a plane of \(X\) and \(Y\), with \(Z\) as the height. The nullcline would be the contour line when \(Z=0\). This can be easily done with the contour function in base R or other plotting packages, such as geom_contour from ggplot2. In this approach, \(Z\) values for all 2D grid points are required to compute the contour. This can be computational intensive when high resolution grid points are needed. The contour method is usually robust, as long as the \(Z\) values are smooth enough for interpolation.
X_all = seq(0, 600, by=10) # all X grids
Y_all = seq(0, 650, by=10) # all Y grids
nX = length(X_all)
nY = length(Y_all)
XY_all = expand.grid(X_all, Y_all) # all combinations of X and Y
results = t(apply(XY_all, MARGIN=1, function(Xs) {return(c(Xs, derivs_ts(0, Xs)))})) # generate all vector field data
z_X = array(results[,3], dim = c(nX, nY))
z_Y = array(results[,4], dim = c(nX, nY))
contour(X_all, Y_all, z_X, levels = 0, col = 2, xlab = "X", ylab = "Y", drawlabels = F)
contour(X_all, Y_all, z_Y, levels = 0, col = 3, add = TRUE, drawlabels = F)
Another approach is the numerical continuation, which we have already applied in the bifurcation analysis in Part 2E. In that previous analysis, we identify the bifurcation curve that satisfies \(f(X, k) = 0\). Similarly, we will identify the nullcline curve that satisfies \(f_X(X, Y) = 0\) (and the other one that obeys \(f_Y(X, Y) = 0\)). Here, we can choose \(Y\) as the control parameter and find \(X(Y)\) (that’s what we will do below). Of course, it is also fine to choose \(X\) as the control parameter and find \(Y(X)\).
f_null1 <- function(X, Y) { # f_X(X, Y)
return(5 + 50 * hill_inh(Y, 100, 4) - 0.1 * X)
}
dfdX_null1 <- function(X, Y) { # Jacobian (1st order partial derivatives) for f_X(X, Y)
return(-0.1)
}
dfdY_null1 <- function(X, Y) {
y_frac = (Y/100)**4
return(-200/Y*y_frac/(1+y_frac)**2)
}
f_null2 <- function(X, Y) { # f_Y(X, Y)
return(4 + 40 * hill_inh(X, 150, 4) - 0.12 * Y)
}
dfdX_null2 <- function(X, Y) { # Jacobian for f_Y(X, Y)
x_frac = (X/150)**4
return(-160/X*x_frac/(1+x_frac)**2)
}
dfdY_null2 <- function(X, Y) {
return(-0.12)
}
# Newton's Method (single variable)
find_root_Newton <- function(X, func, dfunc, X_range, error = 10^-3, ...) {
#X: Initial guess of X
#func: function f(X,...)
#dfunc: df/dX
#X_range: lower and upper limits of root X. If X is outside of the range, the algorithm stops.
f = func(X, ...)
while(abs(f) > error){
X = X - f/dfunc(X, ...)
if((X-X_range[1])*(X-X_range[2]) > 0) break # Check if X is in within X_range;
# This would avoid potential infinite loop; When this occurs, the Newton's method doesn't converge.
f = func(X, ...)
}
return(X)
}
# 4th order Runge-Kutta (RK4) for 1D
# This function issimplified to
# (1) only output steady-state X or final X,
# (2) derivative is not the function of t in this implementation, thus some t-dependent parts are simplified
RK4_1D_steady <- function(derivs, X0, t.total, dt, ...){
# derivs: the function of the derivatives
# X0: initial condition
# t.total: total simulation time, assuming t starts from 0 at the beginning
# dt: time step size
error = 10^-8
t = 0
X = X0
while(t <= t.total){
df = derivs(X, ...)
if(abs(df) < error) break
k1 = dt * df
k2 = dt * derivs(X + k1/2, ...)
k3 = dt * derivs(X + k2/2, ...)
k4 = dt * derivs(X + k3, ...)
X = X + (k1+2*k2+2*k3+k4)/6
t = t + dt
}
return(X)
}
nullcline_numcont <- function(func, dfdX, dfdY, X_range, Y_range, X_init, Y_init, direct = 1, step_arc) {
nmax_cycle = 10000
results = matrix(NA, nrow = nmax_cycle, ncol = 2)
cycle = 1
step_Y_previous = direct
step_X_previous = 0
X_new = X_init
Y = Y_init
results[cycle,] = c(X_new, Y)
while(((X_new - X_range[1])*(X_new - X_range[2]) <= 0) & # boundary check
((Y - Y_range[1]) * (Y - Y_range[2]) <= 0) &
(cycle < nmax_cycle)) {
h = -dfdY(X_new, Y)/dfdX(X_new, Y)
step_Y = step_arc/sqrt(1+h**2)
step_X = step_Y*h
if((step_Y_previous * step_Y + step_X_previous * step_X) < 0){
step_Y = - step_Y
step_X = - step_X
}
step_Y_previous = step_Y
step_X_previous = step_X
Y = Y + step_Y
X_init = X_new + step_X
X_new = find_root_Newton(X=X_init, func=func, dfunc=dfdX, X_range = X_range, Y = Y)
cycle = cycle + 1
results[cycle,] = c(X_new, Y)
}
results = na.omit(results)
colnames(results) = c("X","Y")
return(results)
}
X_range = c(0, 600)
Y_range = c(0, 650)
# the initial X,Y need to be provided, here we manually selected them with an initial simulation
X_init_guess_null1 = 500
Y_init_guess_null1 = 1
X_init_null1 = RK4_1D_steady(derivs = f_null1, X0 = X_init_guess_null1, Y = Y_init_guess_null1,
t.total = 1000, dt = 0.01) # start from ~ (550, 1)
null1_cont = nullcline_numcont(f_null1, dfdX_null1, dfdY_null1, X_range, Y_range,
X_init = X_init_null1, Y_init = Y_init_guess_null1, direct = 1, step_arc = 3)
X_init_guess_null2 = 1
Y_init_guess_null2 =366.6
X_init_null2 = RK4_1D_steady(derivs = f_null2, X0 = X_init_guess_null2, Y = Y_init_guess_null2,
t.total = 1000, dt = 0.01) # start from ~ (8,366.6)
null2_cont = nullcline_numcont(f_null2, dfdX_null2, dfdY_null2, X_range, Y_range,
X_init = X_init_null2, Y_init = Y_init_guess_null2, direct = -1, step_arc = 3)
plot(NULL, xlab="X", ylab="Y", xlim=X_range, ylim=Y_range)
points(null1_cont[,1], null1_cont[,2], pch = 1, cex = 0.3, col = 2)
points(null2_cont[,1], null2_cont[,2], pch = 1, cex = 0.3, col = 3)
legend("topright", inset=0.02, legend = c("dX/dt=0", "dY/dt=0"), col=2:3, lty=1, cex=0.8)
The core code is equivalent to the code for 1-variable bifurcation analysis in Part 2E. It looks lengthy and somewhat cumbersome, as the boundary of two variables and the initial direction of nullclines need to be specified. Again, the readers can try the algorithm without the correction method by find_root_Newton. One issue about the above code is the lack of divergent check of the functions and Jacobians.
save(null1_cont, null2_cont, file = "./extra/data/03A/null_03A.RData") # we save the data for 03B