Mathematical Pharmacology

Introduction

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 is the case in which the response curve is the output map evaluated at a fixed point of a compartmental flow, linear or Michaelis-Menten regardless of its compartment count; Section II.6 works the simplest, one-compartment instance. QSP is the case in which and are left fully nonlinear on the pharmacokinetic and pharmacodynamic sides at once, 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 II.1.

Two anchors run through everything that follows.

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.

I.1. The state-space object

The three literatures the introduction named 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 concentration of drug in compartment , an absolute amount (mass or moles) divided by a constant volume of distribution , the passage Section I.6 derives from conservation of mass; 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 I.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 I.4 and Section I.5 need to say precisely what it does and does not affect. Sections I.2 through I.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.

I.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 I.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 I.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 II.4.

I.3. Controllability

The free response of Section I.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]

I.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 I.3's theorem transposes without change.

I.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 I.7 and Part IV return to under the name of Anchor 2, observability at scale. Section II.3 shows this criterion collapses to a single positivity check once the reversible/dissipative split of Section II.2 is in place.

I.6. Metzler matrices and positivity

Compartmental kinetics 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, derives it from nothing more than conservation of mass, and shows it forces positivity on the whole flow, for every .

Column of collects every rate at which compartment loses concentration: to other compartments (the off-diagonal entries, each nonnegative by [eq:metzler-def]) and out of the system entirely (elimination or clearance, the excess by which [eq:mass-balance] falls short of equality). Equality in [eq:mass-balance] for a given column means compartment only exchanges with other compartments and is never itself a route out of the system; this is exactly the reversible, non-eliminating structure Section II.2 isolates algebraically.

I.7. Structural (generic) observability

The rank test of Section I.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 IV.2 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 III.1, where the wiring diagram is a reaction network rather than a compartment diagram.

II.1. Why the fields fragmented

The introduction opened on three literatures that read as unrelated. Sections I.1 through I.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 I.1's linear specialisation and Section I.2's exponential solution almost always, at a small, fixed compartment count. It relies on Section I.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 I.3) is rarely asked at all, since the dosing route is fixed by the clinical protocol rather than designed; observability (Section I.4) is asked implicitly, in the form of which compartment can actually be sampled, but rarely by name and never structurally (Section I.7).

Pharmacodynamics, as the introduction noted, is typically not posed as a dynamical system at all. Section II.6 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 III.3). 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 I.6 has already exposed.

II.2. The reversible/irreversible split

Section I.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 II.4 checks it for the first worked example below. The two names carry distinct content: 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 towards the origin. Section II.3 proves both claims at once.

The reversible part has a name in another language too. Antisymmetric matrices are exactly the generators of rotations: under the standard identification of antisymmetric matrices with bivectors (grade-2 elements of a geometric algebra, each a space of dimension ), corresponds to a bivector , and is exactly the rotor sandwich for , in the sense this site's geometric-algebra treatment of classical mechanics develops in full for rigid-body rotation. Section II.4's is the simplest possible instance: a single bivector generating planar rotation at rate , with exactly the ordinary rotation matrix at that rate.

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.

II.3. BIBO stability as a corollary

Section I.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 II.2 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 I.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 II.4 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.

II.4. Worked example, fixed k=2

Every tool built in Part I and specialised in Sections II.2 and II.3 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 I.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 II.2'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.

II.5. Controllability and observability Gramians for the k=2 example

Section I.3 and Section I.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 I.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.

II.6. PD as a gradient flow

The introduction 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 II.2'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 II.4's two-compartment flow both a rotational and a contracting part; a single bound complex has nowhere to rotate into, only towards 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 II.1 promised: the pharmacodynamics literature's independently fitted four-parameter curve is this equilibrium, arrived at without ever writing down [eq:receptor-binding-ode].

II.7. First look at scale

Section II.6's receptor-binding state was held apart from Section II.4'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 the introduction'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 III.5 carries out this construction with its own rate constants and its own Gramians, in the same style as Sections II.4 and II.5, 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 II.5 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.

III.1. 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 . The order of the reactions changes only the shape of ; this equation stays fixed.

Nothing about the sign structure Section I.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 I.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.

III.2. The GENERIC bracket formalism

