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

.gitignore

@@ -59,3 +59,10 @@ 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/
+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 three
+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,59 @@
+# Glucose-limited chemostat series (Xia et al. 2022)
+
+The analysis of this repository run on three 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] |
+|---|---|---|---|
+| mu010 | 0.10225 | 0.371 | – |
+| mu028 | 0.28425 | 0.391 | – |
+| mu033 | 0.33375 | 0.419 | 2.19 |
+
+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 --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.12 --uptake-flex 1.08 --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/` holds the tables and figures of this run. The models (`models/`) are not committed.
+
+### Departures from the recipe of this study
+
+- `--rate-tolerance 0.12` instead of 0.08. At 0.08 mu010 has no solution: holding glucose, CO2 and O2 together needs at least 11.6%, even with protein unlimited. As in this study, the protein objective puts predicted glucose at the top of its band: 1.12× measured here, 1.08× there.
+- GAM is the model's (55.3), as in this study. The series' respiratory rates imply a lower GAM (about 32 with NGAM 0.64), which was not used.
+- 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 built measured; do python $X/comparison/respiratory_limit.py ThisStudy $c $v; done; done
+for c in mu010 mu028 mu033;      do for v in built 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 so that the model could grow (minimum usage, relaxation to the dilution rate) 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. Caps the rate fit raised are kept as built (`built`), or put back to the measured abundance (`measured`). 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.
+
+| 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 (caps as built) 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–6 (also as `table<n>.tsv`): conditions and fit, capacity usage per system, enzyme mass used and measured capacity, ATP and redox budget, respiratory limit, and the derived respiratory-chain kcats. |