Diagnosing Timescale Separation with Jacobian Spectra

Why Inspect the Jacobian Spectrum?

Nondimensionalization can suggest which processes may be fast or slow from their characteristic rates. A complementary question is whether the coupled model dynamics actually have distinct local modes. The Jacobian matrix describes how small state perturbations evolve near a specified state, and its eigenvalues give a local first indication of the modes’ rates of growth, decay, or oscillation.

In a QSP model, that diagnosis is local: a spectrum calculated at one state and parameter vector need not describe the entire trajectory or another virtual patient. It is evidence for candidate timescale separation, not by itself a proof that a reduced model will be accurate.

Definition: For an ODE model \(\dot{x}=f(x;\theta)\), where \(x\) is the state vector and \(\theta\) is a fixed parameter vector, the state Jacobian at state \(x_*\) is \(J(x_*;\theta)=\partial f/\partial x\) evaluated at \(x_*\). Its eigenvalues \(\lambda_i\) describe the local linearized modes. For a stable mode with \(\operatorname{Re}(\lambda_i)<0\), its local e-folding decay timescale is \(\tau_i=1/|\operatorname{Re}(\lambda_i)|\). The decay-rate magnitude is \(|\operatorname{Re}(\lambda_i)|\), not generally \(|\lambda_i|\).

If a mode has a complex eigenvalue \(\lambda=a+ib\), its local amplitude has an exponential envelope governed by \(a=\operatorname{Re}(\lambda)\), while \(b=\operatorname{Im}(\lambda)\) sets its angular oscillation frequency. When \(a<0\), the envelope timescale is \(1/|a|\); if \(b\ne 0\), the oscillation period is \(2\pi/|b|\). These are distinct quantities.

From Eigenvalues to Candidate Timescales

To interpret a spectrum, first check the signs of the real parts. If all are negative, the frozen linearized system is locally decaying. If any real part is positive, that mode grows in the linearization; assigning it a decay timescale would be inappropriate. A zero real part gives no finite exponential decay time, and may indicate a conserved, neutral, or marginal mode that needs separate interpretation.

For stable modes, the mode with the largest decay-rate magnitude is the fastest, and the one with the smallest positive decay-rate magnitude is the slowest. A convenient local ratio is

R_{\mathrm{local}}= \frac{|\operatorname{Re}(\lambda_{\mathrm{fast}})|} {|\operatorname{Re}(\lambda_{\mathrm{slow}})|} = \frac{\tau_{\mathrm{slow}}}{\tau_{\mathrm{fast}}}.

A large ratio indicates a candidate gap between those local exponential decay times. There is no universal cutoff that guarantees a useful reduction: the acceptable gap depends on the desired accuracy, the outputs of interest, the duration of interest, and whether the relevant modes can be interpreted and separated in the full model.

Worked Example: Two Real Stable Modes

Consider a hypothetical two-state model at a specified state and parameter vector with local Jacobian

J= \begin{pmatrix} -10 & 0\\ 0 & -0.2 \end{pmatrix} \ \mathrm{h}^{-1}.

The eigenvalues are \(\lambda_1=-10\ \mathrm{h}^{-1}\) and \(\lambda_2=-0.2\ \mathrm{h}^{-1}\). Both are negative, so both modes decay in the local linearized system. Their e-folding times are

\tau_1=\frac{1}{10\ \mathrm{h}^{-1}}=0.1\ \mathrm{h}, \qquad \tau_2=\frac{1}{0.2\ \mathrm{h}^{-1}}=5\ \mathrm{h}.

The first mode is the candidate fast mode and the second is the candidate slow mode. Their timescale ratio is \(R_{\mathrm{local}}=5/0.1=50\). Thus, at this local linearization, the slow mode’s e-folding time is 50 times the fast mode’s. This calculation alone does not show that the same ratio holds elsewhere in the model.

Justification: Why the Real Part Sets Decay

The timescale interpretation follows directly from the local linearization. For a small perturbation \(\delta x\) about a fixed point \(x_*\), the first-order approximation is \(\dot{\delta x}=J(x_*;\theta)\delta x\). If \(v\) is an eigenvector with eigenvalue \(\lambda\), a perturbation in that eigenvector direction has the form \(\delta x(t)=c e^{\lambda t}v\), where \(c\) is its initial modal amplitude. Writing \(\lambda=a+ib\) gives