Section II.2'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 II.2 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 III.1's tensor generalises Section I.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 II.2 and II.3 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.

III.3. Nonlinear controllability and observability

Sections I.3 and I.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 III.1'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 III.1'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 why Section IV.3 replaces this symbolic computation with a numerical surrogate once grows large enough.

III.4. The one geometric remark

Section III.1'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 it is stated once, clearly bounded, before returning to ordinary vector-space language for the rest of the article.

III.5. Worked example extended, k=3

Section II.7 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)^[7], the standard QSP extension of Section II.4'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 II.6 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 III.1'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 III.1'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 II.4's conclusion to .

Section II.2's split, applied to the linearisation [eq:qsp-k3-matrices] exactly as it was to [eq:pk-k2-matrices] in Section II.4, gives

's characteristic polynomial is (, the sum of its principal minors , ), with roots

All three strictly positive: is positive definite, extending Section II.4's conclusion to and reconfirming [eq:qsp-k3-eigenvalues]'s stability verdict by the dissipative route [bibo-generic] supplies, independent of the root-finding above.

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 II.7's promise of a third compartment carrying its own controllability and observability certificate is discharged without solving either Lyapunov equation by hand. Section IV.2 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.

IV.1. The question

Section III.5 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 II.7'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 (the introduction'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 III.5's until the overload appears.

IV.2. Structural observability, applied

Before any Gramian is formed, Section I.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 III.5'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 IV.6 engineers and only a Gramian detects.

IV.3. Empirical Gramians

Section III.3'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 I.3 and I.4 exactly when the dynamics happen to be linear.

This construction, due to Lall, Marsden, and Glavaski^[8], is what a structural identifiability analysis of a large nonlinear QSP model runs in place of Section III.3'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 IV.6 family is linearised, its coincide with the analytic Gramians, and the numbers reported there are computed from the Lyapunov equations directly.

IV.4. 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 II.7 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, and it names the 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^[11], and reading it off the spectrum is the quantitative form of the overload question Section IV.1 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.

IV.5. Computability

Forming and the way Section II.5 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 I.6, Section I.7, and Section III.1 established is what a scalable solver exploits.

This closes the loop the article has been building towards. The sparsity that Section I.6's positivity and Section I.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 IV.1'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.

IV.6. Scaling family

We now grow Section III.5'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.

IV.7. Verdict

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 IV.1'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.

Every object built so far holds the numbers fixed: a matrix whose entries are settled rate constants, a Gramian assembled from them, a Hankel spectrum read off once those numbers are chosen. This Part lets the numbers vary. We collect the rate constants into a parameter vector and ask what varying it buys. Section V.1 builds a Gramian for parameters exactly as Part I built one for states, and shows the observability Gramian of Section I.4 is the single instance in which the parameter is the state itself. Section V.2 reads that parameter Gramian as the Fisher information of a sampled, noisy experiment and imports its Riemannian meaning from this site's probability and statistics article. Section V.3 separates the failure the numbers-free structural test of Section I.7 already catches from the failure only the spectrum sees, and closes the sloppy-model loop Section IV.4 left open. Sections V.4 through V.7 then fit many such systems at once, one per subject, and ask what changes when a population of trajectories replaces a single one.

V.1. Sensitivity equations and the parameter Gramian

Every theorem of Part I holds the system matrices fixed and varies the input or the initial state. We now hold a trajectory fixed and vary the matrices, letting depend on a parameter vector that collects the rate constants a modeller estimates. The dosing matrix and the sensor matrix stay fixed throughout this article's worked examples, since the route of administration and the choice of assay are settled by the experiment rather than fitted. Differentiating the state equation in each parameter produces a linear system with the same dynamics as the original, so Section I.2's machinery transfers without a new idea.

Stack the sensitivities into a single matrix , whose th column records how the whole state trajectory moves when the th parameter is nudged. The sensitivity of the output is then , and the weight each parameter direction carries at the sensor, integrated over the experiment, is a Gramian built from exactly as Section I.4's was built from .

The construction has cost nothing new: it is Section I.4's observability Gramian with replaced by the sensitivity matrix . The replacement is exact, and one parametrisation makes the two objects coincide.

