Mathematical Pharmacology

0. Why PK, PD, and QSP fragmented

Pharmacokinetics is usually presented as a small, self-contained piece of linear algebra. A drug is administered into a compartment, a well-stirred pool in which every molecule of drug is equally likely to leave by any of the pool's available routes, so that a single number (its concentration) together with a handful of rate constants describes it completely. Drug moves between compartments by first-order exchange and leaves the body by first-order elimination, and the resulting differential equation is linear from the outset. The pharmacokinetics literature already writes this in matrix notation whenever it writes it formally at all: a vector of compartment concentrations, a matrix of rate constants, a solution built from a single exponential.

Pharmacodynamics is taught as something else. Rather than a differential equation, its standard object is a curve to fit: the Hill equation or the Emax model, four parameters estimated by nonlinear regression against a scatter of measured effect against dose. The dynamical system that produced that curve, if there is one, is left implicit or omitted outright, and the pharmacodynamics literature rarely gestures back at the linear machinery next door.

Quantitative systems pharmacology (QSP) looks like a third subject again, and a much larger one. A QSP model may carry tens to hundreds of state variables (free drug, bound drug, receptor, ligand-receptor complex, and a cascade of downstream signalling intermediates), coupled by nonlinear mass-action and Michaelis-Menten terms. Each new pathway is typically diagrammed and coded as its own bespoke construction, with little of the structure of the last model built the same way carried over.

These three literatures cite one another rarely and share mathematics less. They read as three different subjects because they are taught, reviewed, and published as three different subjects. They are not. Each is an instance of one object: a state that evolves under a control input and is read out through a measurement ,

Pharmacokinetics is the case in which and are linear. Pharmacodynamics, read correctly, is the case in which the reported curve is the output map evaluated at a fixed point of a one-compartment flow, not a value read off a moving trajectory (Section 13). QSP is the case in which and are left fully nonlinear and the state dimension is allowed to grow large. Nothing about the object in [eq:state-eq-preview] changes across these three cases; only the fragment of its theory that each field happens to import changes, a point we return to explicitly in Section 8.

Two anchors run through everything that follows, and they are worth naming before a single formal definition appears.

Anchor 1: control theory before physics. Every tool built in this Part (controllability, observability, BIBO stability, a positivity structure specific to compartmental sign patterns, and a structural test for observability that needs no numbers at all) was developed by control engineers across the 1950s and 1960s to answer questions about feedback systems, well before anyone asked those questions of a compartmental drug model. We import the finished theory rather than reconstructing pharmacology-specific analogues of it. The theory does not care what measures, and neither do we, until Part II.

Anchor 2: observability at scale. The technical stake of the article is what happens to controllability, and especially observability, once the state dimension grows from a handful of compartments into the range QSP models actually occupy. Does every state coordinate remain identifiable from the small number of outputs a laboratory can measure? Once forming a matrix by hand stops being the bottleneck, can the answer even be computed? Part IV answers both questions; every tool built before it exists to make that answer possible.

The construction proceeds from the ground up, in one direction only. Part I (this Part) builds the skeleton: the state-space object itself, its closed-form solution, controllability, observability, BIBO stability, the positivity structure that compartmental sign patterns force on the flow, and a structural observability test that needs only the wiring diagram of the system, never its numbers. Every statement is made for an arbitrary state dimension ; nothing here is fixed at any particular compartment count. Part II specialises the skeleton to the reversible/dissipative split that compartmental kinetics carries for free, recovers the pharmacodynamic dose-response curve as a fixed point rather than an independent object, and introduces a single worked pharmacokinetic example that the rest of the article returns to and grows, compartment by compartment. Part III lifts the same skeleton from a matrix to a tensor, the setting a QSP reaction network actually lives in. Part IV puts every tool to work at the scale QSP demands and closes with a numeric verdict on when an added compartment buys real dynamical structure and when it buys only a direction no measurement can see.

We assume the reader's linear algebra: vector spaces, matrices, rank, eigenvalues, and the ordinary differential equations built from them. We assume no prior pharmacology whatsoever; every pharmacological term, starting with compartment above, is defined in prose at the point it first does any mathematical work.

Part I: The control-theoretic skeleton

1. The state-space object

