Multistability in the lactose ultization network
The equation in the OP has been updated to $(y-a)^2(y-a\theta)$ which yields the correct expansion. Therefore, by equating the coefficients, we obtain
\begin{align}
\alpha &amp;= (2+\theta)a\tag{1}\\
\rho/\beta^2 &amp;= a^2(1+2\theta)\tag{2}\\
\alpha/\beta^2 &amp;= a^3\theta\tag{3}
\end{align}
Using equations (1) to (3), we get the desired $\rho$ and $\alpha\beta$. Now the authors plotted $\frac{\alpha\beta}{\rho}$ and $\frac{1}{\rho}$.

By plotting $\frac{\alpha\beta}{\rho}$ and $\frac{1}{\rho}$, we can re-create the plot in Mathematica with

ParametricPlot[{1/((1 + 2 x) (1 + 2/x)), (2 + x)^(3/2)/(
  Sqrt[x] (1 + 2 x) (1 + 2/x))}, {x, 0.01, 50 \[Pi]}, 
 PlotRange -&gt; {{0, .15}, {0, 0.6}}, AspectRatio -&gt; 3/4]




If you dont have Mathematica, here is some Python code that will plot it as well.

import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
t = np.linspace(0.01, 25*np.pi, 5000)
x = 1/((1+2*t)*(1+2/t))
y = (2+t)**1.5/(np.sqrt(t)*(1+2*t)*(1+2/t))
plt.plot(x,y)

plt.show()


Here is the Python plot which has smoother lines compared to Mathematica's plot.