The proof is a single line. Under the solution [eq:variation-of-parameters-eq] reads , and only the first term carries , so column by column , giving ; substituting into [eq:parameter-gramian-def] is [eq:parameter-gramian-is-Wo-eq].

The rate-constant parametrisation that Section V.2 builds its Fisher information on runs the same construction with and the forcing term of [eq:sensitivity-solution] switched back on. It is the general [parameter-gramian], read at a different parametrisation of the same object, rather than a new device introduced for parameters. Everything Part I proved about , its Lyapunov characterisation, its kernel, its place in the Hankel spectrum, is a statement about [eq:parameter-gramian-def] evaluated at .

V.2. The Fisher information matrix as a Riemannian metric

Section V.1's parameter Gramian integrates over a continuous horizon. A real experiment samples the output at finitely many times and reads each sample through measurement noise. Restoring both features turns the parameter Gramian into the Fisher information matrix, and the statistical geometry this site's probability and statistics article built for the Fisher information transfers to it wholesale.

Add a Gaussian observation model: the output is sampled at times and each sample is read through independent noise,

This article's worked examples all read a single compartment, , so is a scalar variance; the vector-output case replaces by the inverse noise covariance at every occurrence, with no other change.

The discrete sum, rather than [eq:parameter-gramian-def]'s integral, is the object we want here for two reasons. A continuous-time integral would carry an extra factor of that spoils the dimensional statement of the remark below; and a real experiment samples at finitely many times, so [eq:fisher-information-def] is the quantity a population-pharmacokinetic optimal-design calculation actually forms.

V.3. Practical vs. structural identifiability

Section V.1's Gramian and Section V.2's metric grade parameters on a continuous scale. Two failures of that grading separate cleanly, because one is caught by the numbers-free graph test of Section I.7 and the other only by the spectrum. The distinction is the parameter-space form of the state-space distinction Section IV.4 already drew.

