Structure is sparsity

A biochemical model’s structure is not a different kind of information from its parameters. It is the support of the parameter vector.

Fangzhou Xiao · Westlake University · a working paper, August 2026. Companion to Three schools, one question and The biomachine perspective.

In 2001 James Bailey published a two-page commentary called “Complex biology with no parameters”. It made a claim that the next twenty-five years of computational biology were built on. The information you need to predict what a cell will do splits in two. There is systems-structural information, meaning which components are in the system and which of them interact. And there is parametric information, meaning the rate constants and the binding constants. Bailey’s argument was that the first half is nearly free, because the genome sequence gives it to you, and that the second half is unobtainable and, for most questions worth asking, unnecessary.

The first half of that claim is testable, and Bailey tested it himself, in the same paragraph in which he made it. He noted that the 1999 attempt to build a genome-scale Escherichia coli network from the full genome sequence could not reproduce the organism’s already-known physiology without the addition of 47 chemical reactions to the network based on biochemical information, which was not derived from genome sequence. He filed this as a difficulty in implementation. This document argues it is the whole story.

Two things follow, and they point in opposite directions. The first is empirical: the genome does not give you the network, it gives you a first draft that has to be repaired against biochemistry, and the repair does not converge. The second is conceptual, and it is the reason the first one is not a scandal. Structure and parameters are the same object. A reaction is in a model exactly when its rate constant is nonzero, so the structure of a model is the support of its parameter vector, and the split Bailey drew is a split between the zeros and the nonzeros of one list of numbers. Everything that is true of sparse estimation is therefore true of network discovery, including the parts nobody wants: the support is not identifiable, the threshold that defines it is set by measurement noise rather than by chemistry, and it moves when the cell changes.

1Two kinds of information

Bailey's split is the organising idea of computational cell biology. It is worth stating in his own words before arguing with it.

The question Bailey asked is prior to any particular model. Suppose you want to predict something about a cell. What do you need to know? His answer was that the required knowledge divides in two.

"The information required to predict or calculate such characteristics can also be divided into two categories: first, mechanistic or systems-structural (here we use the term “systems structure” to refer to a particular set of components that are involved in the system, and the network of interactions or connections among these components); and second, quantitative or parametric (e.g., what is the equilibrium constant for reversible association of two proteins, or what equation reasonably describes the rate of an enzyme-catalyzed reaction, and what numerical values should be employed for the parameters in such a kinetic formula?)." James E. Bailey, Complex biology with no parameters, Nature Biotechnology 19:503, 2001¹. It was his last paper; he died on 9 May 2001 while it was in press.

Two claims follow, and they are separable. The first is about supply.

"In principle, genome sequence may provide all of the information necessary to define the systems structure of a biological system of interest. For example, knowing all of the enzymes in a cell, as well as all of the substrates that each one accepts and all of the products that each one can make, it is possible to formulate a master global reaction network that represents the complete repertoire of possible biochemical reaction systems within that cell." Bailey 2001¹

The second is about demand. Bailey argued that you do not need the parametric half for most questions, because chemical engineering has a century of results that extract qualitative answers from structure alone: which steady states exist, whether a system can oscillate, whether it can switch. He cites Wei and Prater on reaction networks, Horn and Jackson on general mass action, Feinberg and Horn on deficiency²¹. These theorems take a reaction network and return conclusions that hold for every choice of rate constants. If you have the network, you get the conclusion for free.

Put together, the argument is an argument about where to spend effort. Sequencing was becoming cheap and rate constants were not, so the strategy that scales is the one that leans on the half that is arriving in bulk. It was the right call in 2001, and it was continuous with the discipline Bailey had named a decade earlier². The constraint-based programme that acted on it produced two decades of results that no kinetic programme came close to matching¹³.

One fact belongs on the record before the argument starts, because it changes who the argument is with. Bailey did not reach for this dichotomy in the absence of an alternative. In 1996 he co-authored a formulation in which the presence or absence of each regulatory interaction is an integer variable, optimised jointly with the continuous strengths, over a candidate set the paper calls a regulatory superstructure⁴³. So the two categories of 2001 are a position taken after building the machine that treats them as one. §10 returns to this.

2Forty-seven reactions

The counterexample to the supply claim is in the sentence immediately after it.

Bailey did not assert that the genome gives you the network and stop. He checked, and he reported what he found.

"However, reaching this state of enlightenment remains a goal for the future; it was not possible to find all of the enzymes involved in even amino acid biosynthesis in Escherichia coli in 1998, and efforts to formulate such a global master reaction network for E. coli from the full genome sequence in 1999 could not describe the already known physiological and metabolic properties of E. coli without the addition of 47 chemical reactions to the network based on biochemical information, which was not derived from genome sequence." Bailey 2001¹, citing Bono et al. 1998 and Edwards & Palsson 2000³

He filed this under difficulties of implementation. Read it instead as a measurement, because it is one. The reference is to the first genome-scale metabolic reconstruction of E. coli³, and that paper prints its own provenance. Its Table 1 lists every gene in the model, by pathway category, and marks with a citation each entry that came from somewhere other than the annotated genome.

Parsing that table gives the numbers directly. There are 706 entries, covering 671 distinct gene symbols, since a gene serving two categories is listed twice. Of those entries, 46 carry a literature citation rather than resting on the genome annotation. Bailey counted 47 reactions, and entries are not reactions, so the two counts are the same measurement with different units.

The more telling number is inside the 46. Sixteen of them are not gene symbols at all.

These sixteen are the sharp end of the point. The other thirty are genes that the genome annotation had missed and a paper had found, which is a gap the annotation could in principle close. The sixteen are different. They are reactions with no gene, and no amount of better annotation of the sequence you have will produce one, because the thing you are looking for is not in the sequence in a form the search can see.

3How a reconstruction is actually built

The procedure is published, in detail, as a protocol. Reading it settles the question of whether the network comes from the sequence.

Thiele and Palsson wrote the method down as a Nature Protocols paper: ninety-six numbered steps, in four stages, with a stated duration of several months to a year for one organism⁴. Stage one is the draft reconstruction, which is the part a genome annotation supplies. Stages two through four are refinement, conversion to a computable form, and network evaluation, and they are where the work is.

The step that matters here is called gap-filling, and it has its own literature. The situation it addresses is this. A draft network built from annotation typically cannot produce biomass, or cannot grow on a carbon source the organism demonstrably grows on. Something is missing, and nothing in the sequence says what. So you search a universal reaction database, drawn from biochemistry across all organisms, for a small set of reactions whose addition restores the observed phenotype.

Reed and colleagues ran the loop end to end and reported what it produced⁵. The model disagreed with growth phenotyping data. The algorithm proposed missing reactions. Homology and context methods proposed genes that might carry them. Experiments then tested the proposals, verifying five cases and assigning function to eight ORFs: yjjLMN, yeaTU, dctA, idnT and putP, comprising two new enzymatic activities and four transport functions.

Notice the direction of that arrow. The genome did not supply the network. The network, plus phenotype data, supplied the genome annotation. Three of those eight ORFs begin with y, the prefix E. coli genetics uses for a gene of unknown function, which is the point: they were unannotated precisely because sequence alone had not been able to say what they did.

4Why the gap does not close

The gap is not a backlog. Three mechanisms keep it open, and all three are properties of enzymes rather than of annotation software.

It would be easy to read the last two sections as a story about 1999 being early. It is not, and the reason is that the inference the annotation performs has a known failure mode in both directions.

