How long does it take for a chemostat to reach equilibrium? A chemostat is a device to grow a cell culture in equilibrium.
If we denote by $x$ the cell density and by $s$ the density of nutrients, we can write down the following equations (1):
$$\frac{ds}{dt} = D(s_R - s) - \mu x /Y$$
$$\frac{dx}{dt} = (\mu - D)x$$
where $t$ denotes time. Here $D$ and $Y$ are constants called dilution rate and yield constant respectively. $\mu$ is the specific growth rate and it can be a complicated function of $s$. However, two characteristics of $\mu(s)$ are generally accepted: 1. it is monotonically increasing; and 2. it saturates, that is, there exists an asymptote limit $\mu_\text{max} = \lim_{s\rightarrow\infty} \mu(s)$.
Usually, a so-called Monod function:
$$\mu(s) = \mu_\text{max}\frac{s}{k+s}$$
with suitable parameter values for $\mu_\text{max}$ and $k$ gives a good approximation.
This dynamical system has a non-trivial point of equilibrium (the trivial equilibrium point being $x=0,s=s_R$):
$$\tilde{s}=\mu^{-1}\left(D\right),\quad\tilde{x}=Y\left(s_{R}-\tilde{s}\right)$$
Here one assumes that $\mu(s_R)>D$. This equilibrium point is stable.
Suppose the chemostat is initially in equilibrium, and suddenly the parameters $D$ and $s_R$ change. What's a good estimate of how long it takes the chemostat to reach equilibrium again? I accept answers in the form of references to the literature.
(1) Szilard, L. (2001). Nonlinear population dynamics in the chemostat. Computing in Science & Engineering, 48–55.
 A: This sort of thing is usually calculated using linear stability theory. Essentially one replaces a set nonlinear equations with linear ones that approximate them in the region of the fixed point. The result is a set of dynamical equations of the form
$$
\mathbf{\dot y} = J \mathbf{y}, \tag{1}
$$
where $x_1 = s-\tilde s$, $y_2 = x - \tilde y$, and $J$ is the Jacobian matrix of the original system, given by
$$
J = \left( \begin{array}{cc}
\frac{\partial \dot s}{\partial s} & \frac{\partial \dot s}{\partial x} \\
\frac{\partial \dot x}{\partial s} & \frac{\partial \dot x}{\partial x} 
\end{array}\right).
$$
I hope you will forgive me for not working through the algebra, but it should be a simple matter to calculate the components of $J$ in terms of $x$ and $s$.
The point of doing this is that the stability and return time of a linear system can easily be calculated by examining the eigenvalues of $J$. If any eigenvalue has a positive real component then the fixed point will be unstable, but in your case (as long as the parameters are sensible) they will all have negative real parts.
Now, the solutions to equation $(1)$ have the form
$$
\mathbf{y}(t) = \sum_j A_j \mathbf{v}_j e^{\lambda_j t},
$$
where and $\lambda_j$ and $\mathbf{v}_j$ are eigenvalues and eigenvectors of $J$, respectively and the constants $A_j$ are determined by the initial conditions. The $\lambda_j$ will sometimes be complex, and if we express them in the form $\lambda_j = r_j + i\omega_j$ the solution becomes
$$
\mathbf{y}(t) = \sum_j A_j \mathbf{v}_j e^{r_j t} e^{i\omega_j t} = \sum_j A_j \mathbf{v}_j e^{r_j t} (\cos \omega_j t + i\sin \omega_j t).
$$
Complex eigenvalues always come in conjugate pairs, so that the $i\sin\omega_j t$ terms cancel out.
We're left with a sum of exponentially decaying (since each $r_j$ is negative) and possibly oscillating trajectories. The ones with smaller $r_j$'s will decay more rapidly, so after enough time the trajectory will be dominated by a term that decays as $e^{rt}$, where $r$ is the real part of the eigenvalue with largest real part.
We may therefore say that, if the perturbation is small enough, the original nonlinear system will decay towards equilibrium with a characteristic time given by $-1/r$, where $r=\max_j \{ \Re(\lambda_j) \}$. This is (approximately) the time taken for the distance from equilibrium to decay to $1/e$ of its original value, and it's what's usually quoted as the return time.
You should find all of the above in any decent text book on dynamical systems theory. Since your system is two-dimensional you can get an analytic expression for $r$ by writing down the characteristic polynomial for $J$ and using the quadratic formula to solve it. Whether this will result in a nice expression or not I don't know.
Now, there's a couple of things left to consider. Firstly, you asked about a change in the parameters rather than a small perturbation to the system. However, if the change in parameters is small then this doesn't make any difference: the steady-state configuration of the old system can be seen as a perturbation to the steady-state configuration of the new system, and everything follows as above. (Just make sure you use the new values of $D$ and $Y$ when calculating the eigenvalues of $J$.)
The other thing is that this theory isn't necessarily very accurate if the perturbation (or change to the parameters) is large. However, your system is not massively nonlinear, so I would expect it to be a decent approximation. If you're worried about that, the best thing is probably to numerically integrate the system with a few different parameter values and check it.
