# A simplified 3-variable model

We saw in chapter 13 that a two-component negative feedback circuit (the auto-repressive gene system) possesses a unique steady state, which is always stable. This underlies homeostasis; the tendency of a negative feedback system to return to its equilibrium when perturbed (for small perturbations, at least).

However, negative feedback circuits can also generate sustained oscillations. One way this can occur is if we extend our two-variable model to a model including at least three variables. This was first proposed in the context of regulated gene expression in the 1960s, and the prototype three-component model is often referred to as the Goodwin oscillator (see chapter references). There are a number of ways of thinking about what might underlie a model of auto-regulatory gene expression that would contain more than two variables. One example would be that in fact the protein product of the gene is translated/produced in one place – the cytoplasm – but regulates transcription in another – the nucleus. We could therefore consider three variables: mRNA, cytoplasmic protein and nuclear protein.

Here, we will look at a slightly simplified model to Goodwin’s original, but the principles remain the same. Rather than focusing on a specific example of gene expression, we will consider a simple model of three interacting components, $X$, $Y$, and $Z$, which form a negative feedback loop by regulating each other. In particular:

• Increasing $X$ leads to increased production of $Y$,
• Increasing $Y$ leads to increased production of $Z$,
• Increasing $Z$ leads to decreased production of $X$.

Put together there is therefore a negative feedback in the system (since increasing $X$ ultimately leads to decreased production of $X$). For simplicity we will assume that the degradation rates of all three components are the same, so our model can be written as,

\begin{align} \frac{dX}{dt} &= f_1(Z) - \mu X\\ \frac{dY}{dt} &= f_2(X) - \mu Y\\ \frac{dZ}{dt} &= f_3(Y) - \mu Z. \end{align}

Each of the functions $f_i$ are non-negative and bounded as we assumed before. Given our assumptions in the bullet points we also have that $df_1/dZ\lt0$ and $df_2/dx,df_3/dY\gt0$.

Equilibria of this system occur where the following three conditions are simultaneously satisfied:

\begin{align} & \mu X = f_1(Z)\\ & \mu Y = f_2(X)\\ & \mu Z = f_3(Y). \end{align}

How many equilibria can we expect? If we let $F_i=f_i/\mu$, then we have $X=F_1(F_3(F_2(X)))=G(X)$. We know that $G(X)$ must be a decreasing function of $X$ (intuitively it must be because we know we have a negative feedback loop, but we can also show it by the chain rule). By the same argument as in the two-component auto-repression model, then, we can say that there must be a single equilibrium (because $X$ is an increasing straight line and $G(X)$ is some decreasing function, they can only cross once).

## Linear stability analysis

The Jacobian for this system may be 3 dimensional but doesn’t look too horrid,

\begin{align*} J=&\left( \begin{array}{ccc} -\mu & 0 & \phi_1\\ \phi_2 & -\mu & 0 \\ 0 & \phi_3 &-\mu \\ \end{array} \right), \end{align*}

where $\phi_i$ are the derivatives of $f_i$ at the equilibrium point. We have seen before that 3-dimensional Jacobians cannot be solved as easily as other models we have met. However, because we kept this system (relatively) simple we can actually check its stability directly from the characteristic equation, a bit like when we looked at the HIV model. To find this we take,

$\begin{equation} \det(J-\lambda I)=(-\mu-\lambda)^3+\phi_1\phi_2\phi_3. \end{equation}$

We then look for solutions for $\lambda$ to where this equation equals 0. That is,

$\begin{equation} (\mu+\lambda)^3=\phi_1\phi_2\phi_3=-\chi^30. \end{equation}$ Here we have made the product of our three gradients equal to a new parameter $-\chi^3$. The minus sign is to remind us that this product is negative (so in fact $\chi\gt0$) and the cubed is to remind us that this is really the composite of three parameters. The solution to this is then, $\begin{equation} \lambda=(-1)^{1/3}\chi-\mu. \end{equation}$ Somewhat susprisingly then, we find ourselves needing to find the cubed root of -1 in order to understand a real-world biological problem. Recalling that we can write $e^{-i\pi}=-1$ (I bet you didn’t expect to see Euler’s relation in this textbook!), we have that the three possible solutions are,

• $\lambda_1 = -\chi-\mu$,
• $\lambda_2 = e^{i\pi/3}\chi-\mu$,
• $\lambda_3 = e^{-i\pi/3}\chi-\mu$.

We can see that $\lambda_1$ is definitley negative. However the other two eigenvalues are complex, and we can find that they have real part $Re(\lambda_{2,3})=-\mu+\chi/2$. Putting this all together we find that,

