Proper Orthogonal Decomposition for QSP Trajectory Reduction

Why Reduce Trajectories?

A QSP model may track many interacting states over time: cell populations, receptor complexes, cytokines, and downstream signals. Simulating this model repeatedly across parameter values and input conditions can be computationally expensive. Proper orthogonal decomposition (POD) offers a data-driven way to represent the state trajectories using a smaller set of coordinated patterns. The aim is not simply to keep a few original states. Instead, POD finds directions in the state space that capture much of the variation present in a chosen collection of simulations.

POD is a representation built from selected simulations, not an automatic guarantee that a reduced model will be accurate in every new setting. Its quality depends on which trajectories were sampled, how state variables were scaled, how many basis vectors were retained, and whether the projected dynamics preserve behavior important to the application.

From Simulations to a POD Basis

Start with the full model

Run it for a collection of sampled parameter and input conditions. At selected times, record each state vector as a column of a snapshot matrix. The time points need not be equally spaced, but the collection should represent the time regions and conditions relevant to the intended use. If a long trajectory is sampled at many nearly identical times, those snapshots can dominate the basis simply by their number.

Because QSP states can have different units and magnitudes, an unscaled Euclidean distance can give disproportionate influence to large-magnitude states. Choose a reference state \(x_{\mathrm{ref}}\), and a fixed positive-definite weighting matrix \(W\). A diagonal \(W\) can encode state scaling; for example, \(W_{ii}=1/s_i^2\) gives a state with characteristic scale \(s_i\) a comparable contribution to the squared distance. Let \(X\) contain columns \(x_j-x_{\mathrm{ref}}\), where \(x_j\) are the sampled states. The weighted squared reconstruction error for a matrix \(A\) is

Here \(A_{\cdot j}\) is column \(j\); this norm sums the weighted squared errors across snapshots. If all states are already on comparable scales, one may use \(W=I\), the identity matrix. The reference and scaling must be chosen consistently when building and using the basis.

To construct the basis, factor \(W=L^T L\), form the scaled snapshot matrix \(Y=LX\), and compute its singular value decomposition (SVD):

The columns of \(V\) are orthonormal directions in scaled state space. For a chosen dimension \(r\), map the first \(r\) directions back to the original state coordinates:

The columns of \(U_r\) are the POD basis vectors. A sampled state \(x\) is represented by \(x_{\mathrm{ref}}+U_ra\), where \(a\in\mathbb{R}^r\) contains its reduced coordinates.

A short justification makes clear what this result optimizes. In scaled coordinates, any candidate \(r\)-dimensional subspace captures the squared lengths of projections of the columns of \(Y\). The matrix \(YY^T\) has eigenvectors equal to the columns of \(V\) and eigenvalues \(\sigma_i^2\). For a projection onto any \(r\)-dimensional subspace, the total captured squared length is a weighted sum of these eigenvalues, with weights between zero and one whose sum is at most \(r\). Since the eigenvalues are ordered from largest to smallest, that sum cannot exceed \(\sum_{i=1}^{r}\sigma_i^2\); choosing the first \(r\) eigenvectors attains the bound. Total squared length is \(\sum_i\sigma_i^2\), so subtracting the maximum captured amount gives the stated minimum residual.

Projecting the Governing Equations

A basis can compress saved states, but a reduced dynamical model also needs equations for its coordinates. Approximate the state as \(x(t)\approx x_{\mathrm{ref}}+U_ra(t)\). In Galerkin projection, require the difference between the full model’s derivative and the derivative represented by the reduced model to be orthogonal to every retained basis vector in the weighted inner product. In symbols, the residual is \(f(x_{\mathrm{ref}}+U_ra,u,\theta)-U_r\dot a\), and its weighted inner product with the basis is zero. Since \(U_r^TWU_r=I_r\), the resulting reduced system is

The initial condition is the weighted orthogonal projection of the initial state onto the retained basis. If a different basis normalization is used, the general equation includes the matrix \(U_r^TWU_r\) on the left of \(\dot a\); it cannot be omitted unless that matrix is the identity. The projection gives a reduced state with \(r\) coordinates, but evaluating its right-hand side may still require evaluating the original nonlinear function \(f\). POD reduces the state dimension; it does not by itself make every function evaluation inexpensive.

Assessing Reconstruction Across Conditions

After selecting \(r\), assess how well the basis reconstructs sampled states. For a snapshot \(x_j\), its reduced coordinates and reconstruction are

A useful aggregate measure is the relative weighted Frobenius error across a specified set of snapshots:

This is a measure of state reconstruction, not necessarily a measure of error in an output of interest, nor a measure of the discrepancy between trajectories generated by integrating the full and reduced equations. Report the conditions and time points contributing to the calculation. A single aggregate can hide poor reconstruction for a particular parameter set, input, or short but important response feature.

Diagnostics, Alternatives, and Limitations

The singular values and reconstruction errors help diagnose whether a low-dimensional representation is plausible for the sampled trajectories. A sharp drop in singular values can indicate that a small basis captures most snapshot variation. A gradual decline suggests there may be no clear low-dimensional cutoff. Neither pattern alone certifies predictive accuracy: the answer depends on the simulations used to build the snapshot matrix.

Useful checks include computing reconstruction errors by condition and by time region, inspecting residuals in individual states, and confirming that the reference and scaling are appropriate. If errors concentrate in a small-magnitude state, reconsider whether its scale makes that error scientifically important. If errors rise sharply for some sampled inputs or parameters, the snapshot collection may not represent that region well, or the chosen dimension may be too small. Errors in a rapid transient can be obscured by many snapshots from a slowly changing period.

POD is not the only possible reduction strategy. Mechanistic aggregation combines states according to biological structure or conserved totals; it can preserve a meaningful interpretation but may not minimize the snapshot error. A POD basis is selected to minimize the specified aggregate squared reconstruction error, but its directions may combine states in ways that are less directly interpretable. Keeping selected original states is another simple option, but does not generally give the best low-dimensional approximation under the POD error criterion. These approaches answer different questions; low reconstruction error alone is not a reason to prefer one in every QSP application.

A linear POD subspace may also fail to respect physical restrictions such as nonnegative concentrations, bounded fractions, or conserved totals. The projection formula does not automatically enforce those restrictions. Reduced trajectories can leave the region represented by the snapshots, especially when the reduced dynamics are integrated over time. Such behavior calls for further model assessment or a different reduction design, not merely a larger number of reported digits in the snapshot error.

A practical sequence is to define the intended parameter and input range, select and scale representative trajectories, inspect the singular-value spectrum, choose a candidate dimension, and compute reconstruction errors both in aggregate and by condition. Then integrate the projected equations and compare the resulting behavior with the full model for the intended outputs and conditions. The latter comparison is a distinct check: the POD truncation theorem concerns the snapshot approximation, not the accuracy of reduced-model simulations.