Section 0's three literatures are one object once its state dimension is left free. This section fixes that object once, first in the nonlinear generality it needs to hold PK, PD, and QSP simultaneously, then specialised to the linear form the rest of Part I works with.

In the compartmental reading that occupies the rest of the article, is the amount or concentration of drug in compartment , is the administered input (an infusion rate, a bolus, a dosing schedule), and is whatever a laboratory can actually measure, almost always a small number of compartments against a state dimension that may be large. The off-diagonal entry for is the rate at which drug leaves compartment for compartment ; the diagonal entry collects everything that leaves compartment by any route, including out of the system entirely (elimination, or clearance, the volume of a reference fluid fully cleared of drug per unit time). This sign structure, every off-diagonal entry nonnegative with the diagonal absorbing the total outflow, is not an accident of any particular model; Section 6 shows it holds of every compartmental system, for every , and derives its consequence for the sign of the whole flow.

plays no role in most compartmental models, since there is rarely a direct algebraic path from dose rate to a measured concentration that bypasses the state entirely, but we carry it throughout because it costs nothing to keep and because Section 4 and Section 5 need to say precisely what it does and does not affect. Sections 2 through 7 develop the linear specialisation [eq:state-eq-linear] in full; Part III returns to the fully nonlinear [eq:state-eq-general] for QSP networks.

2. Solution: the matrix exponential and superposition

Before any structural property of [eq:state-eq-linear] can be checked, we need the trajectory it produces. The linear system already has a closed-form solution, and every section from Section 3 onward is a statement about that formula rather than about the differential equation directly.

The solution in [eq:variation-of-parameters-eq] splits the trajectory into a free response , the state relaxing (or not, Section 5) on its own, and a forced response, the convolution of the input against the same exponential kernel. Superposition, the forced response to a sum of inputs is the sum of the forced responses to each, is immediate from linearity of the integral, and is the fact that lets pharmacokinetics reduce every dosing regimen (bolus, infusion, repeated dosing) to a sum or convolution of exponentials, at any compartment count , without ever re-deriving the underlying differential equation. The same formula, evaluated numerically rather than symbolically, produces the concentration-time curves of Section 11.

3. Controllability

The free response of Section 2 describes what the state does with no input at all. Controllability asks the opposite question, which states the input can reach in the first place, and turns that question into a finite rank computation.^[1]

4. Observability

Controllability asks what the input can do to the state; observability asks the reverse, what the output reveals about the state. The construction is exactly dual, and Section 3's theorem transposes without change.

5. BIBO stability

Every dosing regimen a clinician actually administers is bounded: a finite infusion rate, a finite number of finite doses. BIBO stability asks whether that alone is enough to guarantee a bounded measured response, a weaker and more operationally relevant question than whether the state itself decays.

This is deliberately weaker than asking that every eigenvalue of the full have negative real part, internal, or asymptotic, stability of . A mode that is uncontrollable, unobservable, or both is invisible to : it never appears in the measured output no matter what is asked of it or measured from it, and its eigenvalue cancels out of [eq:transfer-function] regardless of its sign. A system can therefore be BIBO stable while carrying an internally unstable mode that no input reaches and no output sees. This gap between what decays and what is seen to decay is the same gap Section 7 and Part IV return to under the name of Anchor 2, observability at scale. Section 10 shows this criterion collapses to a single positivity check once the reversible/dissipative split of Section 9 is in place.

6. Metzler matrices and positivity

Compartmental kinetics is not merely linear; it is linear with a sign pattern, since mass or concentration can neither go negative nor be created from nothing. This section isolates that sign pattern as a property of alone and shows it forces positivity on the whole flow, for every .

7. Structural (generic) observability

The rank test of Section 4 is decisive, but it needs numbers: a specific and to form and check its rank. At the compartment counts Part IV works at, forming that matrix is undesirable and its rank is fragile: a coincidence among the particular rate constants used can make a system look observable, or not, in a way that a small parameter change would reverse. This section replaces the numerical test with a combinatorial one that depends only on which entries of and are allowed to be nonzero, never on their values.

Fix, rather than a numerical and , only their sparsity pattern: which entries and (for ) are free parameters and which are structurally zero. Such a pair is a structured pair , and a numerical realisation is any choice of real values for its free entries. Associate to a structured pair the directed graph with a vertex for each state coordinate and each output coordinate , an edge whenever is a free parameter, and an edge whenever is a free parameter. is exactly the compartment diagram: an arrow for every rate constant the modeller has allowed to be nonzero, plus an arrow to every sensor.

