Allow substitute() to perform simple reactants scaling and time scaling

It should be possible to create a new model via a coordinate transformation of an existing model.

Currently, this is only possible for the particular transformation: ` = - `, which reduces the number of reactants in the system by 1.

However, one should be able to scale reactants by a parameter value: ` = /`. This would be useful, for example, when taking the limit of the system size `N` to infinity. One could then scale the reactants by `N` to obtain the ODEs that describe how the fraction of each reactant change with time. Currently, this can be achieved by setting $N=1$ and reinterpreting the reactants as fractions of reactants.

Similarly, it would be advantageous if one could scale time in a model by a parameter value. This is common when reducing the number of parameters in a system of ODEs. Currently this can be achieved by setting the scaling parameter to 1 and reinterpreting the other parameters as `/`