Annotation assigns function by transferring it from a characterised homologue. That inference is sound when sequence similarity implies functional identity. It fails when a sequence has an activity its homologues do not, and it fails when an activity has no characterised homologue at all.

Promiscuity breaks the transfer from the sequence side. Enzymes commonly catalyse reactions other than the one they are named for, at rates that are low in absolute terms and entirely sufficient to matter physiologically. Khersonsky and Tawfik's review is the standard statement²⁰. The consequence for a reconstruction is that the true reaction set is strictly larger than the annotated one, by an amount no sequence comparison reports.

This is not hypothetical for E. coli. Guzmán and colleagues used a model that included promiscuous activities to predict which of them would become physiologically load-bearing when a canonical gene was deleted, and confirmed the predictions experimentally¹⁸. The activities were there all along, below the threshold at which anyone had called them reactions. Nam and colleagues had already shown that this is the normal condition rather than an exception: enzyme specificity is a matter of degree, set by network context and selection, not a property the sequence declares¹⁹.

Orphan activities break the transfer from the function side. An orphan is a reaction that biochemists have measured and no one has connected to a gene. The sixteen entries in Table 1 are orphans, and orphans are exactly the reactions a sequence-based method has nothing to search with: there is no query. The only route to closing one is the loop §3 described, where the network proposes the activity and experiments hunt the gene, which is how yjjLMN and yeaTU acquired functions⁵.

Osterman and Overbeek named the problem and gave it its shape¹⁶. Their observation is that the two failure directions are different problems with different remedies. A gene with no function is a hypothetical protein, and there are always many: unassigned genes are 20 to 60% of the proteins in most sequenced genomes. A function with no gene is a missing gene, and the asymmetry is that you cannot even pose the second question without a reconstruction first telling you the function should be there. Their sentence for this is that a reconstruction “provides a rather specific and precise notion of what is actually missing.” The reconstruction is not merely incomplete. It is the instrument that measures its own incompleteness.

They also split the missing genes in two. A locally missing gene has a known form in some other lineage and an unknown one here, which is non-orthologous gene displacement and is what comparative genomics is good at. A globally missing gene has no representative sequence in any organism, so no amount of cross-genome comparison will produce a candidate. Sizing this, they put the shared “central machinery” of life at roughly four thousand functional roles, of which any one organism uses 300 to 3000, and estimated that about 10% remain globally missing. For 2003 that is close to what the reconstruction of the same year shows independently: 6.2% of the internal reactions of iJR904 carry no gene.

And a third of the genome has never been characterised. Ghatak and colleagues built a reproducible workflow for asking which E. coli genes have experimental evidence of function, as opposed to a name inherited by homology. The answer names itself.

5The number

The question is whether the gene-less fraction of a genome-scale network shrinks toward zero as curation improves. It is answerable, because every model in the lineage is public.

The E. coli reconstructions form an unusually clean series: iJR904⁷, iAF1260⁸, iJO1366⁹ and iML1515¹⁰, each superseding the last. Each is a named, versioned, downloadable model, and each reaction in each carries a gene-protein-reaction rule saying which genes are responsible for it, or carries an empty rule saying that no gene is. Counting the empty rules across the series measures exactly the thing Bailey's supply claim is about.

The count below excludes exchange, demand and sink pseudo-reactions, which are boundary bookkeeping rather than chemistry and were never expected to carry genes. What is left is reactions inside the cell for which the reconstruction asserts a chemical transformation and no gene.

Three readings of that figure, in increasing order of interest.

The naive one is that the fraction is falling, so the programme is working. It is falling. From 6.2% in iJR904 to 4.5% in iML1515 is a real gain and it was earned by fourteen years of manual curation, targeted experiments and gap-filling loops of the kind §3 described.

The second is that the absolute count is rising, from 58 to 107. Every reconstruction in the series found and closed gaps, and every reconstruction also expanded into parts of metabolism where the gene assignments were weaker. The frontier moves outward faster than the interior fills in.

The third is the one that settles the question. Follow the individual orphans forward. Of the sixteen activities in Table 1 that had no gene in 2000, most now have one: methylthioadenosine nucleosidase is b0159, cysteine desulfhydrase is b3008 or b3708, the pantothenate steps have their coa genes, tetraacyldisaccharide kinase is lpxK. A few dropped out of the model when the chemistry was revised. And two of them, CMP nucleosidase and NMN glycohydrolase, are still in iML1515, still carrying no gene, seventeen years later.

6Catalysis is not regulation

A stoichiometric matrix says what can be converted into what. It does not say what turns anything on.

Suppose the last five sections were wrong and annotation were perfect. You would then have, for every gene, the enzyme it encodes, and for every enzyme, the reaction it catalyses. Assemble those into a matrix S whose rows are metabolites and whose columns are reactions, and you have a genome-scale reconstruction. What can you compute?

Flux balance analysis answers precisely. Impose steady state, S·v = 0, bound each flux, choose an objective, and solve the linear program¹²¹¹. The answer is a flux distribution, and it is a real answer: growth rates on defined media, gene essentiality, byproduct secretion, the yield ceiling of an engineered pathway. Two decades of it worked¹³.

Now list what is not in S. Not one entry of that matrix records that ATP inhibits phosphofructokinase, that a repressor sits on an operator, that a transcription factor is sequestered by an anti-sigma factor, that two enzymes compete for the same ribosome. Every one of those is a binding event. The point generalises: retroactivity, the reason a downstream module changes an upstream one's behaviour, is binding and nothing else²⁸. Binding does not consume or produce anything on the timescale of metabolism, so it has no column in a stoichiometric matrix, and it is the entire mechanism of regulation.

7Binding and catalysis

The model class this document argues in. Two layers, separated by timescale, with different mathematics and different parameters.

Split every reaction in a cell into two kinds. A binding reaction is a reversible association or dissociation: nothing is made or destroyed, atoms are only rearranged into and out of complexes. A catalysis reaction is everything else: it changes what the cell is made of. The split is chemical, not a modelling convenience, and it is old. Monod's chapter III of Chance and Necessity names exactly these two as the operations a protein performs, the enzyme-proteins as specific catalysts and the noncovalent stereospecific complex²³.

Write the state as x, listing first the d atomic species, meaning those not composed of other tracked species, and then the r complexes. The atomic content matrix records what each complex is made of:

Binding conserves atomic content, so the totals

are unchanged by any binding reaction, and are moved only by catalysis. Because binding is fast, the complexes sit at equilibrium given the totals, and mass action makes that equilibrium a monomial:

Given the totals q, equations (2) and (3) determine every concentration. That system has exactly one positive solution, by Birch's theorem, so the map from totals to state is a well-defined function²¹. Catalysis then runs slowly on top of it:

The timescale separation that licenses this is the only approximation in the class, and it is a controlled one. Binding equilibrates in microseconds to seconds, metabolism in seconds to minutes, expression in minutes to hours. Singular perturbation theory says exactly when the algebraic limit is valid and what the error is when it is not²².

8Bailey's four invariances

Bailey's strongest objection to structural reasoning is a list of four perturbations it cannot see. Every item on the list is a binding phenomenon.

Having argued that structure is nearly free, Bailey then argued against himself, and the argument is better than the thesis. He points out that the structural theories he has been recommending are blind to a specific class of change.