Condition 1 rules out a state coordinate that no measurement can reach even in principle, however the free parameters are chosen; condition 2 rules out two or more state coordinates competing for the same limited routes to the outputs, a structural degeneracy no numerical coincidence can fix generically. Both conditions are decidable in time polynomial in (a saturating matching in a bipartite graph is found by, for example, the Hopcroft-Karp algorithm), which is the entire point: Section 21 and Part IV run this test at a compartment count where forming numerically is the tool of last resort rather than first resort. The same combinatorial test extends unchanged to the tensor networks of Section 15, where the wiring diagram is a reaction network rather than a compartment diagram.

8. Why the fields fragmented

Section 0 opened on three literatures that read as unrelated. Sections 1 through 7 have built one object, at one level of generality, that contains all three; we can now say precisely which fragment of it each field conventionally uses.

Pharmacokinetics uses Section 1's linear specialisation and Section 2's exponential solution almost always, at a small, fixed compartment count. It relies on Section 6's positivity constantly, since a negative concentration is meaningless, but rarely names Metzler matrices or states the constraint as a theorem; the sign pattern is enforced by writing rate constants as manifestly nonnegative symbols rather than by any argument that concentrations stay nonnegative for every input. Controllability (Section 3) is rarely asked at all, since the dosing route is fixed by the clinical protocol rather than designed; observability (Section 4) is asked implicitly, in the form of which compartment can actually be sampled, but rarely by name and never structurally (Section 7).

Pharmacodynamics, as Section 0 noted, is typically not posed as a dynamical system at all. Section 13 shows that its central object, the Hill or Emax dose-response curve, is the output map evaluated at a fixed point of a one-compartment nonlinear flow: the missing differential equation is present all along, solved at equilibrium before pharmacodynamics ever states it.

QSP uses the fully nonlinear generalisation of [eq:state-eq-general], at the large state dimension the other two fields avoid, and needs every tool built in this Part: Metzler positivity generalised to nonlinear mass-action kinetics, and observability generalised from a rank condition on a matrix to a rank condition on Lie derivatives (Section 17). Historically it has built this machinery field by field and model by model, without importing it from control theory as a settled body of results; Part III and Part IV import it.

None of this is a failure of any of the three fields; each solved the problem in front of it. It is an accident of history that the common object underneath was never named. Part II begins naming it, starting from the sign structure Section 6 has already exposed.

Part II: The GENERIC split and the first worked example

9. The reversible/irreversible split

Section 6 isolated a sign pattern in , off-diagonal entries never negative, and showed it alone forces the whole flow nonnegative. That sign pattern does not yet distinguish two processes pharmacokinetics treats as different in kind: drug moving between compartments, and drug leaving the body outright. Non-equilibrium thermodynamics already has an algebra for exactly this distinction, and it applies to unconditionally, before any numbers are chosen.

The extra content in [eq:generic-split-def] is entirely the sign condition . It is not automatic for a compartmental and must be checked case by case rather than assumed; Section 11 checks it for the first worked example below. It is worth pausing on what the two names mean before using them: alone generates a flow that leaves the Euclidean norm exactly unchanged, a pure circulation with no preferred direction of relaxation; alone generates a flow that never increases , pure relaxation toward the origin. Section 10 proves both claims at once.

Resist reading "reversible" here in the everyday pharmacological sense of a route that runs both ways. A compartmental system whose micro-rate-constants happen to be symmetric, for every pair , has already symmetric, so identically and the entire flow is classified as irreversible in the GENERIC sense, even though every route in it runs both ways. This is the correct classification: symmetric exchange is the linear analogue of diffusion, the paradigm case of a process that produces entropy rather than conserving it. instead isolates a net circulation, a persistent asymmetry between the forward and backward rate around some cycle of compartments, and it is the only part of a linear compartmental system with any prospect of not dissipating.

10. BIBO stability as a corollary

Section 5's criterion asks for the spectral abscissa of a restricted matrix , a computation that has to be redone whenever changes. The GENERIC split of Section 9 gives a sufficient condition for the same conclusion that never mentions an eigenvalue directly, phrased instead in terms of dissipation.