e^{\lambda t}=e^{at}e^{ibt}, \qquad |e^{\lambda t}|=e^{at}.

The factor \(e^{ibt}\) rotates or oscillates but does not change the amplitude; the envelope changes as \(e^{at}\). For \(a<0\), one e-fold of decay occurs when \(e^{at}=e^{-1}\), which gives \(t=1/|a|\). This also explains why using \(1/|\lambda|\) as the decay time for a complex eigenvalue is generally wrong: the imaginary part affects oscillation, not the exponential decay envelope.

Worked Example: A Damped Oscillatory Pair and a Slow Mode

Suppose a hypothetical model has the local Jacobian

J= \begin{pmatrix} -2 & -20 & 0\\ 20 & -2 & 0\\ 0 & 0 & -0.1 \end{pmatrix} \ \mathrm{h}^{-1}.

The upper-left block has eigenvalues \(-2\pm20i\ \mathrm{h}^{-1}\), and the third eigenvalue is \(-0.1\ \mathrm{h}^{-1}\). All real parts are negative, so the modes decay locally. The oscillatory pair has envelope timescale

\tau_{\mathrm{osc}}=\frac{1}{2\ \mathrm{h}^{-1}}=0.5\ \mathrm{h}, \qquad T_{\mathrm{osc}}=\frac{2\pi}{20\ \mathrm{h}^{-1}} \approx 0.3142\ \mathrm{h}.

The remaining mode has timescale \(1/(0.1\ \mathrm{h}^{-1})=10\ \mathrm{h}\). Comparing its decay time with the oscillatory pair’s envelope time gives \(R_{\mathrm{local}}=10/0.5=20\). The pair oscillates about 20 times faster in its envelope-decay timescale comparison; its period, \(0.3142\ \mathrm{h}\), is a separate measure. Taking the modulus of either complex eigenvalue would mix damping and oscillation and would not give the decay rate.

State and Parameter Dependence

For a nonlinear model, \(J(x;\theta)\) can change as the state changes, even if \(\theta\) is held fixed. Changing parameters can also alter the eigenvalues, their eigenvectors, or both. A single spectrum is therefore a local snapshot. To investigate a trajectory or a population of parameter sets, evaluate spectra at representative states and parameter vectors, including regions reached after dosing, during transitions, and near relevant outputs or events.

The following hypothetical model makes the dependence explicit. Let \(x\) be a concentration in nM and \(z\) a dimensionless state, with time measured in hours:

\dot{x}=-\left(1+a\frac{z}{1+z}\right)x, \qquad \dot{z}=-0.2z.

Here \(a\) is dimensionless and positive. At \(x=0\), the Jacobian is triangular with eigenvalues \(-[1+az/(1+z)]\ \mathrm{h}^{-1}\) and \(-0.2\ \mathrm{h}^{-1}\). The next calculation compares spectra at different local conditions; it does not assume those conditions persist over time.

Worked Example: Separation Changes With State and Parameter

First take \(a=9\). At the local state \(x=0,\ z=0.1\), the rate for the \(x\) mode is \(1+9(0.1/1.1)=1.818\ \mathrm{h}^{-1}\), while the \(z\) mode’s rate is \(0.2\ \mathrm{h}^{-1}\). The decay times are \(1/1.818=0.55\ \mathrm{h}\) and \(1/0.2=5\ \mathrm{h}\), respectively, giving \(R_{\mathrm{local}}=5/0.55\approx9.091\).

At \(x=0,\ z=9\) with the same parameter, the \(x\)-mode rate is \(1+9(9/10)=9.1\ \mathrm{h}^{-1}\). Its timescale is \(1/9.1\approx0.1099\ \mathrm{h}\); the \(z\)-mode timescale remains \(5\ \mathrm{h}\). The ratio is \(5/0.1099\approx45.5\). The local separation is therefore much larger at this state than at \(z=0.1\).

