Derive the respiratory-chain kcats from the data, and compare with the Xia et al. (2022) series

.gitignore

@@ -59,3 +59,18 @@ dist/
 # Solver
 gurobi.env
 gurobi.log
+
+# Built ecModels of the literature datasets
+literature/**/models/
+
+# Regenerated by literature/xia2022/scripts and comparison
+literature/xia2022/comparison/respiratory_limit/
+# Outputs of the Xia run that the comparison does not read
+literature/xia2022/results/modelSimulation/
+literature/xia2022/results/enzymeUsage/*.pdf
+literature/xia2022/results/enzymeUsage/systemMedians.tsv
+literature/xia2022/results/modelGeneration/average*
+literature/xia2022/results/modelGeneration/modifiedRibosomeSubunits_*
+literature/xia2022/results/modelGeneration/ribosomeSummary.tsv
+literature/xia2022/results/randomSampling/altExchangeFlux.tsv
+literature/xia2022/derived_kcats/respiratory_limit/

README.md

@@ -33,6 +33,7 @@ unavailable solver is an error, not a silent switch.
 | `tests/` | Test suite |
 | `results/` | Output tables and figures |
 | `legacy_matlab/` | The GECKO 2 MATLAB implementation and its results |
+| `literature/` | The analysis run on published datasets, and their comparison with `results/` |
 
 ## Experimental data
 
@@ -71,6 +72,33 @@ Each condition gets `models/ecModel_P_<cond>.yml` and, in `results/`, the
 enzymes whose measured abundance the model could not run on, the enzymes whose
 caps had to be released, and the full flux distribution.
 
+### From proteomics to enzyme caps
+
+Each measured abundance (the filtered mean + 1.96 SD of the replicates) caps its
+enzyme's usage. Enzyme complexes are handled as one unit throughout
+(`overflow.complexes`):
+
+1. **Complex levelling.** Absolute proteomics recovers the subunits of one complex
+   unevenly: subunits present 1:1 can come out tenfold apart while their replicates
+   agree within a few percent. Each complex therefore gets one abundance, the median
+   per copy over its measured subunits, written back to every subunit by its copies in
+   the model; unmeasured subunits of a complex are filled. Enzyme sets that share at
+   least two subunits and half the larger set (isozyme variants) are one complex.
+   Proteins in several complexes, or that also work alone, keep their own measurement:
+   cytochrome c, thioredoxin and LPD1, the E3 subunit of pyruvate dehydrogenase,
+   2-oxoglutarate dehydrogenase and the glycine cleavage system. `--no-complex-fix`
+   turns levelling off.
+2. **Minimum usage.** On a model limited only by the protein pool, a complex the model
+   cannot run at the dilution rate for less than its measured abundance is raised to
+   what it needs.
+3. **Reaching the dilution rate.** One LP finds the least relaxation that lets the
+   model grow at the dilution rate (`relax_to_growth`): it minimises the summed fold
+   increases of the measured caps, one fold per complex, with the protein pool allowed
+   to grow at 1000 times the price of a cap. The result does not depend on the order in
+   which limiting enzymes are found.
+4. **Measured rates.** With `--fit-rates`, the same LP under the measured exchange
+   rates raises the caps further by the least that reproduces them.
+
 Options change what the models are asked to account for:
 
 - `--fit-rates` holds CO2, oxygen and the byproducts within a tolerance of their
@@ -159,6 +187,26 @@ Writes `results/COMPARISON.md`: predicted exchange rates, median capacity usage
 per system and the sampled ATP and redox budget, this pipeline beside the MATLAB
 one, read from `legacy_matlab/results/`.
 
+## Respiratory-chain kcats
+
+With the ecModel's kcats (complex III 290/s and complex IV 1376/s, from BRENDA; ATP
+synthase 85/s), the measured respiratory chain cannot carry the measured respiration
+in several conditions, so the build raises its abundances to what the fluxes need and
+its capacity usage reads 100% by construction. `overflow.respiratory_kcats` derives
+these three kcats from the data instead. In every condition the apparent turnover of a
+complex is its sampled flux over the complex abundance the build uses (the filtered
+mean + 1.96 SD of its measured subunits, levelled to the median per copy). The
+highest apparent turnover across the conditions becomes the complex's kcat, which
+`overflow.build --respiratory-kcats` applies to every isozyme variant before building.
+No other kcat changes.
+
+Two assumptions follow: the condition that sets a kcat runs that complex at full
+capacity, so its saturation is built in and the derived kcat is a lower bound; and a
+dataset's kcats come from its own abundances, so a difference in absolute protein
+quantification between datasets cancels for the respiratory chain. The turnovers and
+kcats are in `results/modelGeneration/respiratoryTurnover.tsv` and
+`respiratoryKcats.tsv`.
+
 ## Running the whole analysis
 
 The analysis is run with `--fit-rates`, implemented in `overflow.build`
@@ -171,20 +219,32 @@ measured physiology. `--uptake-flex 1.08` gives CN4 the glucose it needs, and `-
 proton leak (see the build options).
 
 ```bash
-python -m overflow.build --fit-rates --rate-tolerance 0.08 --uptake-flex 1.08 --block r_2129
+python -m overflow.run_sampling --procs 12
+python -m overflow.respiratory_kcats
+python -m overflow.build --fit-rates --rate-tolerance 0.08 --uptake-flex 1.08 --block r_2129 \
+    --respiratory-kcats results/modelGeneration/respiratoryKcats.tsv
 python -m overflow.build_ribosome
 python -m overflow.analyze_usage
-python -m overflow.run_sampling --procs 12
 python -m overflow.compare
 ```
 
+Sampling uses the conventional GEM and no kcats, so it comes first: the kcats of the
+respiratory chain are derived from its fluxes (above).
+
 The tables and figures under `results/` are from a full run of this recipe. The
 models it builds (`models/ecModel_P_<cond>.yml`, about 67 MB) are not committed and
 rebuild in a few minutes.
 
 Each step reads the models the previous one wrote, so `--models-dir` and
 `--results-dir` keep a run self-contained.
 
+## Published datasets
+
+[`literature/xia2022/`](literature/xia2022/README.md) runs the analysis on nine
+conditions of the glucose-limited chemostat series of Xia et al. (2022) and
+compares them with `results/`, including each proteome's respiratory limit:
+the highest glucose uptake it can respire completely.
+
 ## Tests
 
 ```bash

literature/xia2022/README.md

@@ -0,0 +1,75 @@
+# Glucose-limited chemostat series (Xia et al. 2022)
+
+The analysis of this repository run on nine conditions of the glucose-limited chemostat series of Xia et al. (2022), *Nat Commun* 13:2819, [doi:10.1038/s41467-022-30513-2](https://doi.org/10.1038/s41467-022-30513-2), and compared with the results of this study in `results/`.
+
+| Condition | D [1/h] | Ptot [g/gDW] | Ethanol [mmol/gDW/h] |
+|---|---|---|---|
+| mu003 | 0.027 | 0.299 | – |
+| mu004 | 0.044 | 0.312 | – |
+| mu010 | 0.10225 | 0.371 | – |
+| mu015 | 0.152 | 0.419 | – |
+| mu021 | 0.21425 | 0.320 | – |
+| mu025 | 0.254 | 0.378 | – |
+| mu028 | 0.28425 | 0.391 | – |
+| mu033 | 0.33375 | 0.419 | 2.19 |
+| mu038 | 0.3785 | 0.391 | 5.83 |
+
+How the data were converted is in [`data/README.md`](data/README.md).
+
+## Running
+
+`scripts/run_step.py <step> [args]` points `overflow.config` at this folder's `data/`, condition names, `models/` and `results/`, then runs the step's `main`. For `build_ribosome` it supplies the core ribosome of `data/ribosome_core.tsv`. From the repository root, the recipe of this study:
+
+```bash
+X=literature/xia2022
+O="--models-dir $X/models --results-dir $X/results"
+python $X/scripts/run_step.py run_sampling --procs 12 --rate-tolerance 0.16 --results-dir $X/results
+python $X/scripts/run_step.py respiratory_kcats --results-dir $X/results
+python $X/scripts/run_step.py build --fit-rates --rate-tolerance 0.16 --uptake-flex 1.14 --block r_2129 --respiratory-kcats $X/results/modelGeneration/respiratoryKcats.tsv $O
+python $X/scripts/run_step.py build_ribosome $O
+python $X/scripts/run_step.py analyze_usage $O
+```
+
+The kcats of complex III, complex IV and ATP synthase are derived from this series, as they are from this study's data (see the repository README). `results/` keeps the outputs of this run that the comparison reads, and the sampled fluxes the kcats are derived from; the flux distributions of the built models, the ribosome tables and the pipeline's own usage figures are regenerated by the recipe and not committed, nor are the models (`models/`).
+
+### Departures from the recipe of this study
+
+- `--rate-tolerance 0.16` instead of 0.08, in the build and in sampling (0.05 in this study). Holding growth, glucose, CO2 and O2 together needs at least 15.4% at mu015, 12.7% at mu025 and 11.6% at mu010, even with protein unlimited. The measured O2 uptake of the series exceeds what complete oxidation of its measured glucose consumes by 7–46%, so the band has to absorb the difference. As in this study, the protein objective puts predicted glucose at the top of its band: 1.16× measured here, 1.08× there.
+- `--uptake-flex 1.14` instead of 1.08: at 1.08× measured glucose mu015 and mu025 cannot reach their dilution rate.
+- GAM is the model's (55.3), as in this study. Fitted on this series alone, with glucose fixed and O2 free, it would be about 35 (95% interval 22–48); this study's conditions without ethanol all grow at D = 0.1, so its data cannot fit one. One GAM keeps the ATP budgets of the two datasets comparable. At 55.3 the build's maintenance fit finds 0 for every condition but mu003 and mu004, and no non-negative value at mu010 and mu015; the rate band absorbs the difference.
+- The core ribosome is fixed rather than re-selected; see [`data/README.md`](data/README.md#core-ribosome).
+- `overflow.compare` is not run: it compares against the MATLAB results of this study's conditions only.
+- `--block r_2129` applies to the ecModels. Sampling uses the conventional GEM, where the proton leak stays open, as in this study.
+
+The series attributes all of Ptot to measured proteins, while in this study the measured proteins add up to 72–83% of it. For the respiratory chain this difference in calibration cancels, because each dataset's kcats are derived from its own abundances.
+
+## Comparison with this study
+
+`comparison/` puts both datasets side by side. It needs the ecModels of both runs, so run this study's recipe first.
+
+```bash
+for c in CN4 CN22 CN38 CN75 hGR;                                do for v in lifted measured; do python $X/comparison/respiratory_limit.py ThisStudy $c $v; done; done
+for c in mu003 mu004 mu010 mu015 mu021 mu025 mu028 mu033 mu038; do for v in lifted measured; do python $X/comparison/respiratory_limit.py Xia $c $v; done; done
+python $X/comparison/figures.py
+```
+
+`respiratory_limit.py` finds a condition's respiratory limit: the highest glucose uptake its measured proteome can respire completely at the condition's growth rate. Above it, the model has to secrete overflow carbon. The ecModel is held at the dilution rate with maintenance ATP free, so surplus glucose can be respired until a measured enzyme cap stops it. Caps the pipeline raised (minimum usage, relaxation to the dilution rate, the rate fit) show that the measurement understates what the condition needs, carry no capacity information and are lifted, for the whole complex. Unlisted exits are closed except glycine export, which levelled models need at some growth rates. As a conservative bound, the `measured` variant puts the caps the rate fit raised back to their measured abundance instead of lifting them (`lifted`). The limit is the largest uptake with no more than 0.05 mmol C/gDW/h of forced ethanol, acetate, glycerol and formate. Each run writes its curve and limit to `comparison/respiratory_limit/` (not committed); `figures.py` collects them into Table 5 and figure 5. A condition is predicted to overflow when its measured glucose uptake exceeds its limit; the onset of overflow itself is the measured physiology. As the respiratory-chain kcats are the highest apparent turnover within each dataset, the prediction is a consistency check rather than an independent test: the conditions that set the kcats run their respiratory chain at capacity.
+
+| Output (`comparison/`, figures as PDF) | Contents |
+|---|---|
+| `figures/fig1_capacity_usage` | Capacity usage of glycolysis, TCA cycle, ETC and ribosome, both datasets (`figS1`: the supplementary systems). |
+| `figures/fig2_usage_vs_uptake` | Capacity usage per system (enzyme mass used over the measured capacity of its enzymes) against measured glucose uptake. |
+| `figures/fig3_proteome_reserve` | Enzyme mass used in the simulation against the measured capacity, per system. |
+| `figures/fig4_atp_budget` | Sampled ATP supply and demand, and the respiratory share of ATP. |
+| `figures/fig5_respiratory_limit` | Least overflow the proteome allows against glucose uptake, with the respiratory limits (raised caps lifted) and the measured byproduct carbon. |
+| `figures/fig6_apparent_turnover` | Apparent turnover of complex III, complex IV and ATP synthase relative to the highest in each dataset: sampled flux over measured complex, which needs no kcat. |
+| `tables/tables.md` | Tables 1–7 (also as `table<n>.tsv`): conditions and fit, capacity usage per system, enzyme mass used and measured capacity, ATP and redox budget, respiratory limit, the derived respiratory-chain kcats, and physiology with marker proteins. |
+
+### CN22 and μ0.28: what triggers overflow
+
+Both are the last condition of their series before ethanol, but they reach that point from opposite directions (Tables 3, 5 and 7).
+
+- **Ethanol: glucose uptake outgrows the respiratory capacity.** In this study glucose uptake rises with C/N while the TCA and ETC enzymes shrink with Ptot (TCA capacity 9.4, 7.2, 4.7 and 3.3 mg/gDW from CN4 to CN75). CN22 takes up 94% of a limit set by succinate dehydrogenase; CN38 is past its limit. In the Xia series the respiratory proteome grows with growth rate and peaks at μ0.28 (SDH1, CIT1, TCA capacity), which leaves it 27% below its limit. From μ0.33 the limit stops rising (5.06 at μ0.33 and μ0.38, set by complex IV) while glucose uptake keeps rising, and the surplus leaves as ethanol.
+- **Formate follows nitrogen limitation, not respiratory saturation.** It is secreted at CN22, CN38 and CN75 and rises with C/N (0.05, 0.34 and 0.48 mmol/gDW/h), while hGR secretes the most ethanol of this study and no formate. A C1 imbalance would explain it: C1 units come from serine, made from glycolytic carbon that is in excess, while their main sink, the nitrogen-rich purines, shrinks with nitrogen (ribosome capacity 43 to 22 mg/gDW from CN4 to CN75). The glycine cleavage system, which releases ammonium from glycine and leaves a C1 unit, is 1.6–2.7 times higher in the conditions that secrete formate than at CN4 and hGR, and FDH1, which would oxidise formate to CO2, is not detected in this study. In the carbon-limited series the glycine cleavage system stays level and FDH1 is detected.
+
+Limits of this reading: Xia et al. do not report formate, so a C1 overflow at μ0.28 cannot be excluded; the glycine cleavage system does not scale with the amount of formate; and the model secretes formate only where the measured rate imposes it, so it cannot test the mechanism. The proteomics of the two datasets are calibrated differently, and Xia et al. do not detect several membrane proteins (GAP1, the hexose transporters), so protein levels are compared within each dataset.