For a compartmental system, the connectivity hypothesis has a reading already available from Section 6: it fails only when some subset of compartments exchanges among itself forever without ever reaching a compartment whose column sum in [eq:mass-balance] is a strict inequality, an isolated conservative pocket with no route to elimination. Every compartment diagram in which every compartment has some route, direct or indirect, to elimination satisfies the hypothesis automatically. Section 11 checks the hypotheses of this proposition directly, on numbers, and finds a stronger statement true, strictly positive definite, for which connectivity needs no separate check at all.

11. Worked example, fixed k=2

Every tool built in Part I and specialised in Sections 9 and 10 has so far been stated for an arbitrary , , . This section fixes the smallest case pharmacokinetics actually uses, a central compartment (blood or plasma, the compartment a dose enters and a laboratory samples) exchanging with a single peripheral compartment (a tissue pool with no direct access to either the dose or the assay), and carries every quantity defined so far through on real numbers.

is Metzler ([metzler]): its only off-diagonal entries, and , are both nonnegative. Its columns sum to and ; the first is strictly negative, the elimination rate leaving through the central compartment, and the second is exactly zero, the equality case Section 6 singled out for a compartment that only exchanges and never itself eliminates.

The eigenvalues of follow from its characteristic polynomial, . Here and , giving with discriminant , so

Both eigenvalues are real and strictly negative: alone, with no reference to or , is internally asymptotically stable.

Section 9's split, applied to [eq:pk-k2-matrices], gives

's characteristic polynomial is (, ), with discriminant , so

Both strictly positive: is positive definite, a stronger statement than the that [generic-split] requires and that [bibo-generic] needs, so connectivity is automatic, already. [bibo-generic] reconfirms [eq:pk-k2-eigenvalues]'s conclusion by an entirely different route, through dissipation rather than root-finding, and the two routes agree because they must.

12. Controllability and observability Gramians for the k=2 example

Section 3 and Section 4 stated the Kalman rank tests and the Gramians for an arbitrary and ; this section forms both for [eq:pk-k2-matrices] explicitly and asks what the resulting numbers say about the peripheral compartment nothing directly measures.

The peripheral compartment is never sampled, reads only , and yet says is still uniquely recoverable from alone over any positive time window. The mechanism is visible in [eq:pk-k2-kalman-o] itself: the second entry of is , nonzero exactly because feeds back into , leaving an imprint on the measured trajectory's curvature that a single snapshot of could never carry. Full rank here is the outcome Section 7's structural test already predicts for any nonzero : it is the generic case, holding for every choice of the free rate constants outside the single knife-edge , the one choice that severs 's only route back to a sampled compartment.

13. PD as a gradient flow

Section 0 promised that the pharmacodynamic dose-response curve is a fixed point of a flow rather than an independently fitted object. This section delivers that promise on the simplest receptor-binding mechanism there is: one drug, one receptor, one bound complex.

[eq:receptor-binding-ode] is itself a one-dimensional instance of Section 9's split: a matrix has no antisymmetric part, so identically and carries the whole dynamics. Receptor-binding relaxation is always the fully dissipative endpoint of the same decomposition that gave Section 11's two-compartment flow both a rotational and a contracting part; a single bound complex has nowhere to rotate into, only toward or away from equilibrium.

With , , (so ), [eq:hill-emax-law] reads , . At the half-maximal concentration, , as the name demands; at , above , , strictly between and as [eq:hill-emax-law]'s form requires.

[eq:hill-emax-law] at is a rectangular hyperbola in : saturating, monotonically increasing, and concave throughout, distinct from the S-shaped curve reported for Hill coefficients , which need cooperative multi-site binding that this single-site mechanism does not model. This is the fixed point Section 8 promised: the pharmacodynamics literature's independently fitted four-parameter curve is this equilibrium, arrived at without ever writing down [eq:receptor-binding-ode].

14. First look at scale

Section 13's receptor-binding state was held apart from Section 11's two-compartment flow: entered [eq:receptor-binding-ode] only as an input, never as a coupled state. The moment receptor engagement is tracked as dynamics rather than read off at equilibrium, that separation ends: is one three-coordinate state-space system, , with coupled to exactly as already is. This is the move Section 0's Anchor 2 anticipated: each mechanism added to a model is one more state coordinate, coupled to the ones already built rather than replacing them, and a QSP network is what results from repeating this dozens of times over.

