Bistable systems
We already touched on bistable systems last time and I quickly want to discuss a basic bistable transcriptional system. $$ \frac{dX(t)}{dt} = \frac{x^n}{1+x^n} - \frac{ax}{b+x} $$ where the first term is autorepression and the second term corresponds to an enzyme mediated degradation.
Fixpoints and phase portraits
Consider a general two dimensional dynamical system $$ \frac{dx(t)}{dt} = f(x,y) \quad \mathrm{and}\quad\frac{dy(t)}{dt} = g(x,y) \ . $$ At any time, the future evolution of the system is fully specified by the current position \(x(t_0),y(t_0)\) and the future evolution of the system is a path \(x(t),y(t)\) in the plane. As a corollary, two such path can never cross: If they did, we would encounter the situation where the same state (the point at which the path cross) gives rise to different path. The fact that trajectories can't cross makes qualitative analysis of the phase plane of a two dimensional system a very powerful tool.
In the last lecture, we introduced nullclines as one dimensional objects in the phase plane where the derivative of one of the variables with respect to time vanishes, that is either \(f(x,y)=0\) or \(g(x,y)=0\). On a nullcline, the dynamics of the system is either completely horizontal of vertical in the phase plane. Furthermore, nullclines separate regions in phase space where \(x,y\) are increasing or decreasing.
The intersection of nullclines are fixpoints. We have already analyzed fixpoints of one dimensional systems. In this case, the behavior of the system in the vicinity of the fixpoints is generically exponential growth or decay. In two dimensions, the behavior around the fix point is again exponential, but the behavior is much richer.
Linear stability analysis in two dimensions
Assume the system has a fix point \(x_f, y_f\) and expand the dynamical equations in \(\delta x = x-x_f, \delta y = y-y_f\). $$ \frac{d\delta x(t)}{dt} = f(x_f,y_f) + \frac{\partial f(x_f, y_f)}{\partial x} \delta x + \frac{\partial f(x_f, y_f)}{\partial y} \delta y + \ldots $$ where \(f(x_f,y_f)=0\) by virtue of evaluating the function at the fixpoint. The analogous equation can be written down for \(y_f\) and the two can be combined into a matrix equation $$ \frac{d \delta\vec{x}}{dt} = \begin{pmatrix} \frac{\partial f(x_f, y_f)}{\partial x} & \frac{\partial f(x_f, y_f)}{\partial y} \\ \frac{\partial g(x_f, y_f)}{\partial x} & \frac{\partial g(x_f, y_f)}{\partial y} \end{pmatrix} \begin{pmatrix} \delta x\\ \delta y \end{pmatrix} + \mathcal{O}(\delta x^2) $$ In the vicinity of the fixpoints, we can ignore the higher order terms and the problem reduces to a linear equation with the matrix being the Jacobian of the system at the fix point.
Systems of linear equations
Let's simplify notation and consider $$ \frac{d \vec{x}}{dt} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \vec{x} $$ To solve this, we determine the eigenvalues and eigenvectors. Eigenvalues are the roots of the characteristic polynomial \(det(M - \lambda I)\): $$ \lambda_{1/2} = \frac{a+d\pm\sqrt{(a-d)^2 + 4bc}}{2} $$ Note that depending on the sign of \((a-d)^2 + 4bc\) the eigenvalues are real or complex. Given the eigenvalues, it is straightforward to determine the right eigenvalues from \(M\vec{x} = \lambda_{1/2}\vec{x}\). $$ \vec{v_{i}} = C \begin{pmatrix} \frac{\lambda_{i}-a}{b} \\ 1 \end{pmatrix} $$ Analogously, the left eigenvalues are given by $$ \vec{w_{i}} = D \begin{pmatrix} \frac{\lambda_{i}-a}{c} \\ 1 \end{pmatrix} $$ We will choose \(C\) and \(D\) such that the scalar product \(\vec{v_{1}}\vec{w_{1}}=1\). The product between left and right eigenvectors corresponding to different eigenvalues vanishes as expected. From $$ 4(\lambda_{1}-a)(\lambda_2-a) = 4\lambda_1\lambda_2 - 8a(\lambda_1+\lambda_2) + 4a^2 = (a+d)^2 - (a-d)^2 - 4cb - 4a(a+d) + 4a^2 = 4bc $$ follows that \(\vec{v_{1}}\vec{w_{2}}=0\).
To solve the equation, we expand the \(\vec{x}\) in right eigenmodes $$ \vec{x} = a_1(t)\vec{v_{1}} + a_2(t)\vec{v_{2}} $$ Plugging this into the ODE, we find that $$ a_i(t) = a_i^0 e^{\lambda_i t} $$ where \(a_i^0\) is determined by projecting the initial condition on the left eigenvectors. The general solution of a two dimensional linear system is hence a superposition of two exponentials. Depending on the eigenvalues, this solution has a number of qualitatively different solutions:
Complex eigenvalues \((a-d)^2 + 4bc<0\)
In this case the two eigenvalues are a complex conjugate pair with real part \(\lambda = a+d\) and imaginary part \(\omega = \pm \sqrt{(a-d)^2 + 4bc}\). The solution are hence spirals with a angular velocity \(\omega\). If the real part \(\lambda<0\), these spirals will decay and the fix point is stable. Otherwise they will grow. In the special case \(\lambda=0\) the solution are limit cycles, but in this case the behavior of the system is likely determined by higher order terms.
Real eigenvalues \((a-d)^2 + 4bc>0\)
If eigenvalues are real, the behavior of the systems depends on whether none, one, or both eigenvalues are positive. In the first case the fix point is stable, in the other two cases the fix point is unstable either in one or in both directions.
A concrete example
Let's use the example from the exercises. $$ \frac{dx(t)}{dt} = \frac{a}{1+y^n} - x \quad \mathrm{and}\quad \frac{dy(t)}{dt} = x - \frac{by}{1+y} $$ The Jacobian of this system is (after simplifying using the fixpoint conditions) $$ \begin{pmatrix} -1 & -\frac{any_f^{n-1}}{(1+y_f^n)^2} \\ 1 & -\frac{b}{1+y_f} + \frac{by_f}{(1+y_f)^2} \end{pmatrix} = \begin{pmatrix} -1 & -\frac{nx_f^2}{y_f} \\ 1 & \frac{x_f}{y_f}\left(-1 + \frac{x_f}{b}\right) \end{pmatrix} $$ For general parameters, we'll need to solve for fix point numerically and plug the values into the Jacobian to determine the behavior of the fix point. For \(a=4, b=3, n=3\), we find \(x_f = 1.6\ldots, y_f=1.14\ldots\) and the corresponding eigenvalues are $$ \lambda_{1/2} = -0.86 \pm 2.58i $$
Special case \(n=1\)
For \(n=1\), the fix point is at $$ y_f = \frac{a}{b}\quad \mathrm \quad x_f = \frac{ab}{a+b} $$ Plugging this into the Jacobian we find $$ \begin{pmatrix} -1 & -\frac{ab^3}{(a+b)^2} \\ 1 & -\frac{b^3}{(a+b)^2} \end{pmatrix} $$ Eigenvalues are therefore $$ \lambda_{1/2} = -\frac{1}{2}-\frac{b^3}{2(a+b)^2} \pm \sqrt{(1-\frac{b^3}{(a+b)^2})^2-\frac{ab^3}{(a+b)^2}} $$ This immediately tells us that the fix point is never going to be unstable.
Oscillatory circuits
It is rather difficult to get a generic biologically inspired two-dimensional system to exhibit stable oscillations. Three component systems, however, oscillate rather readily. So what is the crucial ingredient necessary to make a biological system oscillate? Oscillations naturally arise in a delayed negative feed back. Consider for example a protein that represses its own expression. The time delay between transcription and the mature protein implies that the mRNA expression at time \(t\) depends on the mRNA expression at some time \(t-\tau\). $$ \frac{dY(t)}{dt} = \frac{a}{1+Y^n(t-\tau)} - \frac{bY(t)}{1+Y(t)} $$ If the delay \(\tau\) and the parameters \(a,b\) fall into a certain range, this system will oscillate.