"Some of the potential difficulties in applying any of the above theories to biological systems can be illustrated by the following fact: results from all of the theories presented here are completely unchanged by any change in the cell that does not change the set of reactions that occur. More specifically, results from all of these theories are invariant to any and all of the following genetic modifications in the organism: First, insertion or deletion of any DNA fragment or molecule that does not eliminate or knock out a gene for an enzyme or add a gene for an enzyme not already present in the organism's genome; second, up- or downregulation of expression of an enzyme (or any combination of enzymes) or of any other protein; third, any mutation that changes the allosteric regulation of any enzyme; and fourth, altered distribution of isoenzyme expression (including knockouts of all except one). It is well known experimentally that genetic modifications of these types can strongly influence specific growth rates, metabolic flux distributions, and other quantitative and qualitative aspects of cell function." Bailey 2001¹

Read the premise carefully: any change that does not change the set of reactions that occur. The objection is exactly as strong as that premise, and the premise depends entirely on what counts as a reaction. Bailey's theories count catalytic reactions. Count binding reactions too and the list falls into two halves, neither of which survives.

Items three and one become structural changes outright. A mutation that changes the allosteric regulation of an enzyme adds or removes a binding reaction. In the notation of §7 it adds or deletes a column of L₂, which is as structural as a change gets. An inserted DNA fragment that carries a binding site adds a decoy species that titrates a regulator away from its targets, which is another column, and the effect on downstream expression is large and has been measured directly, both for a single site²⁶ and for arrays of them²⁷²⁶. Both changes are invisible only to a model whose reactions are all catalytic. Bailey's premise is false as soon as binding counts.

Items two and four are changes to conserved totals, and the structure predicts their consequences. This half is more interesting, because expression level really is not an edge in any network. It is a total q, which is an input. What the binding structure supplies is the function that maps totals to state, and in particular the set of achievable log-derivatives of that map, computed from L alone with no rate constant. That set is the reaction order polyhedron, and the companion tutorial The biomachine perspective builds it from scratch in its sections 12 to 16²⁵²⁴.

Take the smallest possible instance: one enzyme, one substrate, one complex. The observable is the complex, and the question is how it responds to the two totals.

Every point in that triangle is an exact local power law of the response, and which point you are at is set by the totals. So the structural object is not a single number, it is a geometry, and a change in expression is a walk across it. Raise the enzyme six decades and ∂log C / ∂log q_E falls from 1.000 to 0.001, while ∂log C / ∂log q_S rises to 1.000. That is a qualitative change in behaviour caused by a purely quantitative perturbation, and it was computed from L with no rate constant anywhere.

That much is a calculation. The next three paragraphs are the measurement, because someone did this at genome scale, in an unrelated tradition, and did not read it as an answer to Bailey.

A note on names first, because three literatures converge on this one quantity and it is easy to read them for years without noticing. What this section calls a reaction order, metabolic control analysis calls a scaled elasticity, and biochemical systems theory calls a kinetic order. All three are ∂ log v / ∂ log x. Savageau defined the kinetic order in 1969 as the derivative of the rate with respect to a concentration, times the concentration, divided by the rate, evaluated at an operating point⁴⁰, and Voit's review states the identity with the elasticity outright⁴². The worked example there computes the kinetic orders of a Monod law with product inhibition as K₁/(K₁ + c₁) and −c₂/(K₂ + c₂), which are the substrate and non-competitive-inhibitor elasticities of §6 with the letters changed.

The polyhedron has a relative in that literature as well. Savageau's design space partitions a system's parameter space into regions according to which power-law term dominates each balance, and treats each region as a qualitatively distinct phenotype⁴⁴. The vertices of a reaction order polyhedron are dominance patterns, so the two constructions enumerate the same combinatorial object from opposite ends. The design space asks which parameter values put the cell in a given region. The polyhedron asks which log-derivatives a region permits.

The quantity metabolic control analysis calls the scaled elasticity is the reaction order under another name. For a species x acting on a rate v it is ε = (∂v/∂x)(x/v), which is ∂ log v / ∂ log x. Reznik, Christodoulou, Goldford, Briars, Sauer, Segrè and Noor prove a general result about it²⁹. Write any separable rate law as v = V⁺ κ γ θ(x), where θ is the relative activity of the enzyme and is the only factor the regulator touches. Then

for a cooperative activator or a non-competitive inhibitor with Hill coefficient h, and the result holds for activators and inhibitors alike. Read the two factors separately, because they are the two halves of this document's thesis standing next to each other. The Hill coefficient h is the number of copies of the regulator in the complex, which is an entry of L₂: it is the binding structure, and it bounds the reaction order. The activity θ is a saturation, set by the abundance and the binding constant: it is where the cell is sitting, and it places the reaction order inside that bound. Structure supplies the polytope. Totals supply the point. Reznik and colleagues derive the single-species case of exactly the geometry the figure above computes, and they draw the economic conclusion rather than the structural one: an enzyme cannot be responsive to a regulator without giving up catalytic activity, so regulation has a cost.

They then evaluate ε for every catalogued interaction in E. coli, in thirteen growth conditions, from measured metabolite concentrations and measured binding constants. That is the experiment Bailey's second invariance forbids from mattering. The carbon source changes. The reaction set does not.

The median interaction moves by 0.07, which is small. The distribution has a long tail, which is not: 59 of the 186 move by more than 0.1, 38 by more than 0.2, and the largest by 0.74. And the tail is not noise, because the two largest movers are the two interactions the same paper singles out on physiological grounds. PEP inhibiting phosphofructokinase is near-inert on glycolytic carbon, where it is not needed, and carries a reaction order of 0.70 on succinate. FDP inhibiting glycerol kinase does the opposite. A factor of seven in one, twelve in the other, in the same cell, with the same genome, the same enzymes and the same reactions.

The proteome shows the same thing from the other side, and more recently. Seeger, Pinheiro and Lässig fit a Michaelis-Menten network model to E. coli under graded glucose limitation and ask what sets each enzyme's response³⁰. The answer is saturation. Strongly saturated enzymes keep their efficiency as metabolites are depleted and are downregulated in proportion to flux. Weakly saturated enzymes lose efficiency and are upregulated to compensate. So two enzymes in the same pathway, carrying the same flux, move in opposite directions, and which way is set by their binding constants. Their growth law is q_i ∼ ξ_s (σ_i)^−y with y > 0, tested against measured saturations for 59 enzymes in nine pathways at p < 10⁻⁴, and against expression alone in 18 of 34 pathways at p < 10⁻¹⁵.

Their negative control is the part worth pausing on. Protein complexes with fixed subunit stoichiometry, and genes sharing an operon, are the cases where the structure forbids independent tuning. In those the correlation vanishes: the ribosome, NADH:quinone oxidoreductase, tryptophan and histidine biosynthesis all show y ≈ 0. So the heterogeneity is not scatter. It appears exactly where binding parameters are free to differ and disappears exactly where they are yoked, which is the signature of a structural cause rather than a noisy one.

The same machinery answers item four. The extensions of FBA that reach past pure stoichiometry make the same move from the other side, adding enzyme abundance as an explicit constrained resource¹⁴⁸¹. Two isoenzymes catalysing the same reaction are two species with two rows of L and two sets of complexes. Redistributing expression between them is a change in two totals, of exactly the kind the polyhedron describes. Only a model that has already lumped them into one column of S is blind to it, and the lumping, not the theory, is what caused the blindness.

One consequence of that is worth stating now, because the rest of the document runs on it. What counts as structure is not a property of the cell. It is a property of the model's scope. Allosteric regulation is a parameter to a model of catalysis and a structure to a model of binding. Widen the scope and quantities move across the line. That is the first hint that the line is not where Bailey put it, and §10 is the argument that there is no line.