Section 19 carries out this construction with its own rate constants and its own Gramians, in the same style as Sections 11 and 12, and nothing here anticipates its numbers. What can already be said without them is qualitative. stays the only sampled coordinate; , like , leaves an imprint on and so remains observable by the same structural argument Section 12 made for . A binding rate constant need not sit at the same order of magnitude as an inter-compartmental exchange rate, and that mismatch does not show up in a rank computation at all: stays regardless of how weak 's footprint on becomes, while the observability Gramian's smallest eigenvalue shrinks continuously as that footprint weakens. Whether a coordinate is observable and how well it is observable are different questions, and only the first has a rank test. Part IV exists to answer the second, at a compartment count where the answer is no longer visible by inspection.

Part III: QSP as tensor GENERIC

15. From matrix to tensor

Every reaction network built so far has been linear: each rate constant governs an exchange between exactly two coordinates, or one coordinate and the outside world, and the whole flow is generated by acting once with a fixed matrix on the state. A QSP network is rarely built entirely from such reactions. A drug binding a receptor to form a complex, the paradigm reaction for the rest of this Part, has a rate proportional to the product of two concentrations, not to either alone; a rate law of this shape cannot be written as any matrix acting on , however its entries are chosen, because the rate itself is quadratic in the state rather than linear.

Reaction network theory already has the bookkeeping this needs. A network of species undergoing elementary reactions is fixed by its stoichiometric matrix , whose th column records the net change in each species per firing of reaction , together with a rate vector giving each reaction's instantaneous rate as a function of the current state. The dynamics of every mass-action network, whatever the order of its reactions, take the single form . What changes with the order of the reactions is not this equation but the shape of itself.

Nothing about the sign structure Section 6 isolated is special to either. Mass action assigns every rate constant a nonnegative value, and every reaction removes mass from exactly the species that supply it while adding it to exactly the species it produces; entrywise, this forces every entry with to be nonnegative, the tensor generalisation of [metzler]'s off-diagonal sign condition, and the same argument that proved [perron-frobenius-positivity] at carries over unchanged: the flow of [eq:tensor-dynamics] never leaves , at any order.

Section 7's structural test read the sparsity pattern of a linear pair as a directed graph on compartments and sensors, and anticipated that the same reading would apply unchanged to a tensor network; it does, with each nonzero entry of [eq:tensor-dynamics] read as an edge regardless of how many reactant indices it carries. The sparsity figure below previews this at the scale Part IV needs it, before any of that Part's numbers exist.

16. The GENERIC bracket formalism

Section 9's split was built entirely from a single fixed matrix, and nothing in its construction survives once the generator of the flow is an order- tensor rather than a matrix: there is no longer one object to take an antisymmetric part and a symmetric part of. Non-equilibrium thermodynamics already generalises exactly this split to nonlinear, state-dependent flows, under the name Section 9 already borrowed.

The two degeneracy conditions in [eq:generic-degeneracy] are exactly what makes [eq:generic-bracket-evolution] thermodynamically consistent, and the argument is the same one-line collapse [bibo-generic] already used for a single antisymmetric matrix. Taking , is a scalar equal to its own negation ( antisymmetric), hence zero, and by the second degeneracy condition; so , energy is conserved exactly along every trajectory, independently of the particular chosen. Taking , by the first degeneracy condition, using 's antisymmetry to move it across the inner product, and since ; so , entropy never decreases. Neither law of thermodynamics is assumed anywhere in [generic-bracket]: both fall out of antisymmetry, positive semi-definiteness, and the two degeneracy conditions alone.

and generalise and exactly as Section 15's tensor generalises Section 1's matrix: an antisymmetric generator of conservative circulation and a symmetric positive semi-definite generator of monotone relaxation, now allowed to depend on the state and to act through two independent functionals rather than one. The purely linear setting of Sections 9 and 10 is the furthest possible reduction of this structure: a single quadratic functional stands in for two independent functionals and , and constant matrices stand in for state-dependent . With only one functional in play there is no second, independent functional for either degeneracy condition in [eq:generic-degeneracy] to say anything about; both become vacuous rather than substantive, and it is [bibo-generic]'s direct Lyapunov argument, not the bracket algebra above, that does the work there. Nothing about that argument is superseded here: it is the , single-functional case of the same fact.