• if $\chi\lt2\mu$, all the eigenvalues have negative real part and the equilibrium is a stable spiral (we know it is a spiral because two of the eigenvalues are complex),
• if $\chi\gt2\mu$, two of the eigenvalues have positive real part and the equilibrium is an unstable spiral.

We have seen a case before where stability moves from a stable spiral to an unstable spiral, when we looked at the predator-prey model. You may recall this produced a Hopf bifurcation, where a limit cycle emerges. The same phenomenon occurs here, with sustained oscillations in the expression of the gene arising. Time-courses and phase portraits of the dynamics before and after the Hopf bifurcation are shown in the figure below, with large oscillations in concentrations shown in the lower plots.  Time-courses and phase portraits of the three-variable model using the equations shown in the ‘Explore the model’ box below. In both cases, $\mu=, \theta=0.1, X(0)=0, Y(0)=0, Z(0)=0$. Top: $n=1$, giving a stable spiral into an equilibrium. Bottom: $n=4$, giving a limit cycle.

The three-component feedback loop we have considered here demonstrates that negative feedbacks can do more than homeostasis (stability of a system around a fixed point); it can also result in sustained oscillations. Such oscillations are a central feature of many biological systems. A prominent example is the roughly 24-hour period circadian oscillatory gene expression that occurs in both plants and animals. These oscillations, which are critical for coordinating metabolic activity, etc. with the rotation of the Earth, are generated by a slight elaboration of the mechanism we have studied here. Indeed, the Goodwin oscillator provided an early model for studying the properties of circadian oscillations. Other oscillations, with a wide range of oscillatory periods, are also generated by negative feedbacks.

Explore the model

Using the Python code below, explore how sustained oscillations can emerge in the three variable negative feedback model. The dynamics are governed by the following ordinary differential equations:

$\frac{dX}{dt}=-\mu+\frac{\theta^n}{\theta^n+Z^n}$
$\frac{dY}{dt}=-\mu+\frac{X^n}{\theta^n+X^n}$
$\frac{dZ}{dt}=-\mu+\frac{Y^n}{\theta^n+Y^n}$
Initially we have set $n=1$, and the resulting dynamics are for a slight spiral into the stable equilibrium. Gradually increase the value of $n$ in steps of 0.5 to see the limit cycle emerge and grow.
Click for code

`# Import the necessary libraries`

`import numpy as np`
`import matplotlib.pyplot as plt`
`from scipy.integrate import solve_ivp`
`# Options to make the plots the right size`
`plt.rcParams['figure.figsize'] = [12, 4]`
`plt.rcParams.update({'font.size': 16})`
`# Function for three variable model called 'goodwin'`
`def goodwin(t,x):`
`    # Rename the variables for ease`
`    X=x`
`    Y=x`
`    Z=x`
`    # The ODEs`
`    dX = -mu*X + theta**n/(theta**n+Z**n)`
`    dY = -mu*Y + X**n/(theta**n+X**n)`
`    dZ = -mu*Z + Y**n/(theta**n+Y**n)`
`    return [dX,dY,dZ]`
`# Parameter values`
`mu=1`
`theta=0.1`
`n=1`
`# Initial conditions`
`X0=0`
`Y0=0`
`Z0=0`
`N0=[X0,Y0,Z0]`
`# Time-points to use`
`tc = np.linspace(0, 50, 1000)  `
`# Use 'solve_ivp' to run the model`
`Nc = solve_ivp(goodwin, [tc,tc[-1]],N0, t_eval=tc)`
`# Plotting commands`
`fig, (ax1, ax2) = plt.subplots(1, 2)`
`ax1.plot(tc, Nc.y, "r", label="X") `
`ax1.plot(tc, Nc.y, "k", label="Y")`
`ax1.plot(tc, Nc.y, "b", label="Z")`
`ax1.set(xlabel='Time', ylabel='Concentrations')`
`ax1.legend()`
`ax1.axis([0,50,0,1])`
`ax2.plot(Nc.y,Nc.y,'b')`
`ax2.axis([0, 2, 0, 2])`
`ax2.set(xlabel='X', ylabel='Y')`
`ax2.axis([0,1,0,1])`

Key Takeaways

• We can create a 3 variable negative feedback circuit with two positive regulatory functions and one negative.
• In a simplified system we can calculate the eigenvalues directly.
• We find limit cycles can emerge for steep enough regulation funcions, such that gene expression regularly oscillates.

### Chapter references

• This model is a simplified version of that proposed by Goodwin (1965).
• The content in the Gene networks section is based on the unpublished Mathematical Biology lecture notes developed by Nick Monk.

## Licence 