The condition is necessary only. Its converse fails: a structurally observable system can still carry a Fisher matrix that is nonsingular yet arbitrarily close to singular, and the Cramer-Rao bound (the statistics article's Section 2.5) reads the consequence at once: the covariance of any unbiased estimator of is bounded below by , so a parameter combination lying near has asymptotic variance growing without bound as its Fisher eigenvalue shrinks, however finely the output is sampled. That is practical unidentifiability, the sloppy-model regime of Gutenkunst et al.^[11], and it is a statement about the numbers in that no graph test can see. Structural failure is a knife-edge of the wiring; practical failure is a matter of degree in the spectrum.

The empirical spectra bear the split out. We compute [eq:fisher-information-def] for the rate-constant parametrisations of the networks of [qsp-k3] and [scaling-family], at a single bolus into and a single assay on , and read off the eigenvalues. Every spectrum falls off a pronounced cliff, with one or more eigenvalues collapsing to the numerical floor of the Lyapunov solve, parameter combinations no single-dose, single-sensor experiment can resolve at any measurement precision, the direct hits of [structural-forces-singular-fim]; several further eigenvalues sit orders of magnitude below the leading one, the sloppy directions the Cramer-Rao bound penalises with large but finite variance. The number of floor-level eigenvalues grows as the parameter count grows from to , and the qualitative shape is the rate-constant counterpart of the Hankel cliff of [scaling-family], now grading parameters rather than states.

V.4. The population extension

Everything in Part V so far fits one subject's data to one parameter vector. A population study measures many subjects at once, each metabolising the same drug through the same compartmental wiring but with rate constants of their own, and asks for the distribution of those rate constants across the population rather than any single subject's values. We lift [state-space] to a hierarchy: a population-typical parameter shared by everyone, and a subject-specific deviation drawn around it.

The log-scale placement of the random effect in [eq:population-hierarchy] is forced by Section I.6. Every entry of is a rate constant, and [metzler] together with [perron-frobenius-positivity] requires every one of them to be positive for the flow to stay entrywise nonnegative, so that no subject's model ever predicts a negative concentration. The multiplicative form holds for every real , so the Gaussian random effect on the log scale keeps each rate constant inside the positive cone for free, whatever value it takes. An additive random effect would place positive probability on , violating the Metzler sign pattern of [eq:metzler-def] and forfeiting the positivity guarantee that holds for every compartment count . The positivity structure of Section I.6 is what selects the parametrisation of the population model.

The likelihood of the population parameters marginalises each subject's random effect against its prior,

with the Gaussian data density read off [eq:population-obs] and the density. This integral has no closed form because depends nonlinearly on through the exponential of [eq:population-hierarchy] feeding , even though [eq:state-eq-linear] is linear in : the intractability sits entirely in how enters the system matrix, and none of it in the propagation of the state.

V.5. FOCE and SAEM as approximations to the population likelihood

The marginal likelihood [eq:population-likelihood] cannot be evaluated in closed form, so every population-pharmacokinetic estimator is a way of approximating the same -dimensional integral. The two in routine use, FOCE and SAEM, are the two natural choices: linearise the integrand and integrate it exactly, or leave it alone and sample from it. The first is where Section V.2's Fisher information reappears, since the curvature a Laplace approximation needs at each subject's mode is exactly that Fisher information pulled back to the random-effect coordinate.

The identity [eq:foce-hessian] is the payoff of Part V's framing. The quantity a population fitting routine forms at every iteration, the curvature of each subject's conditional log-likelihood at its mode, is the parameter Gramian of Section V.1 sampled and noise-weighted into the Fisher information of Section V.2, pulled to the log scale by the same transformation law [eq:fim-transform] that made that Fisher information a metric in the first place. The random-effect prior adds the single term , the only ingredient the single-subject theory of Section V.2 lacked.

V.6. Software instantiation: mrgsolve, nlmixr2, and typed units

The two directions of Part V each have a standard instantiation in the population-pharmacokinetics software stack, one package per arrow. The forward map, from a population triple to simulated data, is what Section V.4's hierarchy describes; the inverse map, from data back to the triple, is what Section V.5's algorithms compute.

mrgsolve runs the forward direction. It compiles a user's ODE model specification to a shared library carrying an LSODA integrator, with a C++ backend, and simulates [population-state-space] at whatever compartment count the model declares: fix a population triple and a dosing design, draw each subject's random effect by [eq:population-hierarchy], form each , and integrate [eq:population-obs] forward for every subject.^[14] It is the single-subject object of [state-space] from Part I, replicated across a population and pushed forward in time, evaluated exactly as [eq:population-obs] is written.

nlmixr2 runs the inverse direction. The succession runs as follows: the original nlmixr package is archived, and its maintained successor nlmixr2 carries the estimation forward, built on the compiled-ODE engine rxode2 and adding the nlmixr2est estimation layer.^[15] That layer implements the FOCE, FOCEI, and SAEM approximations of [foce-laplace] and [saem], solving the inverse problem of the marginal likelihood [eq:population-likelihood]: recover from the observed subject data . Where mrgsolve maps a triple to data, nlmixr2 maps data to a triple, along the same population hierarchy read in opposite directions.

The Rust uom crate carries physical dimensions in the type system by way of typenum, so a dimensionally inconsistent combination is rejected before the program runs. Scaling a concentration by a dimensionless number stays a concentration and compiles; adding a concentration to a time does not typecheck, because the two inhabit distinct dimensional types.

use uom::si::f64::{Time, MolarConcentration};
use uom::si::time::hour;
use uom::si::molar_concentration::millimole_per_liter;

fn main() {
    let t = Time::new::<hour>(2.0);
    let conc = MolarConcentration::new::<millimole_per_liter>(0.5);

    // Compiles: scaling a concentration by a dimensionless factor
    // is still a concentration.
    let scaled: MolarConcentration = conc * 0.5;

    // Rejected at compile time: concentration and time are distinct
    // dimensional types, so `uom` refuses their sum, the same tensorial
    // mismatch that `fim-is-metric` forbids at runtime.
    // let bad = conc + t;

    let _ = (t, scaled);
}

The compile error uom raises for conc + t is the same check that [fim-is-metric] makes on 's units , made mechanical: an instance of the abstract theorem realised in the type system, the same rhythm as the worked examples of Sections II.4, III.5, and IV.6, where a general result was read off a specific network.

V.7. Synthesis

One object has carried the whole article. 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 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 imported settled results rather than rebuilding them mechanism by mechanism. Part V then let the numbers vary, and the same object answered a second family of questions with no new construction.

The parameter derivatives of that object are its Fisher information. Differentiating the state trajectory in the rate constants (Section V.1) builds a Gramian for parameters exactly as Part I built one for states, and recovers the observability Gramian of Section I.4 as the single instance in which the parameter under study is the initial state itself. Sampling that Gramian through a noisy assay (Section V.2) turns it into the Fisher information matrix, a -tensor on parameter space whose Riemannian meaning is the Fisher-Rao metric of this site's probability and statistics article, imported wholesale rather than rebuilt.

Structural observability (Section I.7) is necessary but not sufficient for the practical identifiability that machinery grades on a continuous scale. The numbers-free graph test catches a coordinate that no assay on can ever recover; the Cramer-Rao bound, read off , penalises the directions resolvable in principle and sloppy in practice, with variance growing as their Fisher eigenvalue shrinks (Section V.3). That closes the loop Section IV.4 opened: the Hankel cliff there, cited to Gutenkunst et al.^[11] as an analogy, is by [hankel-graded-identifiability] a cliff in Fisher eigenvalues, and the sloppy-model reference is now a proved identity in place of a gesture.

The population layer is that same object copied across subjects, the extension the article's earlier close marked as unfinished business; Part V is exactly that extension, developed in full. The hierarchy of Section V.4 replicates [state-space] once per subject, its log-scale random effect forced by the positivity of Section I.6 rather than chosen for convenience; the estimators that fit it, FOCE and SAEM (Section V.5), are a Laplace approximation and a stochastic expectation-maximisation assembled from the sensitivity-and-Fisher machinery of Sections V.1 and V.2, their per-subject curvature being that Fisher information pulled to the random-effect coordinate by the same transformation law that made it a metric. The software running these two arrows, mrgsolve forward and nlmixr2 back (Section V.6), simulates and estimates the single-subject object of [state-space] pushed across a population and read in either direction, rather than any device new to population work.

One discipline runs under all of it. Correctness at every level of this construction, from a compartment's sign pattern (Section I.6) to a population parameter's asymptotic variance (Section V.3), is a statement about which transformations an object is invariant under. The positivity of Section I.6 is invariance of the state cone under the flow; the coordinate-independence of Section IV.4 is invariance of the Hankel spectrum under a change of state basis; the tensoriality of the Fisher information ([fim-is-metric]) is invariance of the metric under reparametrisation. A type system that tracks physical dimension (Section V.6) is a computational instance of that same discipline: refusing the sum of a concentration and a time before the program runs enforces, at the level of arithmetic, the one invariance the article has been reading off the object at every scale.

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. R. Hermann, A. J. Krener (1977). Nonlinear controllability and observability. IEEE Transactions on Automatic Control, 22(5), 728-740. DOI
5. H. J. Sussmann, V. Jurdjevic (1972). Controllability of nonlinear systems. Journal of Differential Equations, 12(1), 95-116. DOI
6. M. Feinberg (2019). Foundations of Chemical Reaction Network Theory. Springer. DOI
7. 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
8. 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
9. B. Moore (1981). Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE Transactions on Automatic Control, 26(1), 17-32. DOI
10. K. Glover (1984). All optimal Hankel-norm approximations of linear multivariable systems and their L-infinity error bounds. International Journal of Control, 39(6), 1115-1193. DOI
11. 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
12. M. J. Lindstrom, D. M. Bates (1990). Nonlinear mixed effects models for repeated measures data. Biometrics, 46(3), 673-687. DOI
13. E. Kuhn, M. Lavielle (2005). Maximum likelihood estimation in nonlinear mixed effects models. Computational Statistics & Data Analysis, 49(4), 1020-1038. DOI
14. K. T. Baron (2026). mrgsolve: Simulate from ODE-Based Models. CRAN. Link
15. M. Fidler, J. J. Wilkins, R. Hooijmaijers, T. M. Post, R. Schoemaker, M. N. Trame, Y. Xiong, W. Wang (2019). Nonlinear mixed-effects model development and simulation using nlmixr and related R open-source packages. CPT: Pharmacometrics & Systems Pharmacology, 8(9), 621-633. DOI