9There is no gene for a binding event

The catalysis layer at least has a draft. The binding layer has never had one, and the reason is not effort.

A gene encodes a protein. A protein is a node. A binding event is an edge, and edges are not encoded anywhere as such. The information that two proteins associate, or that a metabolite sits in an allosteric pocket, is distributed across both partners' folds and is not locally readable from either sequence.

So there is no analogue of the gene-to-reaction rule for binding, and consequently no analogue of the draft reconstruction. Where the catalysis layer starts from an annotation and gets repaired, the binding layer starts from nothing and gets assembled interaction by interaction.

How incomplete is it? The question is answerable, because someone measured a whole subnetwork exhaustively rather than waiting for the literature to accumulate. Diether, Nikolaev, Allain and Sauer took E. coli central carbon metabolism, which is the single best-studied metabolic subnetwork in biology, and screened it by ligand-detected NMR³¹.

There is a genome-scale count to set beside that, and the two together give the shape of the problem. Reznik and colleagues assembled what they call the small molecule regulatory network of E. coli by mining the whole biochemical literature through BRENDA and BioCyc and mapping the result onto iJO1366²⁹. This is the binding layer's nearest thing to a draft reconstruction, and it is worth being precise about what it contains and what it rests on.

Three readings of that table, in increasing order of how much they should worry a modeller. First, the coverage is not nothing: half the enzymes in the model have at least one catalogued regulator, which is a real body of knowledge and considerably more than any kinetic model uses. Second, the support is asserted far more often than it is measured: 89% of the entries have no binding constant at all, so the network is a list of which b_j are nonzero with almost no information about how large they are. Third, and worst, three quarters of the entries rest on a single paper. Reznik and colleagues are explicit that the two-report subset is the defensible one, and it is 325 edges.

The screen's raw matrix is worth reading for its shape as well as its size, because a support that is concentrated is a far easier thing to recover than one spread evenly. Recomputing from Diether's Dataset EV2: the 98 interactions fall on 24 of the 29 enzymes and only 28 of the 55 metabolites, so half the panel bound nothing at all. A binding enzyme holds a median of four metabolites and at most eleven, and the five busiest enzymes carry 43% of the interactions. The busiest metabolite is GTP, at thirteen.

That is the empirical answer to a worry the sparsity view invites. If every protein bound every metabolite a little, the support would be dense and useless as a structure. It is not. It is sparse, uneven and heavy-tailed, which is the regime in which support recovery works at all.

Now put the two measurements together, because they disagree in an informative way. The catalogue lists 1669 interactions over 364 enzymes, which is 4.6 per regulated enzyme. Diether's exhaustive screen finds 98 over the 24 enzymes that bound anything, which is 4.1 per enzyme, against a panel of only 55 metabolites. The densities agree to within twenty percent. The memberships do not: 76 of Diether's 98 were absent from the literature the catalogue was mined from. So the accumulated record has roughly the right number of edges per enzyme and largely the wrong edges. That is precisely a support-recovery failure, and §18 is about why it is the expected one.

That is the state of the art for a subnetwork of 29 enzymes. Genome-scale, the ratio is worse, and the methods that scale are the ones that trade confidence for coverage: affinity-purification and proteome-scale interaction mapping for protein-protein binding⁶⁶⁶⁷, limited proteolysis mass spectrometry for protein-metabolite binding³², and dynamic metabolite data fitted against a kinetic model to select among putative interactions, which is how Link, Kochanowski and Sauer tested 126 candidate allosteric interactions in glycolysis and identified the ones that actually govern the switch between glycolysis and gluconeogenesis³³.

Note what that last method is. It proposes a library of candidate binding reactions, fits, and selects the subset the data support. That is structure discovery by sparsity, done by hand, in 2013, for one pathway. Hold onto it. It is the same procedure as gap-filling in §3, in the other layer. It had also been automated seventeen years earlier and then left alone, which is what §24 is about.

10The statement

In a binding-and-catalysis network there is exactly one formal difference between having a reaction and not having it, and it is not a difference of kind.

Return to equation (3). Complex j sits at c_j = b_j ∏_i x_i^(L₂)_ij. Now ask what it means for that complex not to be in the network. It means c_j = 0 at every state, and since the monomial is positive whenever the concentrations are, that happens exactly when b_j = 0.

The statement is not deep. That is the point of it. Once the candidate columns are written down, there is no separate discrete object called the network, and no separate act of choosing one. There is a vector of nonnegative numbers, and the network is the list of places where it is not zero.

Everything Bailey called structural is therefore a statement about which coordinates vanish, and everything he called parametric is a statement about the values at the coordinates that do not. Qualitative and quantitative are not two kinds of information. They are the zeroth and the higher-order digits of one kind.

That much has been said before, and the honest thing is to say by whom and in what form. Biochemical systems theory has treated a network as a support since the 1960s. Savageau's power-law formalism writes every process as a product of powers of the concentrations, and the exponent on each concentration, its kinetic order, is the log-derivative of the rate with respect to that concentration⁴⁰. The model-building rule printed in every BST text is the one this section just derived: a positive effect gives a positive kinetic order, an inhibitory effect a negative one, and no effect gives exactly zero⁴¹.

So the network of a BST model is the pattern of nonzeros in its kinetic-order matrix, and Voit's review says as much in one sentence: for canonical models the transition between parameter estimation and structure identification is fluid, because the structure of the inferred system changes when a kinetic order becomes zero⁴². Forty years of practice stands behind that clause.

The method was built as well as stated. Hatzimanikatis, Floudas and Bailey attached a binary variable to every candidate regulatory elasticity and solved for the support and the values in one mixed-integer linear program⁴³. Their name for the candidate set is the one §11 is about to need: a regulatory superstructure, in which every metabolite may potentially regulate every enzyme. Run on the aromatic amino acid pathway of E. coli, which carries eight feedback inhibition loops, the program found that inactivating three of them and overexpressing three enzymes raises phenylalanine selectivity by 42%, and that a structure with two loops, neither of them present in the wild type, does better than all eight.

What survives as new, then, is not the continuum. It is what the support is a support of. A kinetic order and an elasticity are log-derivatives evaluated at an operating point. A formation constant is a constant.

The difference is not cosmetic, and §14 is where it bites. An interaction whose site is saturated has an elasticity of zero and a formation constant that has not changed, so it leaves the support of the kinetic-order matrix while remaining in the support of b. Reznik's trade-off | ε | = h (1 − θ) from §8 is the exact conversion between the two²⁹: the Hill coefficient is the structure and it does not move, the activity is the saturation and it does. So the object BST calls structure is the resolvable structure of §14 rather than the chemistry, which is why its practitioners had to prune small kinetic orders to zero by hand⁴². The statement boxed above is the condition-independent version, and it is available only because binding and catalysis were separated first.

11What the claim costs

The statement is conditional on an enumeration, and the enumeration is where the difficulty went. It did not disappear, and pretending otherwise would be the interesting way to be wrong.

The phrase doing the work in §10 is with all candidate columns enumerated. Setting b_j = 0 removes a complex that the model already has a column for. It does not remove a complex the model never wrote down, and it does not create one.

So the reduction of structure to sparsity is a reduction relative to a superstructure: a fixed, over-complete list of candidate complexes and candidate catalytic steps, of which the real network is a subset. Three things follow, and the first two are costs.

