Intermediate expressions in the `@ODEmodel` interface

Many models have subexpressions appearing in several places like

x1' = A * x2,
x2' = (A + b) * x1

where A = b^2 + b. A more practical example is given by SI-type models when one often has N being equal to the sum of all other variables.
Currently, MTK interface allows giving names to such subexpressions (via the observables mechanism) but the @ODEmodel interface would require the to type them explicitly each time. Would be great to have another way of doing this. Possible options:

1. Add some extra annotation to the @ODEmodel macro to introduce such quantities. I do not have an elegant solution off the top of my head. Could be something like

ode = @ODEmodel(
    x1'(t) = A * x2(t),
    x2'(t) = (A + b) * x1(t),
    y(t) = x1(t),
    where
     A = b^2 + b
)

1. Add some substitute functions which would allow to substitute some parameters/inputs in models with given expressions. This could look like

ode = @ODEmodel(
    x1'(t) = A * x2(t),
    x2'(t) = (A + b) * x1(t),
    y(t) = x1(t)
)
ode = substitute(ode, Dict(A => b^2 + b)).

In this case, it is not clear how to handle nicely the case when these expressions contain parameters/inputs which did not appear in the original model.

Overall, the desired features are:

1. simplicity and intuitivness
2. support of scalar and time-varying intermediate quantities
3. ability to add new symbols (parameter/inputs) during the substitution

Inspired by this issue

Start by tracing the @ODEmodel macro and the existing observables mechanism, then compare the proposed intermediate-quantity syntax with the substitute approach. Done should support scalar and time-varying quantities, keep the interface simple, and allow substitution to introduce new parameter or input symbols.