Parameter changes matter as well. At \(z=1\), changing \(a\) from 1 to 19 changes the \(x\)-mode rate from \(1+1/2=1.5\) to \(1+19/2=10.5\ \mathrm{h}^{-1}\). The corresponding ratios to the \(z\)-mode rate are \(1.5/0.2=7.5\) and \(10.5/0.2=52.5\). Thus a model may have modest or substantial local spectral separation depending on both state and parameter regime.

Diagnostics, Alternatives, and Limitations

A practical spectral diagnosis should be repeated across the regions and parameter sets that matter. Check the following before treating a gap as meaningful:

Stability and near-zero modes: identify positive or nearly zero real parts rather than labeling every eigenvalue “fast” or “slow.” A growing mode is not a decaying fast process.
State and parameter coverage: track the eigenvalues along relevant trajectories and across plausible parameter regimes. A gap that closes in an important region is not a uniform separation.
Numerical reliability: verify the Jacobian calculation, units, and eigenvalue computation. Compare analytic derivatives with finite-difference checks when appropriate, and check that the chosen perturbation sizes are numerically sensible.
Mode interpretation: inspect eigenvectors to see which state combinations participate in each mode. A rate gap is less informative if the apparent modes are strongly mixed or do not align with the biological quantities of interest.

Eigenvalues alone can also miss transient behavior. A Jacobian is nonnormal when it does not commute with its transpose (for a real matrix), and its eigenvectors may be nonorthogonal. In that case, perturbations can temporarily grow or shift substantially even when every eigenvalue has a negative real part. When this is plausible, inspect short-time perturbation responses or singular-value amplification in addition to the eigenvalues.

The main alternatives answer related but different questions. Nondimensionalization, as in the earlier material on identifying fast and slow processes, uses scaled equations and characteristic rates to clarify units and candidate balances. Jacobian spectra include local coupling and provide modal rates, but are state-specific. Directly perturbing states or parameters and simulating trajectories can show whether outputs of interest respond on distinct time windows, but it may be computationally heavier and depends on the selected perturbations and outputs. These checks complement rather than replace one another.

A local eigenvalue gap does not establish that a reduced model will preserve important outputs, remain accurate over a dose regimen, or behave correctly across patients. The linearization is a first-order description near its evaluation state; nonlinear trajectories may move into regions with different spectra. Spectral diagnosis is best used to identify and scrutinize candidate modes, not as a stand-alone reduction rule.

Common Pitfall

Common mistake: For a complex eigenvalue, treating \(|\lambda|\) as the decay rate and calculating a decay time of \(1/|\lambda|\). The correction is to use \(1/|\operatorname{Re}(\lambda)|\) for a stable mode’s exponential envelope, and to use \(|\operatorname{Im}(\lambda)|\) separately for its oscillation frequency. If the real part is zero, there is no finite exponential decay time.

Check Your Understanding

Use the definitions and calculations above to answer these questions. Brief solutions follow the problems.

A local stable mode has eigenvalue \(-4\ \mathrm{h}^{-1}\). What is its e-folding timescale?
A mode has eigenvalue \(-0.5+3i\ \mathrm{h}^{-1}\). State its envelope timescale and oscillation period.
Two stable real modes have eigenvalues \(-8\ \mathrm{h}^{-1}\) and \(-0.4\ \mathrm{h}^{-1}\). Find their timescale ratio and identify the candidate fast mode.
Why is a spectrum calculated at one state insufficient to establish separation throughout a nonlinear QSP trajectory?
Can all-negative eigenvalue real parts rule out transient amplification in a nonnormal system? Explain.

Solutions. (1) \(\tau=1/(4\ \mathrm{h}^{-1})=0.25\ \mathrm{h}\). (2) The envelope timescale is \(1/0.5=2\ \mathrm{h}\); the period is \(2\pi/3\approx2.094\ \mathrm{h}\). (3) The timescales are \(1/8=0.125\ \mathrm{h}\) and \(1/0.4=2.5\ \mathrm{h}\), so the ratio is \(2.5/0.125=20\); the \(-8\ \mathrm{h}^{-1}\) mode is the candidate fast mode. (4) The Jacobian, and hence its eigenvalues, may change with state and parameters, so the local ratio can change along the trajectory. (5) No. Nonorthogonal eigenvectors can permit temporary perturbation growth even when every eigenvalue has a negative real part.