First, the superstructure has to contain the truth. If a real complex has no column, no value of b produces it. This is not a technicality that better computation removes. It is the same assumption gap-filling makes when it searches a universal reaction database, and it fails in the same way: the database is a catalogue of reactions someone has seen, so a genuinely novel chemistry is outside it by construction.

Second, the superstructure is large. The number of candidate columns grows combinatorially in the number of atomic species. For iML1515, which accounts for 1,192 unique metabolites¹⁰ and roughly 1,500 gene products, the pairwise enzyme-metabolite candidates alone number in the hundreds of thousands. Search over subsets of a set that size is not a search anyone will perform.

Third, and this is the part that is a gain rather than a cost: the superstructure is where prior knowledge enters, and it enters in the right place. Everything the genome does supply, everything in a curated database, everything a structure predictor proposes, is a statement about which columns are plausible. None of it has to be a statement about which are real. Bailey's structural information is not discarded by this view. It is demoted from a specification to a prior, which is what the evidence in §5 says it always was.

12Deletion is a limit, not an edit

There is a second, independent route to the same conclusion, and it comes from the geometry of model fitting rather than from chemistry.

Work in log parameters, which is the natural coordinate for rate and equilibrium constants because they span orders of magnitude and enter the equations multiplicatively. Then b_j = 0 is not a point in parameter space at all. It is log b_j → −∞, a boundary at infinity.

This is exactly the object the sloppy-models literature found by a different route. Fit a multiparameter model to data and the achievable predictions form a model manifold in the space of possible data. Gutenkunst and colleagues observed that this manifold is extremely anisotropic: parameter combinations span many decades of sensitivity, and most of them are unconstrained by any feasible measurement³⁴. Transtrum, Machta and Sethna developed the geometry³⁵, and Transtrum and Qiu turned it into an algorithm.

"We propose a new approach by translating the model reduction problem for an arbitrary statistical model into a geometric problem of constructing a low-dimensional, submanifold approximation to a high-dimensional manifold. When models are overly complex, we use the observation that the model manifold is bounded with a hierarchy of widths and propose an approximation by successively removing the thin widths, i.e. by finding the manifold boundaries." Transtrum & Qiu, Model reduction by manifold boundaries, Phys Rev Lett 113:098701, 2014³⁶

Their method walks to a boundary of the model manifold and reads off the limit that reaches it. Those limits are always of the same kind: a rate constant to zero or infinity, two species merging, a timescale separating. And every one of them is a structural simplification. A reaction disappears, two nodes become one, a fast step becomes algebraic.

One consequence is worth noticing because it explains an old irritation. Modellers frequently find that a simpler model fits as well as a complex one, and treat this as a nuisance of model selection. In this picture it is the generic situation: the data localise you to a neighbourhood of the manifold, that neighbourhood touches boundaries, and every boundary it touches is a simpler structure that fits the data. §16 measures how close those boundaries are.

13Structure is a threshold

A support is a discontinuous functional of a continuous vector. No finite measurement determines one. What a measurement determines is which coordinates exceed a threshold, and the threshold is a property of the instrument.

Here is the awkward corollary of §10. If the network is the set of j with b_j > 0, then to know the network you have to distinguish b_j = 0 from b_j = 10⁻⁹. Data never does that. What data does is bound b_j above, and everything below the bound is indistinguishable from absence.

That bound is computable, so let us compute it. Take the smallest interesting case: a scaffold with two binding sites for a ligand, strongly cooperative, so that the doubly bound complex is well populated and the singly bound intermediate may or may not be. The question is how tightly the intermediate must bind before a binding assay can tell that it exists.

The numbers are worse than intuition suggests. On this system, at 5% measurement error, an intermediate binding more weakly than K_d = 1.1 µM leaves no trace a 25-point titration can find. Push the precision to 0.2%, which is beyond most binding assays, and the floor moves only to 12.4 µM. A factor of twenty-five in measurement precision buys a factor of eleven in the weakest affinity you can see, and the returns are worse than that outside the linear regime.

And this is the easy version. The intermediate here competes for its partners with a much stronger complex, which absorbs it. The same calculation for a candidate that competes with nothing gives a floor roughly an order of magnitude weaker: for a species at 1 µM total, measured to 5% on three replicates, K_d = 69 µM. Interactions weaker than that are, at that abundance and that precision, not in any network the data supports.

14The threshold moves with the cell

The floor is proportional to abundance. Abundance is a property of the condition. So the resolvable structure is a property of the condition, and there is no such thing as the network of a cell independent of the state it was measured in.

Read the closed form again: K_d,floor = 2 q √n / s. Everything on the right except q belongs to the experiment. And q, the total abundance of the species, is set by the cell.

Raise the abundance tenfold and the floor moves out tenfold: interactions that were invisible become visible, because a weak constant times a bigger total is a bigger complex. Lower it and they vanish again. The same chemistry, measured in the same organism with the same instrument, yields a different structure in a different growth condition.

The consequence has been measured, and it is the same measurement §8 used for a different purpose. Reznik and colleagues' elasticities are computed from metabolite concentrations, so they move when the concentrations move: across thirteen carbon sources, 59 of 186 catalogued interactions shift their reaction order by more than 0.1, and the largest by 0.74²⁹. An interaction whose reaction order is 0.06 in one condition is, for practical purposes, not in the network of that condition. The same interaction at 0.75 in another condition plainly is. Nothing about the chemistry changed.

This is not a small effect in bacteria. Protein abundances in E. coli move over orders of magnitude across nutrient and stress conditions, and they do so systematically rather than randomly, because the proteome is allocated⁸². So the set of complexes above the detection floor is reorganised by exactly the perturbations physiologists spend their time on.

The practical consequence is a rule about what to fit. A model whose parameters are dissociation constants and turnover numbers has a chance of being condition-independent, because those quantities are. A model whose parameters are effective interaction strengths fitted per condition does not, because it has absorbed the abundances into the coefficients. The two look similar on the page and behave differently the moment you extrapolate, which is the only reason to have built a mechanistic model at all.

15Two networks, one dynamics

The floor of §13 is a limit of precision. This is worse: even with perfect, noiseless, complete observation, the network is not determined.

Suppose you could measure every species in a cell, continuously, exactly. Could you then write down the reaction network? The answer has been no since 1981, when Hárs and Tóth posed the inverse problem of reaction kinetics³⁷, and it was made precise by Craciun and Pantea.

"We show that there exist reaction networks R for which the reaction rate constants are not uniquely identifiable, even if we are given complete information on the dynamics of concentrations for all chemical species of R. Also, we show that there exist reaction networks R₁ ≠ R₂ such that their dynamics are identical under appropriate choices of reaction rate constants." Craciun & Pantea, Identifiability of chemical reaction networks, J Math Chem 44:244, 2008³⁸

The mechanism is elementary once seen. A mass-action vector field is a sum of terms, one per reaction, each a reaction vector times a rate times a monomial:

Two networks agree as vector fields whenever, monomial by monomial, the rate-weighted reaction vectors sum to the same thing. And a vector can be split. Take any reaction and replace it by two reactions from the same source complex whose reaction vectors add up to the original.

The two networks are different biology. Network 1 says A decomposes into B and C. Network 2 says A converts to B and, separately, A catalyses the production of C without being consumed. One is a lysis, the other is an enzyme. No trajectory distinguishes them.

What breaks the tie? In practice, sparsity does. Szederkényi turned the choice into an optimisation and made the criterion explicit.

