Part 2: Ordinary differential equations
D. Stability
Mingyang Lu
10/26/2021
Stability of steady states
We start the discussion with the circuit of a constitutively expressed gene, whose dynamics is governed by the rate equation (Equation (1) from Part 2A):
\[ \frac{dX}{dt} = f(X) = g - kX \tag{1}\]
The steady states of a dynamical system must satisfy the condition \(\frac{dX}{dt} = 0\). In this case, there is one and only one steady state: \(X_s = \frac{g}{k}\). Now, let us evaluate the stability of the steady state. In general, a steady state is stable, if, after a small perturbation \(X \rightarrow X + \Delta X\) , the system goes back to the same steady state; otherwise, the steady state is unstable.
We separate the right hand side of Equation (1) into \(g\) and \(kX\) and plot them as the function of \(X\) below.
g = 50; k = 0.1 # Define parameters
X_all = seq(400, 600, 1) # All X values to be sampled
plot(X_all, array(g,dim=length(X_all)), type = "l", col =1,
xlab="X (nM)", ylab="Rates", xlim=c(400,600), ylim=c(40,60)) # Plot the first term g
lines(X_all, k * X_all, col=2) # Add the second term kX
legend("bottomright", inset=0.02, legend = c("g", "kX"),
col=1:2, lty=1, cex=0.8)
The intersect of the two curves represents the steady state, where the level of the steady state \(X_s = 500\) nM. Now, consider a small perturbation to increase \(X\) to \(X + \Delta X\) (\(\Delta X > 0\)). The black curve \(y=g\) stays constant, while the red curve \(y=kX\) slightly increases. That means \(\frac{dX}{dt} < 0\). Therefore, \(X\) level will decrease in next time steps. Similarly, if we make a small perturbation to decrease \(X\) to \(X - \Delta X\). We find \(\frac{dX}{dt} > 0\) and \(X\) levels will increase right after. Therefore, the steady state \(X_s\) is a stable state. In math,
\[\begin{equation} X_s \text{is} \begin{cases} stable & \frac{df(X)}{dX} < 0 \\ unstable & \frac{df(X)}{dX} \geq 0 \end{cases} \tag{2} \end{equation}\]
We will discuss about unstable steady states later. Here, in this case, \(\frac{df(X)}{dX} = -k\), which is negative. Thus, \(X_s\) is a stable steady state. Indeed, in Part 2A, we have shown that the system reaches to the same steady state \(X_s = 500\) nM, no matter which initial condition is selected (First plot in Part 2A).
Now, let us consider the circuit with a self-inhibiting gene, whose dynamics is governed by
\[\begin{equation} \frac{dX}{dt} = g(X) - kX \tag{3} \end{equation}\]
,where \(g(X) = g_0 + g_1 H^{inh}(X)\) (see Equations (10 - 11) in Part 2A).
We separate the right hand side of Equation (2) into \(g(X)\) and \(kX\) and plot them as the function of \(X\) below.
g0 = 10; g1 = 60; k = 0.1; X_th = 200; n = 4 # Define parameters
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))
}
X_all = seq(0, 600, 1) # All X values to be sampled
plot(X_all, g0 + g1 * hill_inh(X_all, X_th, n), type = "l", col =1,
xlab="X (nM)", ylab="Rates", xlim=c(0,600), ylim=c(0,70)) # Plot the first term g(X)
lines(X_all, k * X_all, col=2) # Add the second term kX
legend("top", inset=0.02, legend = c("g(X)", "kX"),
col=1:2, lty=1, cex=0.8)
The above plot shows an example of such two curves, who intersect exactly once. Therefore the circuit has only one steady state. It is possible to show that the circuit can only have one steady state, disregarding the choice of the model parameters. Similar to the first example, we can find that the steady state is also stable from either plotting or evaluation of Equation (2). Note that since \(H^{inh}\) monotonically decreases as \(X\) increases, thus its derivative with respect to \(X\) is always negative.
Let us also compute \(\frac{df(X)}{dX}\). Note that we use vector operations for high efficiency.
f_1g_self_inh <- function(X, g0, g1, X_th, n, k) { # Calculate f(X)
return(g0 + g1 * hill_inh(X, X_th, n) - k*X)
}
cal_dfdx <- function(X_all, dt, derivs, ...) { # Calculate df(X)/dt
# X_all is a vector of all X values to be sampled
# derivs takes the name of the derivatives function
# ..., ellipsis to take model parameters
f_all = derivs(X_all,...) # f(X) for all X values in a vector
f_all_shift_by_plus_one = c(f_all[-1], NA) # f(X+dx), achieved by shifting the vector by one to the left
f_all_shift_by_minus_one = c(NA, f_all[1:(length(f_all)-1)]) # f(X-dx), achieved by shifting the vector by one to the right
return((f_all_shift_by_plus_one - f_all_shift_by_minus_one)/2/dt) # df/dX = (f(X+dx) - f(X-dx))/2/dt
}
X_all = seq(0,600,0.1) # All X values to be sampled
dfdX = cal_dfdx(X_all = X_all, dt = 0.1, derivs = f_1g_self_inh,
g0 = 10, g1 = 60, X_th = 200, n = 4, k = 0.1)
plot(X_all, dfdX, type = "l", col =1,
xlab="X (nM)", ylab="df/dX", xlim=c(0,600), ylim=c(-0.5,0.1)) # Plot df(X)/dX
abline(a = 0, b = 0, lty=2, col = 2) # add a horizontal line for zeros
In this case, \(/frac{df(X)}{dX} < 0\) at the steady state X (around 250 nM). Thus, the steady state is stable.
Multi-stability
The third example is the circuit with a self-activating gene, whose dynamics is governed by
\[\begin{equation} \frac{dX}{dt} = g(X) - kX \end{equation}\]
,where \(g(X) = g_0 + g_1 H^{ex}(X)\) (see Equations (7 - 8) in Part 2A).
Again, we plot \(g(X)\) and \(kX\). In the plot below, we vary the value of \(k\) and fix the rest of the parameters.
g0 = 10; g1 = 45; X_th = 200; n = 4 # Define parameters
hill_ex <- function(X,X_th,n) { # Excitatory Hill function
# X_th: Hill threshold, n: Hill coefficient
a = (X/X_th)**n
return(a/(1+a))
}
X_all = seq(0, 600, 1) # All X values to be sampled
plot(X_all, g0 + g1 * hill_ex(X_all, X_th, n), type = "l", col =1,
xlab="X (nM)", ylab="Rates", xlim=c(0,600), ylim=c(0,60)) # Plot the first term g(X)
lines(X_all, 0.1 * X_all, col=2) # Plot the second term kX, k = 0.1
lines(X_all, 0.15 * X_all, col=3) # Plot the second term kX, k = 0.15
lines(X_all, 0.2 * X_all, col=4) # Plot the second term kX, k = 0.2
legend("topleft", inset=0.02, legend = c("g(X)", "kX (k = 0.1)", "kX (k = 0.15)", "kX (k = 0.2)"),
col=1:4, lty=1, cex=0.8)
When \(k = 0.1\) or \(k = 0.2\), there is one intersection of the two curves \(g(X)\) and \(kX\), so the circuit has one steady state. When \(k = 0.15\), there are three intersections, so the circuit has three steady states. From stability analysis, we can show the leftmost and the rightmost steady states are stable, while the middle steady state is unstable. We call this property of dynamical systems multistability.
We can intuitively understand the states of the self-activating circuit as follows. When \(X\) is low, there is no sufficient \(X\) product for transcription factor binding. Thus, the transcription rate is low, and \(X\) is kept at a low level. When \(X\) is sufficiently high, \(X\) binds to the promoter of itself most of time. Thus, the transcription rate is high, and \(X\) stays at a high level. This explains why the circuit can have two stable states at the same time.
We can also simulate the ODE with different initial conditions. The plot below shows that the time trajectories converge to the two stable steady states. We can also clearly observe that \(X\) of the unstable state is somewhere between 150 nM to 200 nM.
library(deSolve) # Use deSolve for ODE simulations here
parameters = c(g0 = 10, g1 = 45, X_th = 200, n = 4, k = 0.15) # Define all parameters for deSolve
dy_deSolve_1g_self_act<- function(t, y, parameters) { # Modified derivative function for deSolve
with(as.list(c(y, parameters)), {
return(list(g0 + g1 * hill_ex(y, X_th, n) - k*y)) # Output has to be a list type for deSolve
})
}
t_all = seq(0, 80, 0.1) # Time sequence for which output is needed
X0_all = seq(50,400,50) ## all initial conditions, each leads to a different simulation
legend_text = paste0("X0 = ", X0_all)
plot(NULL, xlab="t (Minute)", ylab="X (nM)", xlim=c(0,80), ylim=c(0,500)) # Initialize plotting
ind = 1
for(X in X0_all){
# ODE simulation per initial condition
ind = ind + 1
results = ode(y = X, times = t_all, func = dy_deSolve_1g_self_act, parms = parameters, method = "rk4")
lines(results[,1], results[,2], col=ind)
}
legend("bottomright", inset=0.02, title="Initial conditions",
legend = legend_text, col=2:9, lty=1, cex=0.8)
Again, we can also learn the stability of steady states from the \(\frac{df(X)}{dX}\) curve.
f_1g_self_act <- function(X, g0, g1, X_th, n, k) { # Define f(X) for a circuit of one gene with self-activation
return(g0 + g1 * hill_ex(X, X_th, n) - k*X)
}
X_all = seq(0,500,0.1) # All X values to be sampled
dfdX = cal_dfdx(X_all = X_all, dt = 0.1, derivs = f_1g_self_act,
g0 = 10, g1 = 45, X_th = 200, n = 4, k = 0.15)
plot(X_all, dfdX, type = "l", col =1,
xlab="X (nM)", ylab="df/dX", xlim=c(0,500), ylim=c(-0.2,0.2))
abline(a = 0, b = 0, lty=2, col = 2) # add a horizontal line for zeros
Effective potential
The concept of potential energy is widely used in science. In particular, for a molecular system, one can write down the potential functions to describe molecular interactions. However, a dynamical system describing a gene regulatory circuit is non-equilibrium in nature. For those systems, energy is not well defined. Nevertheless, it is possible to define effective potential for a gene circuit.
For a one-variable system, we can make an analogy between circuit dynamics governed by the rate equation
\[ \frac{dX}{dt} = f(X)\] with an overdamped particle in one dimension (1D). For the overdamped particle, \(X\) represents its position in 1D, \(-f(X)\) is the force, \(\frac{dX}{dt}\) represents the velocity of the particle. Here, velocity is proportional to the force for an overdamped system. Therefore, we can define the potential function
\[ U(X) = -\int_{X_0}^X f(x)dx + U(X_0) \] We will illustrate how to obtain the effective potential with numerical integration for the circuit with one self-activating gene. We set \(U(X_0 = 0) = 0\). A common numerical method to integrate a function is by the trapezoidal rule.
\[ \int_{X_0}^X f(x)dx = \sum_{x=X_0}^{X-\Delta x}\frac{f(x) + f(x+ \Delta x)}{2}\Delta x\] or
\[ U(X+\Delta x) = U(X) - \frac{f(X) + f(X+ \Delta x)}{2}\Delta x \] The first method uses a for loop to obtain U(X)
# Integrate f(X) with the trapezoidal rule
cal_int <- function(Xmin,Xmax,dX, f, ...) {
# Xmin: minimum X values
# Xmax: maximum X values
# dX: X step size
# f: function name
# ... : ellipsis to pass all model parameters to f(X, ...)
X_all = seq(Xmin, Xmax, dX)
nX = length(X_all)
# Initialize a vector of potential U
U_all = numeric(nX)
for(i in 1:(nX-1)) {
X = X_all[i]
U_all[i+1] = U_all[i] - (f(X,...) + f(X+dX,...))/2*dX
}
return(cbind(X_all,U_all))
}
results_1 = cal_int(Xmin = 0, Xmax = 500, dX = 0.1, f = f_1g_self_act,
g0 = 10, g1 = 45, X_th = 200, n = 4, k = 0.15)
plot(results_1[,1], results_1[,2], type = "l", col =1,
xlab="X (nM)", ylab="Effective Potential", xlim=c(0,500))
From the potential, there are two basins at around 100 nM and 300 nM, which are exactly the two stable steady states, and a barrier at around 200 nM, which is the unstable steady state. This provides an intuitive way to characterize a dynamical system with multistability. In addition, it also gives some ideas about how stable each state is. Here, the right basin has lower effective potential than the left basin. When there is sufficient gene expression noise, the system may transit from one stable state to another stable state. The effective potential also illustrates how easy it is to make such a transition. More details will be discussed in Parts 4 (Stochastic differential equations) & 6 (Monte Carlo simulations).
Now, let us do the numerical integration again in R, but using only vector operations.
cal_int_vector <- function(Xmin,Xmax,dX, f, ...) {
# Xmin: minimum X values
# Xmax: maximum X values
# dX: X step size
# f: function name
# ... : ellipsis to pass all model parameters to f(X, ...)
X_all = seq(Xmin, Xmax, dX)
f_all = f(X_all, ...) # generate all f(X)
f_all_shift_by_plus_one = c(f_all[-1],NA) # another array containing f(X+dx), except for the last point
f_all = (f_all + f_all_shift_by_plus_one)*dX/2 # each integrated part
U_all = -cumsum(f_all)
return(cbind(X_all, U_all)) # perform cumulative sum
}
results_2 = cal_int_vector(Xmin = 0, Xmax = 500, dX = 0.1, f = f_1g_self_act,
g0 = 10, g1 = 45, X_th = 200, n = 4, k = 0.15)
plot(results_2[,1], results_2[,2], type = "l", col =1,
xlab="X (nM)", ylab="Effective Potential", xlim=c(0,500))
We benchmark the efficiency of the two methods. The vector approach is more concise and efficient than the one with for loop.
library(microbenchmark) # use microbenchmark library for timing
my_benchmark = microbenchmark(
integration_normal = {r1 = cal_int(Xmin = 0, Xmax = 500, dX = 0.1, f = f_1g_self_act,
g0 = 10, g1 = 45, X_th = 200, n = 4, k = 0.15)},
integration_vector = {r2 = cal_int_vector(Xmin = 0, Xmax = 500, dX = 0.1, f = f_1g_self_act,
g0 = 10, g1 = 45, X_th = 200, n = 4, k = 0.15)},
times = 10, unit = "s")
microbenchmark:::autoplot.microbenchmark(my_benchmark)## Loading required namespace: ggplot2## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Finally, we show the effective potential for \(k = 0.2 \text{Minute}^{-1}\).
results_3 = cal_int_vector(Xmin = 0, Xmax = 500, dX = 0.1, f = f_1g_self_act,
g0 = 10, g1 = 45, X_th = 200, n = 4, k = 0.2)
plot(results_3[,1], results_3[,2], type = "l", col =1,
xlab="X (nM)", ylab="Effective Potential", xlim=c(0,500))
And the effective potential for \(k = 0.1 \text{Minute}^{-1}\) .
results_4 = cal_int_vector(Xmin = 0, Xmax = 500, dX = 0.1, f = f_1g_self_act,
g0 = 10, g1 = 45, X_th = 200, n = 4, k = 0.1)
plot(results_4[,1], results_4[,2], type = "l", col =1,
xlab="X (nM)", ylab="Effective Potential", xlim=c(0,500))