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

[renapagli]

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: <reactant> = <system size> - <reactants>, which reduces the number of reactants in the system by 1.

However, one should be able to scale reactants by a parameter value: <reactant> = <reactant>/<parameter>. 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 <parameter>/<scaling parameter>