"A numerical procedure for finding the sparsest and densest realization of a given reaction network is proposed … The problem is formulated and solved in the framework of mixed integer linear programming, where the continuous optimization variables are the nonnegative reaction rate coefficients, and the corresponding integer variables ensure the finding of the realization with the minimal or maximal number of reactions." Szederkényi, Computing sparse and dense realizations of reaction kinetic systems, J Math Chem 47:551, 2010³⁹

Read the sentence structure. The continuous variables are the rate constants. The integer variables count how many are nonzero. That is the thesis of §10, written down in 2010, as an algorithm, for the catalysis layer. The structure is not an input to the problem. It is the cardinality of the support of the solution, and the modeller chooses it by choosing an objective.

16What a titration does fix

Non-identifiability is not total. Some structural facts survive the data and some do not, and which is which is measurable rather than a matter of opinion.

The two negative results so far could be read as saying structure is unknowable, which would be both wrong and useless. So here is the positive version, on the binding layer, measured.

Generate a titration from a known two-complex truth. Fit four candidate structures to it: the truth, a one-complex model keeping only the doubly bound species, a one-complex model keeping only the singly bound species, and an over-specified three-complex model. Compare each fit's residual to the measurement noise.

Three things come out of that, and they are different in kind.

Some structure is fixed. Keeping only the singly bound complex gives a residual of 1.7 × 10⁻¹, more than six times the noise at 5% error. The data refutes it decisively. Cooperativity, in the sense of “the dominant complex carries two ligands, not one”, is recoverable from a dose-response curve, which is why Hill coefficients have been useful for a century.

Some is not. Keeping only the doubly bound complex gives 3.6 × 10⁻³. That is below the noise at 1% relative error and far below it at 5%. The number of complexes, which is the most basic structural fact there is, is not determined by a 25-point titration at ordinary precision.

And over-specification is harmless when the estimator prefers sparsity. The three-complex model fits exactly and sets the third constant to 1.4 × 10⁻¹⁷. Given a superstructure that contains the truth and clean data, the extra columns switch themselves off. This is the empirical case for the strategy of §11, and it is also, honestly, the easy case: the data here is noiseless and the candidates are nested. §18 is about what happens when neither holds.

17Where structure alone still decides

Bailey's second claim, that structure gets you a long way without any parameter values, is not merely defensible. In its strongest form it is a theorem.

It would be easy to leave this document having argued that structure is soft, uncertain and instrument-dependent, and to conclude that structural reasoning is therefore weak. That conclusion is wrong, and the counterexample is the sharpest result in the field Bailey was pointing at.

Shinar and Feinberg ask when a network guarantees that some species sits at the same concentration in every positive steady state, no matter how the totals are supplied. They call it absolute concentration robustness, and the conditions are combinatorial.

The authors state the force of it plainly: these are structural attributes that will impart ACR to any mass-action network having them, regardless of the values that the rate constants take. Counting nodes and taking a difference gives a conclusion about behaviour that no amount of parameter measurement could overturn. This is Bailey's argument, vindicated in its strongest available form, and it descends directly from the Horn and Jackson theory he cited²¹, whose modern algebraic form is the theory of toric dynamical systems⁴⁵²¹.

Two things are worth saying about how it fits with the rest of this document, because they cut in opposite directions.

It cuts for structure, because it says the support is not merely a scaffold for parameters. It carries theorems. Every bit of support you actually know buys conclusions that hold for a whole open set of parameter values, which is why §11 insists that prior structural knowledge is valuable even after being demoted to a prior.

And it cuts against the easy reading of Bailey, because the networks the theorem applies to include their binding steps explicitly. Deficiency is computed from the complete node set, and the two non-terminal nodes differing in one species are, in the canonical examples, a free enzyme and an enzyme-substrate complex. The theorem is a theorem about a network that has a binding layer in it. Compute the deficiency of a stoichiometric reconstruction, which has lumped every complex away, and you are computing the index of a different network.

18The statistics of support recovery

Sparse estimation has sharp results about when the support can be recovered. Read them as statements about biology and they are not encouraging.

If structure is a support, then the relevant theory is the theory of support recovery, and that theory is mature. This is also, quietly, what sparse identification of nonlinear dynamics does when it is applied to biology: propose a library of candidate terms, fit with a sparsity penalty, keep the survivors⁵⁰⁵¹. Two results matter here.

Sample complexity is a threshold, not a trend. Wainwright analysed ℓ₁-constrained quadratic programming, which is the lasso⁴⁷, and showed that support recovery undergoes a sharp transition: below a sample size scaling like 2k log(p − k) for a support of size k among p candidates, the probability of recovering the correct support goes to zero, and above it, to one⁴⁸. There is no partial credit regime in between. Put the numbers of §11 into that expression and the required number of independent, informative conditions is large.

And recovery requires a condition on the candidates that biology is designed to violate. Zhao and Yu identified what the lasso needs in order to select the right variables: the irrepresentable condition, which says roughly that the irrelevant candidates must not be too correlated with the relevant ones⁴⁹. Violate it and the estimator is consistent for prediction and inconsistent for selection: it fits beautifully and picks the wrong variables.

The parameter-side version of the same difficulty is well mapped. Profile likelihood separates the parameters a dataset constrains from those it does not⁵², and structural identifiability analysis decides in advance which combinations could ever be determined⁵³. The topology-side version is the one Babtie, Kirk and Stumpf named. They showed that in systems-biology model fitting, uncertainty in the network topology typically dominates uncertainty in the parameters, and that many topologies fit the data comparably well⁵⁴. That is the empirical shadow of the two theorems above, and their recommendation, to work with sets of topologies rather than one, is the recommendation of §26.

19What network inference achieves

The theory says support recovery should be hard. The benchmarks say it is, and they say by how much.

Two community efforts have measured this properly, with held-out gold standards rather than self-reported performance.

DREAM5 put the field's methods on E. coli, S. cerevisiae, S. aureus and a simulated benchmark⁵⁷. The headline conclusion was that no single inference method performs optimally across all datasets, and that integrating predictions across methods is more robust than any of them. Two details matter more than the headline. Performance on S. cerevisiae was poor for every method, which the authors attribute partly to the lower coverage of the yeast gold standard and partly to the weak correlation between mRNA and the regulation being inferred. And the community network for E. coli, the best case in the study, was tested: 53 novel predicted interactions were assayed experimentally and 23 were supported, a hit rate of 43%.

BEELINE did the same for single-cell methods, scoring against synthetic networks, curated Boolean models and experimental data⁵⁸. The scale is the AUPRC ratio, the area under the precision-recall curve divided by a random predictor's. On synthetic networks twelve algorithms exceeded a ratio of 2.0 and seven exceeded 5.0. On the curated biological models the best median ratio was 1.4, several methods had median early-precision at or below a random predictor, and some scored worse than random on particular models.

There is a longer record than the benchmarks, and it points the same way. Structure identification on kinetic-order matrices has been an active subfield since the 1980s, using genetic algorithms, particle swarms, simulated annealing, Kalman filters, collocation, branch-and-reduce, mixed-integer programming, and an ℓ₁ penalty with a pruning step⁴². The methods are the ones §18 describes, applied to biology two decades before the benchmarks existed. Voit's summary of the outcome, written from inside the programme, is that the algorithms worked fairly well for some applications and failed for others, and that the estimation community is still awaiting a truly exceptional solution.

That verdict is worth more than a benchmark score, because it is not a statement about any one method. Forty years of method development against a fixed obstacle is evidence that the obstacle is not methodological. §18 says what it is.

