How to determine new dimensionless variables when non-dimensionsionalizing a system of ODEs? I am following this a research paper called "Evolution of Within-Host Antibiotic Resistance in Gonorrhea", which can be found here. I have one question regarding the method of non-dimensionalization. How do they come up with the new dimensionless state and time variables? 
Here's what the parameters and variables mean: $S$, $R_l$, and $R_h$ are the bacteria strains. $b_s$, $b_l$, and $b_h$ are birth rates of the bacteria strains. $\mu$ is the natural death rate. $i$ is the immune system clearing rate. $T$ is the treatment clearing rate. $m_1$, $m_2$, $p$ are mutation rates. $\alpha$ is the level of drug resistance.
The system is as follows,
\begin{align}
\frac{dS}{dt} &= S(b_s - \mu - i - T - m_1) - pS(R_l + R_h) \\
\frac{dR_l}{dt} &= R_l(b_l - \mu - i - \alpha T - m_2) + m_1 S + pR_l(S - R_h) \\ 
\frac{dR_h}{dt} &= R_h(b_h - \mu - i - \alpha T) + m_2R_l + pR_h(S + R_l) 
\end{align}
Their new dimensionless state and time variables are:
\begin{array}{cc} 
s = S\frac{p}{m_1} & l = R_l\frac{p}{m_1} & h = R_h\frac{p}{m_1} & \tau = tm_1 \\ 
\end{array}
With these new variables, the new system is 
\begin{align}
\frac{ds}{d\tau} &= sA - s(l+h) \\
\frac{dl}{d\tau} &= l(B-M) + s + l(s-h)\\
\frac{dh}{d\tau} &= hC + lM + h(s+l)
\end{align}
where 
\begin{array}{cc}
A = \frac{a}{m_1} - 1 & B = \frac{b}{m_1} & C = \frac{c}{m_1} & M = \frac{m_2}{m_1} \\
\end{array}
\begin{array}{cc}
a = b_s - \mu - i - T & b = b_l - \mu - u -\alpha T & c = b_h - \mu - i - \alpha T
\end{array}
I know how to end up with the new model by using the fact that $$\frac{dS}{dt} = \frac{dS}{d\tau} \frac{d\tau}{dt}$$
How do they know to make $s = S\frac{p}{m_1}$ and so on?
 A: Take a look at the Buckingham PI theorem, which provides a way for computing sets of dimensionless parameters from the given variables, even if the form of the equation is still unknown. 
A: I don't think there is a hard and fast rule, but here is my take on what the authors have done.
The first thing you want to do when non-dimensionalising a problem is to think about what time scale you are interested in, or what time scale you can measure against.
In the paper referred to by the OP the authors have, for whatever reason, selected $m_1$ to be the base time scale. 
The transfer rate of $S$ to $R_l$ might be easily measurable or estimated.
The next thing the authors have done is scale the state variables based on the  nonlinear interaction rate $p/m_1$ so that the nonlinear terms have coefficient $1$. One reason why the authors might want to do this scaling is to push all the parameters into the linear terms of the ODE. This could be useful for looking at the linearisation of the trivial solution and any bifurcations that may involve it.