17. Nonlinear controllability and observability

Sections 3 and 4 built controllability and observability by iterating a single linear map, , and stacking the results into a rank condition. Once the flow is nonlinear in earnest, as Section 15's case forces it to be, there is no linear map left to iterate: "apply again" has to be replaced by an operation that manufactures new directions from nonlinear vector fields, and the operation control theory supplies, developed independently of any pharmacological question, is the Lie bracket.

For smooth vector fields , their Lie bracket is ( the Jacobian matrices of ), the infinitesimal difference between flowing along then and flowing along then : a direction unavailable to either field alone whenever it is nonzero. For a control-affine system (drift , control directions ), the accessibility distribution at is the span, evaluated at , of every vector field obtainable from by repeated Lie bracketing.

[eq:nonlinear-obs-rank] is, entrywise, the tool structural identifiability analysis of a real QSP model actually runs: a state coordinate or parameter direction is structurally unidentifiable exactly when it lies in the common kernel of every , however many derivatives are taken, and software auditing a nonlinear model before it is ever fit to data checks precisely this condition by differential-algebraic elimination rather than by simulation. Section 15's tensor structure is what keeps computable at all for beyond the first two or three: without it, the symbolic expressions in a network of any size explode long before a rank can be read off, which is exactly why Section 22 replaces this symbolic computation with a numerical surrogate once grows large enough.

18. The one geometric remark

Section 15's stoichiometric matrix constrains in a way that has nothing to do with the nonlinearity of ; this constraint has a name in reaction network theory, and stating it once, clearly bounded, is worth doing before returning to ordinary vector-space language for the rest of the article.

19. Worked example extended, k=3

Section 14 promised a third coordinate, coupled rather than held apart, and called it ; this section builds it, gives it its own rate constants, and calls it . The mechanism is target-mediated drug disposition (TMDD)^[5], the standard QSP extension of Section 11's two-compartment model for a drug whose target is abundant enough, or whose affinity for it is high enough, that binding measurably depletes the free-drug pool rather than merely being read off it at equilibrium as Section 13 did. The bound complex is a full third state coordinate, coupled to by mass action rather than by the one-way input of [eq:receptor-binding-ode]: exactly the order-2 term Section 15's tensor formalism was built to carry.

The linearisation is a direct Jacobian computation. Differentiating [eq:qsp-k3-x1] and evaluating at , reduces at to ; everywhere; reduces at to . [eq:qsp-k3-x2] is already linear, contributing the row unchanged by linearisation. Differentiating [eq:qsp-k3-x3], reduces at to , and reduces at to ; every entry of [eq:qsp-k3-matrices] is accounted for. In Section 15's language, the term , appearing with coefficient in [eq:qsp-k3-x1] and in [eq:qsp-k3-x3], is exactly an order-2 rate tensor contribution, coupling two reactant indices rather than one; it is precisely this term that a linearisation, by construction, discards: carries only the order-1 part of the true dynamics, valid near and nowhere claimed to hold far from it.

's off-diagonal entries, (and two structural zeros), are all nonnegative: is Metzler ([metzler]). Its column sums are , , and , exactly as claimed. This is not a coincidence to be checked numerically so much as a prediction to be confirmed: every exchange between and contributes once with each sign to the first two columns, and every binding/unbinding step between and contributes once with each sign to the first and third, so both cancel exactly in the column sums, column 1 collapsing to and column 3 to . Only the two true eliminations, drug clearance and complex internalisation , remove mass from the system at all; mass conservation forces the column sums to equal independently of the specific values of , and the computed matrix confirms exactly this.

's characteristic polynomial has , sum of principal minors , and , giving

All three roots are real and strictly negative: the linearised network is asymptotically stable at the drug-free state, extending Section 11's conclusion to .

The Kalman controllability matrix is built from , the first column of , and , giving

expanding the determinant along the first column, so and is controllable ([kalman-rank-c]). Dually, is built from , the first row of , and , the first row of , giving

expanding along the first row, so and is observable ([kalman-rank-o]). Both determinants nonzero means both Gramians , by the equivalence already established in [kalman-rank-c] and [kalman-rank-o]: Section 14's promise of a third compartment carrying its own controllability and observability certificate is discharged without solving either Lyapunov equation by hand. Section 21 confirms the observability half of this conclusion by an entirely different route, reading it directly off 's sparsity pattern as a graph ([structural-observability]), before any of these numbers are fixed.