20Sequence to catalysis: four links, three of them biochemistry

“The genome gives you the network” compresses a chain of four inferences. Only the first is about sequence.

It is worth laying out the chain explicitly, because the compression is what makes the claim sound stronger than it is.

Link 1 is sequence and it is finished. Link 2 is an inference from homology whose failure modes were §4. Link 3 is not an inference at all, it is a lookup in a table that biochemists filled in. And link 4 is an optimisation against phenotype data, which is to say it is fitting.

So the honest form of Bailey's supply claim is: sequencing solves one of the four links, and the industrialisation of link 3 into databases is what made the other three feel automatic. That is a real achievement and it is not the same as the network being given.

21Sequence to binding structure

This is where the answer has genuinely changed since 2001, and where the published accuracy numbers describe a different task from the one a virtual cell needs.

Bailey could take for granted that nothing predicts a complex from sequence. That is no longer true. AlphaFold solved single-chain structure prediction⁶⁰, AlphaFold-Multimer extended it to assemblies⁶¹, AlphaFold 3 extended it to protein-ligand, protein-nucleic-acid and modified residues⁶², and the question became quantitative rather than categorical.

Three results bound the current answer, and they should be read in this order.

Given a pair that interacts, a model is often right. Bryant, Pozzati and Elofsson found that AlphaFold2 with optimised alignments produces models of acceptable quality for 63% of heterodimers, and that a score derived from the predicted interface distinguishes interacting from non-interacting pairs well enough to identify 51% of interacting pairs at a 1% false positive rate⁶³.

A proteome-scale screen is a different problem, and the difference is arithmetic. A 1% false positive rate is an excellent number on a balanced benchmark and a ruinous one on a candidate set of every pair.

So the real screens run at much stricter thresholds, and pay for it in recall. Humphreys and colleagues screened 8.3 million yeast protein pairs and called 1,505 as likely interacting, which is 0.018% of the candidates⁶⁴. That threshold is 55 times stricter than a 1% false positive rate, and it produced 106 previously unidentified assemblies, which is a substantial result. Burke and colleagues took the other approach and modelled 65,484 already-reported human interactions, of which 3,137, or 4.8%, ranked as highly confident⁶⁵. They also give the baseline: fewer than 5% of the hundreds of thousands of known human protein interactions have ever been structurally characterised.

22Sequence to parameters

Three parameters, three different answers, and the pattern in them is not the one Bailey's split predicts.

Equilibrium constants are obtainable, and almost exactly. The equilibrium constant of a reaction is fixed by the Gibbs energies of its metabolites, which are properties of small molecules rather than of the enzyme. Free energies are additive over chemical groups, so they can be estimated by decomposition and reconciled against measurements. Noor and colleagues' component contribution method does this consistently across a genome-scale network⁷⁰, and it is the substance behind eQuilibrator⁷¹. This is a parameter, in Bailey's sense, that is obtainable from structure. It is the cleanest counterexample to the split in the whole document, and it is a counterexample in the direction nobody expects.

Turnover numbers are not obtainable from sequence, and the strongest published claim to the contrary does not survive its own evaluation. DLKcat predicts k_cat from an enzyme sequence and a substrate structure⁷². Kroll and Lercher re-analysed it against the obvious baseline⁷³.

Read carefully, that is not a criticism of one model. It is a measurement of how much information about a turnover number a sequence carries, and the answer is: about as much as its nearest neighbour's identity, and no more. Turnover is a property of a transition state, and a transition state is not a homology-transferable feature.

Michaelis constants sit in between. They are part affinity and part turnover, and a model over sequence and structural features reaches genome scale with usable but not high accuracy⁷⁵. Structural features do carry some signal about turnover⁷⁷, which is consistent with the nearest-neighbour reading above rather than an alternative to it. And the distribution of the parameters is itself well characterised: enzymes are moderately efficient and cluster far from any diffusion limit⁷⁸, which is exactly why a prior over k_cat helps and a prediction of it does not.

Turnover numbers are, however, measurable at scale in vivo. Davidi and colleagues combined proteomics with computed fluxes to obtain the maximal observed catalytic rate of each enzyme inside cells, and found it agrees with in vitro k_cat at r² = 0.62 on a log scale, root-mean-square difference 0.54, which is 3.5-fold⁷⁶. That is the useful number, and it comes with a coverage limit: the method yielded estimates for 436 E. coli enzymes⁷⁴, against the 1,515 gene products in iML1515. Curated repositories are the other route and their coverage is the binding constraint: BRENDA⁷⁹ and SABIO-RK⁸⁰ hold on the order of a hundred thousand kinetic parameters across all organisms, of which only a few thousand are E. coli.

23The obtainability table

Every object in the model class, the data that constrains it, and what that data is worth. This is the direct answer to the question the document opened with.

Three patterns in that table are worth naming, and none of them is the pattern Bailey's split predicts.

The split does not run between structure and parameters. It runs between quantities determined by small molecules and quantities determined by proteins. Equilibrium constants are obtainable because they are set by metabolite thermodynamics. Turnover numbers and binding affinities are not, because they are set by a fold. That the obtainable one is a parameter and the unobtainable ones include a structure is the whole argument of this document in one row of a table.

Coverage and fidelity fail separately, and confusing them is how the field talks past itself. A k_cat predictor has unlimited coverage and no fidelity. An in vivo measurement has good fidelity and 29% coverage. A structure predictor has good fidelity on a curated set and unusable precision on the real candidate set. Any claim that something is “obtainable” has to say which of the two it means.

And everything unobtainable is unobtainable in the same way. It has to be inferred from a fit, against a candidate set, from data that does not determine it uniquely. Which is to say: the objects Bailey put on opposite sides of his split turn out to need the same procedure, and Part VI is about what that procedure has to look like.

24Gap-filling is sparse structure learning

The procedure this document has been arguing for has been in production since 2006. At genome scale it runs on the catalysis layer and nowhere else, but the regulatory version was written first, on small pathways, and then left.

Put the pieces of Part I next to the pieces of Part III and they are the same object.

Szederkényi's realisation algorithm makes the identity exact rather than analogical: continuous variables are the rate constants, integer variables count the nonzeros, and the objective selects the sparsest network consistent with the dynamics³⁹. That is §10, as a solver, for mass-action kinetics, published in 2010.

It was not the first such solver, and the earlier one ran on the layer this document says has no algorithm. Fourteen years before, Hatzimanikatis, Floudas and Bailey put a binary variable on every candidate regulatory elasticity of a pathway and solved for the support and the strengths together⁴³. The regulatory superstructure they postulated is the complete bipartite candidate set, every metabolite against every enzyme, and they said plainly that enumerating it leads to a large combinatorial problem. Their instance was one pathway with eight loops and seventeen binary variables, solved in a quarter of a second.

So the missing ingredient was never the algorithm. It was the superstructure, which for a pathway of six reactions can be written down by hand and for a genome cannot. That is what changed in 2017: Reznik's catalogue of 1669 interactions over half the EC numbers in §9 is a measured draft of exactly the object the 1996 formulation needed and had to invent. Nobody has yet run the 1996 program on the 2017 superstructure, and the reason is not that it would be hard to set up.

The field also already knows the correct response to non-identifiability, and has implemented it. Biggs and Papin observed that draft reconstructions admit many alternative network structures, all equally consistent with available data, and built EnsembleFBA to make predictions from the whole set rather than from one member⁸⁷. Bernstein and colleagues surveyed the same problem across the reconstruction pipeline⁸⁸. An ensemble of gap-filled networks is a sample from a posterior over supports, arrived at empirically because the alternative did not work.

