Stable Parameterizations for Multi-Pool Turnover Models

Why Parameterization and Scaling Matter

A linked turnover model can be biologically sensible and still be numerically difficult to fit. A solver may struggle when one pool changes on a minutes-long timescale and another on a days-long timescale; an optimizer may propose a negative conversion rate; or a single absolute error tolerance may be appropriate for one concentration but too loose for another. These computational problems can alter fitted concentrations and parameter estimates even when the biological model itself has not changed.

The goal of a stable parameterization is not to make a model fit better by changing its meaning. It is to represent the same model in variables and coordinates that are better suited to numerical integration and optimization. The key choices are state scales, positive-rate parameterizations, and solver tolerances. Each should be checked against the model’s units, anticipated concentration ranges, and observed predictions.

Scaling States Without Changing the Model

Consider \(m\) linked pool concentrations collected in \(x(t)\), where each \(x_i\) is measured in concentration units. A common model has an input vector \(u(t)\) and linear turnover and transfer terms:

Here, entries of \(u(t)\) have units of concentration per time, entries of \(A(k)\) have units of inverse time, and \(k\) denotes the model’s rate parameters. For a compartment-style system, transfer terms between pools have nonnegative coefficients, while diagonal terms generally represent the negative total rate of loss from each pool. Scaling is exact: it changes coordinates, not the predictions, provided the equations and initial conditions are transformed consistently.

Choose a positive concentration scale \(S_i\) for each pool and a positive time scale \(T\). Define \(z_i=x_i/S_i\) and \(\tau=t/T\). In vector notation, let \(D\) be the diagonal matrix with entries \(S_i\). Then \(x=Dz\), and the transformed model is:

Every term in this equation is dimensionless. In particular, a transfer coefficient from pool \(j\) to pool \(i\) is multiplied by \(T S_j/S_i\). This factor matters: scaling the states but forgetting the ratio \(S_j/S_i\) changes the model. Scales should be fixed reference values chosen to put typical states in a numerically useful range, not quantities that vary during integration.

Enforcing Positive Rates

Turnover, formation, conversion, and degradation rates are ordinarily nonnegative, and a rate that is assumed strictly positive should not become negative merely because an optimizer takes a step outside the feasible region. One simple approach is to estimate a real-valued parameter \(\eta\) and define \(k=k_{\mathrm{ref}}\exp(\eta)\), where \(k_{\mathrm{ref}}>0\) is a rate with the same units as \(k\). This guarantees \(k>0\) for any finite \(\eta\). Its derivative is \(dk/d\eta=k\), so changes in \(\eta\) act on a relative rather than additive rate scale.

When several outgoing pathways share a total loss rate, a useful alternative is to parameterize the total rate and the fractions assigned to each pathway. For positive total rate \(K\), define pathway fractions using a softmax transform:

Each fraction is positive and the fractions sum to one, so the pathway rates are positive and sum to \(K\). This can be convenient when the data are more informative about total loss than about its decomposition. It does not, however, create information that the data do not contain: weak identification of individual pathways remains a scientific and statistical limitation.

A log transform and a softmax transform are not interchangeable in every model. Use a log transform for a rate that is independently positive. Use a total-rate-and-fractions representation when the rates genuinely partition one total. If a parameter can equal zero under the scientific model, a strictly positive exponential transform excludes that boundary; the parameterization must reflect the intended model rather than impose positivity automatically.

Choosing Solver Tolerances and Checking Numerical Behavior

A numerical solver approximates the solution to the differential equations. In a common error-control convention, a step is judged against a scale of the form \(\mathrm{atol}_i+\mathrm{rtol}|x_i|\), where \(\mathrm{atol}_i\) is an absolute tolerance for state \(i\) and \(\mathrm{rtol}\) is a relative tolerance. The precise implementation varies by solver, and this expression is an error-control scale, not a guarantee that the actual error is always smaller at every reported time.

If states are scaled, the tolerances must be interpreted on the scaled scale. For \(z_i=x_i/S_i\), a dimensionless absolute tolerance \(\mathrm{atol}_{z,i}\) corresponds to a concentration-scale tolerance \(S_i\mathrm{atol}_{z,i}\) in \(x_i\). Thus a shared absolute tolerance for all dimensionless states can be appropriate when their scales were chosen sensibly. With unscaled states, a single absolute tolerance may be much too loose for a small pool or needlessly strict for a large one.

Justification: Scaling Preserves Predictions

The coordinate transformation preserves the model because it is invertible for positive scales. Let \(D\) be diagonal with positive entries \(S_i\), and set \(x=Dz\), \(\tau=t/T\). By the chain rule, \(dx/dt=(D/T)(dz/d\tau)\). Substitution into \(dx/dt=u(t)+A(k)x\) gives:

Recovering \(x=Dz\) therefore returns the original differential equation, including its initial condition \(z(0)=D^{-1}x(0)\). The same transformation applies to an observation equation, for example \(y=c^\mathsf{T}x=c^\mathsf{T}Dz\). A fit can nevertheless change numerically because optimization and integration use finite precision and stopping criteria. A changed fit after scaling is a reason to check implementation and numerical convergence, not evidence that scaling changed the biology.

Diagnostics, Alternatives, and Limitations

For a linear turnover system with nonnegative inputs and nonnegative initial concentrations, nonnegative transfer coefficients, and nonpositive diagonal loss terms, the exact solution should not produce negative concentrations. To see the role of the assumptions, the derivative of one pool at \(x_i=0\) is its nonnegative input plus nonnegative inflows from other pools; its outflow term is proportional to \(x_i\) and is then zero. The vector field therefore does not point out of the nonnegative region at a boundary. A small negative value from numerical integration can still occur, but a substantial or persistent negative prediction calls for investigation rather than silent clipping.

Diagnostics should assess the quantities used in the analysis, not just whether the solver reports success. Compare predictions at observed sampling times after tightening tolerances; inspect small and rapidly changing pools; check units and transformed initial conditions; and verify that predictions are stable under a reasonable change of solver or step controls. In parameter estimation, also check whether the fitted parameters or objective change materially when numerical accuracy is increased. Solver noise can be mistaken for a change in fit, particularly when local sensitivities are evaluated by small parameter perturbations.

A system is numerically stiff when its dynamics include substantially different timescales, such that a method stable for the fastest component may need very small steps to follow slower behavior. Scaling state magnitudes can improve error control, but it does not remove stiffness in the underlying dynamics. An implicit solver may help with stiffness, while suitable state scales and tolerances address a different issue. Neither solver choice nor positive-rate transformations repair an incorrect pathway structure, poor identifiability, or misspecified observations.

Common warning signs include fitted concentrations that change materially when tolerances are tightened, negative states larger than negligible numerical noise, abrupt oscillations unsupported by the model, and optimizer estimates that depend strongly on initial values or solver settings. Do not automatically clip negative states to zero: clipping changes the numerical trajectory and can create artificial discontinuities. Instead, assess their size, investigate the solver and model implementation, and report any remaining numerical sensitivity that affects interpretation.