Part IV: Compartment overload, observability at scale, and computability

20. The question

Section 19 added a third compartment and it paid for itself: the bound complex carried its own controllability and observability certificate, both Gramians came back positive definite, and Section 14's promise of a coordinate that earns its place was discharged. A real QSP model does not stop at three coordinates; it reaches dozens, as each mechanism a modeller believes in becomes one more state (Section 0's Anchor 2). The question this Part exists to answer is what happens to the certificate as that count grows. Each new compartment either adds real dynamical structure, a direction the single sampled trajectory can actually resolve, or it adds only a direction the data cannot separate from ones already present, inflating the state without adding anything an assay could ever see. We call the second outcome compartment overload, and the rest of Part IV builds the machinery that tells the two apart at a scale where [kalman-rank-o]'s determinant is no longer a thing one reads off by hand, then runs it on a family of networks grown from Section 19's until the overload appears.

21. Structural observability, applied

Before any Gramian is formed, Section 7's combinatorial test can be run on the wiring diagram alone, and it is the correct first screen precisely because it costs no numerics: a network that fails it is unobservable for every choice of rate constants, so there is no point solving a Lyapunov equation for one. We run it on Section 19's network, whose graph is read straight off the sparsity pattern of [eq:qsp-k3-matrices].

Both conditions hold, so the pair is structurally observable: observable for every choice of its rate constants outside a measure-zero exceptional set, the full-rank [eq:qsp-k3-kalman-o] being the generic case rather than a numerical accident. The network passes the cheap screen before any Gramian is formed, which is the outcome a well-posed model should give and a licence to proceed to the finer, value-dependent questions the screen cannot answer. The screen sees only whether a route to the sensor exists at all; it is blind to whether two coordinates take near-identical routes, which is exactly the failure Section 25 engineers and only a Gramian detects.

22. Empirical Gramians

Section 17's Lie-derivative rank test is exact and symbolic, and it stops being computable well before a real QSP network's compartment count: the differentials in [eq:nonlinear-obs-rank] grow in size with each order taken, and the differential-algebraic elimination that reads a rank off them chokes on the resulting expressions once is large and is nonlinear. The tool that survives to that scale replaces symbolic differentiation with simulation, and it reduces to the analytic Gramians of Sections 3 and 4 exactly when the dynamics happen to be linear.

This construction, due to Lall, Marsden, and Glavaski^[6], is what a structural identifiability analysis of a large nonlinear QSP model runs in place of Section 17's symbolic test once that test is out of reach. It asks only that the model can be simulated, never that its Lie derivatives can be written down, and it returns a pair of matrices in the same coordinates the analytic Gramians live in, so every question the next two sections ask of can be asked of without changing a line. Because our Section 25 family is linearised, its coincide with the analytic Gramians, and the numbers reported there are computed from the Lyapunov equations directly.

23. Balanced truncation and Hankel singular values

A positive-definite certifies that every initial state is recoverable, and a positive-definite that every state is reachable, but neither number-free verdict says how well. Section 14 already isolated the gap: stays full while the smallest eigenvalue of shrinks continuously as a coordinate's imprint on the sensor weakens. The quantity that grades a coordinate on this continuous scale, and does so in a way no change of state coordinates can distort, combines both Gramians at once.

Each grades one input-output direction by how strongly it is jointly reachable from the dose and visible at the sensor: a direction with large is excited hard by the input and imprints hard on the output, and a direction with near zero is one the dose barely excites, or that barely reaches the assay, or both. This grading has an operational payoff that names the whole method.

[eq:truncation-bound] is what turns the spectrum [eq:hankel-sv] into a verdict. When the Hankel singular values fall off a cliff, a sharp drop of one or more orders of magnitude between and , the tail beyond the cliff is negligible and the model has an effective order well below its nominal : the extra coordinates are carrying a joint reachability-observability weight the data cannot resolve. This is the signature of a sloppy model^[8], and reading it off the spectrum is the quantitative form of the overload question Section 20 posed. The location of the cliff is the number of directions the single sampled trajectory can actually resolve; everything past it is compartment overload measured in decibels.