25Two stages inherit an error one stage would not

The standard construction fixes the structure, then fits the parameters. Under §10 that is not the natural order of work. It is a constraint, and it has a cost.

Every mechanistic whole-cell model has been built the same way: assemble the network from the literature, then find values for its parameters, also from the literature. Karr and colleagues' M. genitalium model rests on a synthesis of over 900 publications and more than 1,900 experimentally observed parameters across 28 submodels⁸³. The E. coli whole-cell model was constructed by populating the network with literature-derived parameters and then using simulation to find where those parameters were mutually inconsistent⁸⁴. The minimal-cell simulation follows the same pattern⁸⁵.

These are serious achievements and the method has a specific structural weakness. Fixing the structure first makes an assertion that the data has not been asked about, and the second stage cannot revise it. If a complex is missing from L₂, no value of b repairs the omission. The fit will still converge, because a flexible model with enough parameters absorbs the discrepancy into whichever coefficient is nearest. What comes out is a model that fits and has the wrong mechanism, which is the failure mode that matters, because the only reason to build a mechanistic model is to extrapolate.

The reason this is feasible for binding-and-catalysis networks specifically, rather than for reaction networks in general, is that the forward map is well behaved. Equations (2) and (3) have exactly one positive solution for every choice of totals, so the model is a genuine function of its parameters rather than a branch of a multivalued one, and it is differentiable, so the penalised loss has a gradient. Those two facts are what make (6) an optimisation rather than a search.

The honest counterpoint is that (6) is not obviously better in practice, and this document does not claim it has been demonstrated at scale. What it claims is narrower: the two-stage construction is a choice, its cost is a first-stage error the second stage cannot see, and §5 establishes that first-stage errors are the normal case rather than the exception.

26The recipe

Six steps, applicable to a pathway or to a cell. The only unusual one is the first.

One. Write down a superstructure, not a network. Enumerate candidate complexes and candidate catalytic steps generously, including the ones you doubt. This is the step that decides what is findable, and it is the only irreversible commitment in the procedure.

Two. Put every piece of prior knowledge into the penalty, not into the topology. A curated interaction gets a small λ. A structure predictor's proposal at 51% recall gets a larger one⁶³. A homology-transferred annotation gets a larger one still. Nothing gets a hard zero, and nothing gets asserted.

Three. Fit both layers at once, using equation (6), against data from as many conditions as you have.

Four. Report the floor alongside the network. For every candidate the fit rejected, compute the smallest constant that would have been detectable at the abundances in the data, exactly as §13 does. “Not in the model” and “below K_d = 50 µM at the abundances we measured” are different claims, and only the second is true.

Five. Keep an ensemble. Refit from many starts, or sample the posterior, and carry forward the set of supports that fit within noise rather than the best one. This is what EnsembleFBA does for the catalysis layer⁸⁷ and what topological sensitivity analysis recommends generally⁵⁴. A prediction on which the ensemble disagrees is the next experiment.

Nothing in those five steps is specific to a small system, and none of them is what the current virtual-cell programmes propose. The published priorities for building a virtual cell with machine learning are about representation, data volume and evaluation⁸⁶. All three matter. None of them addresses what the fitted object is a property of, which is the question §14 says decides whether the model extrapolates.

Six. Perturb, do not merely observe. §18 gives the reason: support recovery fails when candidate columns are correlated, and in a cell they are correlated by chemistry and by descent. Perturbation is what decorrelates them, which is why a hundred perturbation conditions are worth more than a thousand observational ones, and why the argument for perturbation data here coincides with the one causal inference makes⁵⁵.

27Where this is wrong

Four places, in decreasing order of how much they would cost if they turn out to matter.

The superstructure has to contain the truth, and there is no way to check. This is the load-bearing assumption and §11 stated it plainly rather than burying it. If a real complex is not enumerated, the procedure will fit around its absence and report a support with high confidence. The only defence is the one gap-filling already uses, which is to keep the candidate set embarrassingly large, and it is a partial defence: a genuinely novel chemistry is outside any catalogue by construction. This is the failure mode that would make the whole approach quietly wrong rather than visibly wrong, which is the worse kind.

Exact zeros are a modelling convention, not a fact about cells. Any two molecules in solution associate to some degree. So the true b has no zeros at all, and the sparse model is an approximation whose error is the sum of everything below the floor. Usually that sum is negligible. It is not negligible when many weak interactions share a partner, because then the aggregate sequestration is large while every individual term is invisible. Whether that case is rare in cells is an empirical question this document does not answer, and Diether's finding that some enzymes bind eleven different metabolites³¹ is not reassuring.

The timescale separation is an approximation, and it is the only one in the model class. Equations (2) and (3) are exact in the limit of infinitely fast binding. At finite separation the error is controlled by an explicit small parameter with a persistence theorem behind it²², which is a better epistemic position than a fitted smoothing function occupies, and it is still an approximation. Where binding and catalysis run at comparable rates, the algebraic layer has to be integrated instead, and the cheap gradient goes with it.

And the analogy to the lasso breaks in at least one place that matters. Sparse recovery theory is developed for linear models. The map from b to observables here is a monomial system, positive, and strongly nonlinear at the crossovers that §8 showed are where the interesting behaviour lives. The sample-complexity thresholds⁴⁸ and the irrepresentable condition⁴⁹ are quoted here as the right shape of result rather than as theorems that apply as stated. Establishing the corresponding conditions for equilibrium binding networks is open, and it is the piece of theory this argument most needs.

28What the argument comes to

Three claims, in the order the evidence supports them.

The first is a matter of record. Genome-scale networks are not read off genomes. The first E. coli reconstruction needed 46 entries from the biochemical literature, sixteen of which named an enzyme activity with no gene at all, and after seventeen years of the best-resourced curation effort in biology two of those sixteen still have no gene. The gene-less share of the network fell from 6.2% to 4.5% over fourteen years, which is progress at a rate that reaches zero in the 2040s. What closes a gap is an experiment, not a better reading of the sequence.

The second is conceptual, and it is why the first is not a scandal. There was never a separate discrete object to be supplied. A complex is in a network exactly when its formation constant is nonzero, so the structure is the support of the parameter vector, and the qualitative and the quantitative are the zeros and the nonzeros of one list. The same conclusion arrives independently from the geometry of model fitting, where reduction means walking to a boundary of parameter space and every boundary is a simpler structure.

The third follows from the second and is the one with consequences. If structure is a support then it inherits everything true of supports: it is not identifiable from dynamics, it is defined by a threshold that belongs to the instrument, and the threshold moves when the cell changes its abundances. So “the network of a cell” names two different objects that are routinely conflated, and only one of them is observable.

Bailey's split was the right strategic call and the wrong ontology, and the evidence for both halves of that was already in his two pages. He argued that the genome supplies the structure, and reported the counterexample. He argued that structure carries the qualitative answers, which is a theorem, and then listed the four perturbations his theories could not see. All four are binding phenomena. That list is the specification for what a structure has to contain, written by the person who had already built a solver for finding one, and a quarter of a century later it is still the part nobody has supplied.

References

80 entries. Every one is a full text in literature/, read rather than cited from an abstract. Two items that could not be obtained are logged in literature/NEEDED-structure-and-parameters.md with what was tried; no claim above rests on either.