24. Computability

Forming and the way Section 12 did, by solving the Lyapunov equations [eq:lyapunov-c] and [eq:lyapunov-o] densely, costs arithmetic and storage, the same as one dense eigendecomposition. At the of this article that cost is nothing; at the compartment counts a systems-scale QSP model reaches it is the difference between a computation that runs and one that does not, and the sparsity Section 6, Section 7, and Section 15 established is what a scalable solver exploits.

This closes the loop the article has been building toward. The sparsity that Section 6's positivity and Section 7's graph test read as structure is the same sparsity a low-rank Lyapunov solver reads as a preconditioner; the low-rank factor it returns is the same object [balanced-truncation] truncates and Section 20's question interrogates. The cheap structural screen, the affordable Gramian, and the reduced model that reports the overload are three readings of one sparse operator.

25. Scaling family

We now grow Section 19's network in two deliberate ways and watch the Hankel spectrum. First to , adding two near-duplicate peripheral compartments and one slow deep compartment alongside the retained TMDD complex; then to , adding a third near-duplicate peripheral, a second deep compartment, and extending the complex into a two-step cascade that ends in an internalised product. The rates are chosen so that both matrices stay Metzler and mass-conserving, so nothing about their positivity or stability is at issue; the only variable under test is what the Hankel spectrum does as the count climbs.

26. Verdict and synthesis

The spectrum [eq:scaling-hsv8] lost a direction outright, and [balanced-truncation]'s balancing coordinates say which. We form for [eq:scaling-A8], take the eigenvector of belonging to its smallest (zero) eigenvalue, and read off the state coordinates that dominate it.

The two failures are the two answers Section 20's question admits. The internalised product is unobservable exactly, a coordinate no assay on can ever recover, caught by a graph test that costs no arithmetic. The near-duplicate peripherals are unobservable only in practice, a direction whose Hankel weight has fallen five orders of magnitude below the leading one, caught only by grading directions on the continuous scale [hankel-singular-values] supplies. A modeller who added either coordinate believing it bought new structure added instead a direction the single sampled trajectory cannot resolve, and the machinery of this Part is what returns that verdict before a single parameter is fit.

That machinery was never pharmacological. Part I fixed the state-space object and its controllability, observability, and stability once, for arbitrary compartment count. Part II read the mass-action sign pattern as a Metzler structure and the reversible/irreversible split as a decomposition. Part III lifted the generator from a matrix to a rate tensor and the split to the GENERIC bracket, carrying every Part I tool across the lift unchanged. Part IV then spent that generality: the structural screen, the empirical Gramian, the Hankel spectrum, and the low-rank solver are one body of control theory, and running them on a QSP network is importing settled results rather than rebuilding them mechanism by mechanism. The three literatures Section 0 opened on were one object all along; naming it is what let a question about compartment overload be answered by a theorem about Gramians. The same object sits under any population of such models fit across many individuals at once, an ensemble extension we note here only as an open direction.

References

1. R. E. Kalman (1963). Mathematical description of linear dynamical systems. Journal of the Society for Industrial and Applied Mathematics Series A: Control, 1(2), 152-192. DOI
2. C.-T. Lin (1974). Structural controllability. IEEE Transactions on Automatic Control, 19(3), 201-208. DOI
3. H. C. Öttinger (2005). Beyond Equilibrium Thermodynamics. Wiley-Interscience. DOI
4. M. Feinberg (2019). Foundations of Chemical Reaction Network Theory. Springer. DOI
5. D. E. Mager, W. J. Jusko (2001). General pharmacokinetic model for drugs exhibiting target-mediated drug disposition. Journal of Pharmacokinetics and Pharmacodynamics, 28(6), 507-532. DOI
6. S. Lall, J. E. Marsden, S. Glavaški (2002). A subspace approach to balanced truncation for model reduction of nonlinear control systems. International Journal of Robust and Nonlinear Control, 12(6), 519-535. DOI
7. B. Moore (1981). Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE Transactions on Automatic Control, 26(1), 17-32. DOI
8. R. N. Gutenkunst, J. J. Waterfall, F. P. Casey, K. S. Brown, C. R. Myers, J. P. Sethna (2007). Universally sloppy parameter sensitivities in systems biology models. PLoS Computational Biology, 3(10), e189. DOI